% Nutrient Diffusion Simulations Imperial iGEM 2014
% For more information, please refer to our wiki page:
% http://2014.igem.org/Team:Imperial/Modelling

numx = 50;   % x_coordinate
numy = 50;   % y_coordinate
numt = 2000;  % number of time steps to be iterated over 

%parameters
%cell growth
initial_acetate_concentration=1;
acetate_release_time_trigger=5;
top_acetate_feed=1;
cell_growth_time_trigger=30;    %cell growth time interval AND nutrient consumption
acetate_consumption_amount=0.1;    %associated with cell growth time trigger
oxygen_consumption_amount=0.1;
cell_growth_rate=1;   %integers only

%cellulose production
initial_glucose_concentration=1;
glucose_release_time_trigger=5;
top_glucose_feed=1;
cell_density_threshold=5; %cell cellulose production is triggered when cell density is greater than the threshold
glucose_consumption_amount=0.1;     %associate with cellulose production time trigger
cellulose_production_time_trigger=30;
cellulose_production_rate=1; %not necessarily integer

DA=1;    %acetate diffusion coefficient
DO=1;    %oxygen diffusion coefficient
DG=1;    %glucose diffusion coefficient
A_dependent_coefficient=1;  %cellulose bulk effect
O_dependent_coefficient=1;
G_dependent_coefficient=1;


dx = 1/(numx - 1);
dy = 1/(numy - 1);
dt = 0.00005;

x = 0:dx:1;   %vector of x values, to be used for plotting
y = 0:dy:1;   %vector of x values, to be used for plotting

% acetate profile-vineger
A = zeros(numy, numx, numt);    %create a matrix for carbon source
A(2:(numy-1), 2:(numx-1), 1) = initial_acetate_concentration;     %initialise carbon source distribution

% glucose profile-sugar
G = zeros(numy, numx, numt);    %create a matrix for carbon source
G(2:(numy-1), 2:(numx-1), 1) = initial_glucose_concentration;     %initialise glucose distribution

% oxygen profile initial and boundary conditions
O = zeros(numy, numx, numt);    %create a matrix for oxygen
O(1:numy, 1:numx, 1) = 0;  %O(x, y, 0)=0 initial condition
O(numy, 1:numx,1:numt) = 1;    %O(x, numy, t)=1 oxygen: top boundary = 1

% cell profile
D = zeros(numy, numx);

% cellulose production profile
X = zeros(numy, numx);

H = 1/(sqrt(2));

for k=1:(numt-1)
  
    check=check_trigger_time(k, acetate_release_time_trigger);    %check trigger for carbon source release (from the top), time interval=50
    if check==1     %carbon source release (from top)
        A(numy-1, 2:(numx-1), k)=top_acetate_feed;
    end
    
    check=check_trigger_time(k, glucose_release_time_trigger);    %check trigger for carbon source release (from the top), time interval=50
    if check==1     %glucose source release (from top)
        G(numy-1, 2:(numx-1), k)=top_glucose_feed;
    end
    
    check=check_trigger_time(k, cell_growth_time_trigger);    %check trigger for oxygen AND carbon source consumption (by cells) AND cell growth
    if check==1 
        for j=2:(numy-1)
            for i=2:(numx-1) 
                if ((D(j,i)==0)&&(O(j, i, k)>=oxygen_consumption_amount)&&(A(j, i, k)>=acetate_consumption_amount))
                    O(j, i, k) = O(j, i, k)-oxygen_consumption_amount;
                    A(j, i, k) = A(j, i, k)-acetate_consumption_amount;
                    D(j, i) = 1;
                end
                if ((D(j,i)>0)&&(O(j, i, k)>=D(j,i)*oxygen_consumption_amount)&&(A(j, i, k)>=D(j,i)*acetate_consumption_amount))
                   O(j, i, k) = O(j, i, k)-D(j,i)*oxygen_consumption_amount;
                   A(j, i, k) = A(j, i, k)-D(j,i)*acetate_consumption_amount;
                   D(j, i) = D(j, i)+cell_growth_rate; %linear growth rate
                end
            end
        end
    end
    
    check=check_trigger_time(k, cellulose_production_time_trigger);     %check trigger for glucose consumption AND cellulose production
    if check==1
        for j=2:(numy-1)
            for i=2:(numx-1)
                if ((D(j,i)>cell_density_threshold)&&(G(j, i, k)>=D(i, j)*glucose_consumption_amount))
                    G (j, i, k) = G (j, i, k) - D(i, j)*glucose_consumption_amount;
                    X (j, i) = X(j, i)+cellulose_production_rate; %linear cellulose production
                end
            end
        end
    end
    
   for j=2:(numy-1)
        for i=2:(numx-1)
            %diffusion coefficient variation
            if X (j, i)>0
                DA=DA/(X(j, i)*A_dependent_coefficient);
                DO=DO/(X(j, i)*O_dependent_coefficient);
                DG=DG/(X(j, i)*G_dependent_coefficient);
            end
            %carbon diffusion
            A(j, i, k+1) = DA*(A(j, i, k) + (dt/(dy^2))*(A(j+1, i, k) - 2*A(j, i, k)+A(j-1, i, k))+ (dt/(dx^2))*(A(j, i+1,k)-2*A(j, i, k)+A(j, i-1, k))+H*((dt/(dy^2))+(dt/(dx^2)))*((A(j+1, i+1, k)-2*A(j, i, k)) +(A(j+1, i-1, k))+(A(j-1, i-1, k)-2*A(j, i, k))+(A(j-1, i+1, k)))) ;
            %oxygen diffusion
            O(j, i, k+1) = DO*(O(j, i, k) + (dt/(dy^2))*(O(j+1, i, k) - 2*O(j, i, k)+O(j-1, i, k))+ (dt/(dx^2))*(O(j, i+1,k)-2*O(j, i, k)+O(j, i-1, k))+H*((dt/(dy^2))+(dt/(dx^2)))*((O(j+1, i+1, k)-2*O(j, i, k)) +(O(j+1, i-1, k))+(O(j-1, i-1, k)-2*O(j, i, k))+(O(j-1, i+1, k)))) ;
            %glucose diffusion
            G(j, i, k+1) = DG*(G(j, i, k) + (dt/(dy^2))*(G(j+1, i, k) - 2*G(j, i, k)+G(j-1, i, k))+ (dt/(dx^2))*(G(j, i+1,k)-2*G(j, i, k)+G(j, i-1, k))+H*((dt/(dy^2))+(dt/(dx^2)))*((G(j+1, i+1, k)-2*G(j, i, k)) +(G(j+1, i-1, k))+(G(j-1, i-1, k)-2*G(j, i, k))+(G(j-1, i+1, k)))) ;
        end
   end
    
   for i= 1:numx
       A(1, i, k+1) = A(2, i, k);          %acetate: insulated top and bottom boundary
       A(numy, i, k+1)=A(numy-1, i, k);
   end
   
   for j= 2:(numy-1)                         %acetate: insulated side walls
       A(j, 1, k+1)= A(j, 2, k);
       A(j, numx, k+1)=A(j, numx-1, k);
   end
   
   for i= 1:numx
       G(1, i, k+1) = G(2, i, k);          %glucose: insulated top and bottom boundary
       G(numy, i, k+1)=G(numy-1, i, k);
   end
   
   for j= 2:(numy-1)                         %glucose: insulated side walls
       G(j, 1, k+1)= G(j, 2, k);
       G(j, numx, k+1)=G(j, numx-1, k);
   end
   
   for i= 1:numx
       O(1, i, k+1) = O(2, i, k);          %oxygen: insulated bottom boundary
   end
   
   for j= 2:(numy-1)                         %oxygen: insulated side walls
       O(j, 1, k+1)= O(j, 2, k);
       O(j, numx, k+1)=O(j, numx-1, k);
   end
        
end

for k=1:1:numt  %display
    surf(G(:,:,k)); %A:acetate, O:oxygen, G: glucose
    caxis([0 1]);
    view(2);
    drawnow;
    pause(1/20);
end

% End of code
% For further enquiries please contact dk3712 AT imperial DOT ac DOT uk
