% Imperial iGEM 2014 Induction model
% This model was designed to investigate the effect of induction on overall
% bacterial cellulose (BC) yield
% This model was used to assist the wetlab by proposing means for improving
% BC yield within an experiment time frame
% Please refer to our wiki page
% (http://2014.igem.org/Team:Imperial/Induction_Model) for further
% information

clear all; clc; close all;

dt=0.05;
tend=5000; % total time setting
time_scaling_factor = 5000/50;
induction_index = 50; % induction time setting
induction_interval = 1;

Substrate = zeros([1 tend]);
Population_plus = zeros([1 tend]);
Population_minus = zeros([1 tend]);
Cellulose = zeros ([1 tend]);
time_scale (1:tend) = 0.01:0.01:50;
diffusion_coefficient = zeros ([1 tend]); % simulate a diffustion coefficient
induction_time = zeros([4 induction_index]);
first_time = 1;
frame_number = 1;

% parameters
Vp_original = 0.3; % cell population
Vp_induction = 0.04;
Kp = 6;
Vc_original = 0.05; % cellulose production
Vc_induction = 0.1;
Kc = 10;
Vd = 0.05; % cell death rate
Ko2 = 5; % oxygen factor
Vm = 0.08; % mutation rate from +cells to -cells
Kdc = 100; % variation in diffusion coefficient
Oxygen = 5E1;
oxygen_factor = (Ko2+Oxygen)/Ko2;  % oxygen sensing upregulate cellulose production

% initialisation
Substrate(1) = 1E4;
Oxygen (1:tend) = 5E1; % oxygen level, input from pellicle model
Population_plus(1) = 5E-3;
Population_minus(1) = 0;
Cellulose(1) = 0;
diffusion_coefficient(1) = 1;

% code
for index_it = 1:induction_index
    induction_time(1, index_it)=index_it*induction_interval;
end

for index_it = 1:induction_index;
    Vp = Vp_original;
    Vc = Vc_original;
    
    for time = 1:(tend-1)
        % kinetics
        if Substrate(time)>0
            substrate_level_cells = Vp*Substrate(time)/(Kp+Substrate(time));
            substrate_level_cellulose = Vc*Substrate(time)/(Kc+Substrate(time));
            Population_plus(time+1) = Population_plus(time)+dt*Population_plus(time)*(((1-Vm)*substrate_level_cells));
            Population_minus(time+1) = Population_minus(time)+dt*(Population_plus(time)*(Vm*substrate_level_cells));
            Cellulose(time+1) = Cellulose(time)+dt*Population_plus(time)*oxygen_factor*substrate_level_cellulose;
            Substrate(time+1) = Substrate(time)+dt*Population_plus(time)*(-substrate_level_cells-(oxygen_factor*substrate_level_cellulose));
        else
            Population_plus(time+1) = Population_plus(time)-dt*Vd*Population_plus(time);
            Population_minus(time+1) = Population_minus(time)-dt*Vd*Population_minus(time);
            Cellulose(time+1) = Cellulose(time);
            Substrate(time+1) = Substrate(time);
        end
        
        % induction
        if time == 100*induction_time(1, index_it)
            Vp = Vp_induction;
            Vc = Vc_induction;
        end
    end
    
    induction_time (2, index_it) = max(Cellulose(1:tend)); % max cellulose production for each induction time
    induction_time (3, index_it) = max(Population_plus(1:tend)); % max population achieved for each induction time
    
    first_max_cellulose = 1;
    for time = 1: tend
        if Cellulose (time) == induction_time (2, index_it) && first_max_cellulose == 1 % time to achieve max cellulose production for each induction time
            induction_time (4, index_it) = time/time_scaling_factor;
            first_max_cellulose = 0;
        end
    end
    
    axis([0 50 0 10000]);
    plot (time_scale, Substrate, 'k');
    hold on
    axis([0 50 0 10000]);
    plot (time_scale, Cellulose, 'g');
    axis([0 50 0 10000]);
    plot (time_scale, Population_plus, 'b');
    axis([0 50 0 10000]);
    plot (time_scale, Population_minus, 'c');
    legend ('Substrate', 'BC Production','+Cell Population', '-Cell Population');
    xlabel('Time (hour)');
    hold off
    
    if first_time == 1
        axis(axis);
        set(gca,'nextplot','replacechildren','visible','on');
        f=getframe; [im,map]=rgb2ind(f.cdata,256,'nodither');
        im(1,1,1,induction_index)=0;
        first_time = 0;
    else
        axis(axis);
        f = getframe; im(:,:,1,frame_number)=rgb2ind(f.cdata,map,'nodither');
    end
    
    frame_number = frame_number+1;
    
end
imwrite(im,map,'induction_model_optimal_induction_time.gif','DelayTime',0,'LoopCount',inf);

figure (2); % plot for induction time against cellulose production
max_cellulose_production = max(induction_time(2, 1:induction_index));
for index_it = 1:induction_index; % point plot
    if induction_time (2, index_it) == max_cellulose_production
        max_cellulose_production_induction_time = induction_time (1, index_it);
    end
end
max_cell_population = max(induction_time(3, 1:induction_index));
plot ( induction_time(1, 1:induction_index), induction_time(2, 1:induction_index));
hold on;
plot (max_cellulose_production_induction_time, max_cellulose_production, 'kx');
xlabel ('Induction Time (hour)');
ylabel ('Maximum Bacterial Cellulose Yield (arbitrary unit)');
title ('Induction Time Against Maximum Bacterial Cellulose Yield');

figure (3); % plot for induction time against time to achieve maximum cellulose production
plot (induction_time(1, 1:induction_index), induction_time (4, 1:induction_index));
xlabel ('Induction Time (hour)');
ylabel ('Time to Achieve Maximum Bacterial Cellulose Yield (hour)');
title ('Induction Time Against Time to Achieve Maximum Bacterial Cellulose Yield');

figure (4); % plot for time to achieve max cellulose production against max cellulose production
scatter (induction_time (4, 1:induction_index), induction_time (2, 1:induction_index), 'x');
for index_it = 1:induction_index
    label = num2str (induction_time (1, index_it));
    text (induction_time (4, index_it), induction_time (2, index_it), label, 'horizontal','left', 'vertical','bottom');
end
xlabel ('Time to Achieve Maximum Bacterial Cellulose Yield (hour)');
ylabel ('Maximum Bacterial Cellulose Yield (arbitrary unit)');
title ('Induction Time Against Time to Achieve Maximum Bacterial Cellulose Yield');
legend ('induction time (hour)','Location','SouthWest');

% end of the code
% For enquiries please contact dk3712 AT imperial DOT ac DOT uk