%%

%parameters: 
tic
number_of_cells = 10000;
velocity = 14;
cells = zeros(6,number_of_cells);
strain_1_radius = 1000;
strain_2_radius = 2000;
center = [3000,3000];
step = 0.1;
dimension = 6000;
grids_dimension = 15;
grids_scale = dimension/ grids_dimension;
growth_rate = 0.3;
death_rate = 0.3;

ten_minutes = 6000;
one_hour = 36000;

Num_steps = one_hour*16;
beta = 5e-7; %production rate of HSL
omega = 5e-4; %degradation rate of HSL
lambda = 8e-5; %production rate of Mscarlet (depending on HSL)
w = 6e-5; %degradation rate of Mscarlet
alpha = 2e-4; %diffusion rate of HSL

IPTG_level = 1000;
normalized_IPTG = sigmf((4/19000)*(10000-IPTG_level),[2 2]);


cells = initialize_cells(cells,dimension);
for j = 1:Num_steps
    for i = 1:length(cells)
        A = cells(1,i);
        X = cells(3,i);
        Y = cells(4,i);
        direction = cells(2,i);
        temp = (X- center(1))^2 + (Y - center(2))^2;
        if temp <= strain_1_radius^2
            I = 1;
        else
            I = 0;
        end
        stick = cells(5,i);
        new_A = A + step*(10e-3 + I*10e-1 - 10e-3*A);
        cells(1,i) = new_A;
        
        if stick == 0
            new_X = X + step*14*cos(direction);
            new_Y = Y + step*14*sin(direction);
            if new_X > dimension || new_X < 0 || new_Y > dimension || new_Y < 0
                cells(3,i) = X;
                cells(4,i) = Y;
                cells(2,i) = rand*2*pi;
            else
                cells(3,i) = new_X;
                cells(4,i) = new_Y;
                if rand < 0.1
                    cells(2,i) = rand*2*pi;
                end
            end
            p_adsorb = step*10e-5;
            if rand < p_adsorb
                cells(5,i) = 1;
            end
        else
            p_desorb = step*10e-3/new_A;
            if rand < p_desorb
                cells(5,i) = 0;
            end
        end 
    end
    if (rand < 0.0001)
        cells = grow_and_die(cells, 0.3, 0.3);
    end
end
strain_1 = cells;
for i = 1:number_of_cells
    cells(1, i) = 1;
    cells(2, i) = rand*2*pi;
    cells(3,i) = randi(dimension);
    cells(4,i) = randi(dimension);
    while (cells(3,i)- center(1))^2 + (cells(4,i) - center(2))^2 <= strain_1_radius^2
        cells(3,i) = randi(dimension);
        cells(4,i) = randi(dimension);   
    end
    cells(5,i) = generate_zero_one(rand);
end
for j = 1:Num_steps
    for i = 1:number_of_cells
        A = cells(1,i);
        X = cells(3,i);
        Y = cells(4,i);
        direction = cells(2,i);
        temp = (X- center(1))^2 + (Y - center(2))^2;
        if temp <= strain_2_radius^2
            I = 1;
        else
            I = 0;
        end
        stick = cells(5,i);
        new_A = A + step*(10e-3 + I*10e-1 - 10e-3*A);
        cells(1,i) = new_A;

        if stick == 0
            new_X = X + step*14*cos(direction);
            new_Y = Y + step*14*sin(direction);
   
            if new_X > dimension || new_X < 0 || new_Y > dimension || new_Y < 0 || (new_X- center(1))^2 + (new_Y - center(2))^2 <= strain_1_radius^2
                cells(3,i) = X;
                cells(4,i) = Y;
                cells(2,i) = rand*2*pi;
                
            else
                cells(3,i) = new_X;
                cells(4,i) = new_Y;
                if rand < 0.1
                    cells(2,i) = rand*2*pi;
                end
            end
           
            p_adsorb = step*10e-5;
            if rand < p_adsorb
                cells(5,i) = 1;
            end
            
        else
            p_desorb = step*10e-3/new_A;
            if rand < p_desorb
                cells(5,i) = 0;
            end
        end 
    end
end 

strain_2 = cells;
index1 = strain_1(5,:) == 1;
temp = [];
for i = 1:length(index1)
    if ((strain_1(3,i)- center(1))^2 + (strain_1(4,i) - center(2))^2 <= strain_1_radius^2)
        temp = [temp,i];
    end
end
index1 = temp;
index2 = strain_2(5,:) == 1;

%%
figure;
plot(strain_1(3,index1), strain_1(4,index1),'.', 'Color', [0, 0.5, 0]	)

hold on
plot(strain_2(3,index2), strain_2(4,index2),'.', 'Color', [0.90, 0.70, 0.30]	)
hold off
legend("sender strain", 'receiver strain')
%% 
HSL_grids = zeros(grids_dimension);
bacteria_A = strain_1(:,index1);
bacteria_B = strain_2(:,index2);


%%

for k = 1:Num_steps
    %in this for loop, for each step
    %first, update each grid's HSL diffusion rate
    for j = 1:grids_dimension
        for i = 1:grids_dimension
            delta_HSL = alpha * sum_around(i,j,HSL_grids) - omega * HSL_grids(i,j);
            HSL_grids(i,j) = HSL_grids(i,j) + delta_HSL * step;
        end
    end
    
    %second, update bacteria A's HSL production for the HSL grids
    for i = 1:length(bacteria_A)
        delta_HSL = beta * normalized_IPTG;
        x_position = ceil(bacteria_A(3,i)/grids_scale);
        y_position = ceil(bacteria_A(4,i)/grids_scale);
        HSL_grids(x_position, y_position) = HSL_grids(x_position, y_position) + delta_HSL * step;
    end 

    %third, update bacteria B's Mscarlet amount
    for i = 1:length(bacteria_B)
        x_position = ceil(bacteria_B(3,i)/grids_scale);
        y_position = ceil(bacteria_B(4,i)/grids_scale);
        
        M = bacteria_B(6,i);
        delta_M = lambda * HSL_grids(x_position,y_position) - w * M;
        bacteria_B(6,i) = bacteria_B(6,i) + delta_M * step;
    end 
end 

x = bacteria_B(3,:);
y = bacteria_B(4,:);
z = bacteria_B(6,:);

Mscarlet_grids = zeros(grids_dimension);
for i = 1:length(bacteria_B(3,:))
    x_position = ceil(bacteria_B(3,i)/grids_scale);
    y_position = ceil(bacteria_B(4,i)/grids_scale);
    Mscarlet_grids(x_position, y_position) = Mscarlet_grids(x_position, y_position) + bacteria_B(6,i);
end


%%

figure
clf
subplot(2,1,1);
heatmap(HSL_grids);
title("HSL diffusion");

subplot(2,1,2);
heatmap(Mscarlet_grids);
title('Mscarlet diffusion');
toc 

figure
clf
heatmap(HSL_grids, 'CellLabelFormat', '%0.2f');
str = sprintf("HSL concentration level after diffused for %d hours", Num_steps/60/60/10);
title(str);

figure
clf
heatmap(Mscarlet_grids, 'CellLabelFormat', '%0.2f');
str = sprintf("Normalized fluorescence intensity after %d hours", Num_steps/60/60/10);
title(str);


function y = generate_zero_one(x)
if x < 0.5
    y = 0;
else
    y = 1;
end

end 
function z = sum_around(i,j,HSL_grids)
sum = 0;
for x = -1:1
    for y = -1:1
        if i+x >= 1 && i+x <= length(HSL_grids) && j+y >=1 && j+y <= length(HSL_grids)
            sum = sum + HSL_grids(i+x, j+y);
        end 
    end 
end 
z = sum - HSL_grids(i,j);
end 
function y = initialize_cells(cells, dimension)
for i = 1:length(cells)
    cells(1, i) = 1;
    cells(2, i) = rand*2*pi;
    cells(3,i) = randi(dimension);
    cells(4,i) = randi(dimension);
    cells(5,i) = generate_zero_one(rand);
end 
y = cells;
end

function y = grow_and_die(cells, growth_rate, death_rate)
    number_of_cells_divided = floor(growth_rate* length(cells));
    cells_divided = randperm(length(cells), number_of_cells_divided);
    cells_divided_columns = cells(:,cells_divided);
    for k = 1: length(cells_divided_columns)
        cells_divided_columns(5,k) = generate_zero_one(rand);
    end
    
    number_of_cells_died = length(cells_divided);
    cells_died = randperm(length(cells), number_of_cells_died);

    cells(:,cells_died) = [];
    length(cells)
    cells = cat(2,cells, cells_divided_columns);
    y = cells;
end
