Figure 8a:
Distribution of the GFP at t = 3.5 hr.
Code for Figure 8a
Text of the GNU GPL.
main.m
clear all
format compact
% set seeds for reproducible runs
rand ("seed", -1)
randn("seed", -1)
% Initiation
tol=0.01; % maximum allowed change (0.01 = 1%) in Langevin species per time step;
max_step=10; % maximum allowed time step, when no stochastic reaction occurs;
checksize = 100; % species level has to be higher than checksize to switch to Langevin implementation
runs=1000; % number of simulation runs
%% runs=250; % smaller number of run for debugging
tfinal=14400; % simulation time 14400 secs = 4 hours
% Kinetic paramters:
factor=1.29e9; % conversion factor for the parameters found in the literature
% (L/mol/sec --> n/sec): factor = 4/3 * pi * (16e-6/2)^3 * 6.022e23; volume*Avogadro
% (The average cell is assumed to be spheric with a radius of 16 micrometers)
k2=0.0755; % protein translation rate per mRNA
k3=(4e7+7.9e7)/2/factor; % transcription and replication initiation, average of 2 Qbeta (Werner 1991)
k1=[1 .75 .5625 .422 .422 .0633]*k3; % transcritption and replication initiation rate
% Initial condition:
y=[1 zeros(1,2523) 1258 466 1826 1205 0 50 2.962e11]; % initial condition of all states
% max_events sets the maximum number of stored events for delayed reactions
% the smaller these arrays are, the faster you can look up the time stamps of the stored delayed reactions
max_events1 = 5000;
max_events2 = 4500;
max_events3 = 4000;
max_events4 = 3500;
max_events5 = 3000;
max_events6 = 600;
max_events7 = 8000;
max_events8 = 2000;
max_events9 = 3000;
max_events10 = 3000;
max_events11 = 3000;
max_events12 = 3000;
max_events13 = 4000;
max_events14 = 700;
max_events15 = 9000;
max_events16 = 4000;
max_events17 = 1000;
% Delays
nuc_length=[1333;822;838;1672;720;6380]; % nucleotide length
tau=((cumsum(nuc_length)./3.7))+cumsum(240*ones(6,1)); % delays for transcription (time to transcribe the gene + stop time (240 secs) at each gene junction (4 mins))
tau(7:8)=sum(nuc_length)/3.7; % delays for replication (time to replicate the genome, NO stop time at gene junctions)
% Calculations for Poisson distribution
p=zeros(101,1); % will later be assigned with the probability of being in states 0 to 100
facts = zeros(101,1); % factorials that will later be used to calculate the probability distribution
for j=1:101
facts(j) = gamma(j);
end
% Storing time and state vectors
% It would be impossible to store all runs and states
% We select a number of points in time and the states we want to store in order to reduce the amount of memory needed
time_vec=0:1:tfinal; % time vector, we store all states every second of the run
time_vec=time_vec';
for c=1:18
eval(['x',num2str(c),'=zeros(tfinal+1,runs);'])
end
% Start of simulation:
for u=1:runs % Loop for all simulation runs
% First part of simulation
% No fast switching N species, therefore, all states are modeled and only 15000 iterations or about 2.5 hours, whatever comes first
maxi=15000; % maximum number of iterations in the first part of the simulation
tend=10000; % maximum simulated time until switch into quasi-steady state implementation of the fast switching species
% v1: stoichiometric matrix of the non-delayed reactions
% v2: stoichiometric matrix of the delayed reactions
v1=zeros(2531,2531);
v2=zeros(15,2531);
v1(1,2530)=-1;
v1(2:1259,2525)=-1;
for j=2:1259
v1(j,j-1)=-1;v1(j,j)=1;
end
v1(1260,2530)=-1;
v1(1261:2518,2525)=-1;
for j=1261:2518
v1(j,j-1)=-1;v1(j,j)=1;
end
for j=2519:2524
v1(j,2530)=-1;
end
for j=2525:2530
v1(j,j)=1;
end
v1(2525,2531)=-444;
v1(2526,2531)=-274;
v1(2527,2531)=-279;
v1(2528,2531)=-557;
v1(2529,2531)=-240;
v1(2530,2531)=-2127;
for j=2519:2524
v2(j-2518,j)=1;
end
v2(7,1)=1;
v2(8,1260)=1;
v2(9:15,2530)=1;
% event times et of all delayed reactions are initially set to inf
et1 = Inf * ones (max_events1, 1);
et2 = Inf * ones (max_events2, 1);
et3 = Inf * ones (max_events3, 1);
et4 = Inf * ones (max_events4, 1);
et5 = Inf * ones (max_events5, 1);
et6 = Inf * ones (max_events6, 1);
et7 = Inf * ones (max_events7, 1);
et8 = Inf * ones (max_events8, 1);
et9 = Inf * ones (max_events9, 1);
et10 = Inf * ones (max_events10, 1);
et11 = Inf * ones (max_events11, 1);
et12 = Inf * ones (max_events12, 1);
et13 = Inf * ones (max_events13, 1);
et14 = Inf * ones (max_events14, 1);
et15 = Inf * ones (max_events15, 1);
et16 = Inf * ones (max_events16, 1);
et17 = Inf * ones (max_events17, 1);
% Initialization (has to be done for each simulation);
T = zeros(maxi,1); % time vector
x = zeros(maxi,2531); % state vector for all states
g = zeros(maxi,2); % g(:,1) = sum of all (-)RNA genomes but the naked strand, g(:,2) = sum of all (+)RNA genomes but the naked strand
x_tau = zeros(7,2531); % states of the past, only needed for delayed reactions that will be implemented with Langevin equations (continuous part)
genomes = zeros(maxi,1); % all (-)RNA genomes that are not fully encapsidated (still produce mRNAs and not used for replication)
alpha = zeros(maxi,1); % mean of fast switching species N (Poisson distribution)
T(1) = 0; % initial time
x(1,:) = y; % initial state
current_time = 0; % current simulation time (needed for storage of delayed reactions)
a = zeros(2531,1); % all reaction rates
rd = zeros(2531,1); % all continuous reaction rates
r_con = zeros(2531,1); % all continuous reaction rates that do not take part in the stochastic time step / reaction calculation
sd = zeros(2531,1); % all Langevin reaction rates
s_con = zeros(2531,1); % all Langevin reaction rates that do not take part in the stochastic time step / reaction calculation
%checks for fast and slow reactions Default to slow.
flag = ones(2531,1); % flags for non-delayed reactions, 1 for stochastic, 0 for Langevin implementation
flagp= ones(2531,1); % flags for delayed reactions, 1 for stochastic, 0 for Langevin implementation
pos=ones(17,1); % index of all stored delayed reactions (some reactions have more than one product at different times)
% Start of first part of the simulation until quasi-steady state of N is reached
for i=1:maxi-1
genomes(i)=sum(x(i,1:1258)); % all (-)RNA genomes that are not fully encapsidated (still produce mRNAs and not used for replication)
% Calculation of all reaction rates:
a(1) = 50*k1(1)*x(i,2518)*x(i,2530);
a(2:1259) = k3*x(i,1:1258).*x(i,2525);
a(1260) = k1(1)*x(i,1259)*x(i,2530);
a(1261:2518) = k3*x(i,1260:2517).*x(i,2525);
a(2519:2524) = k1(1:6).*genomes(i).*x(i,2530);
a(2525:2530) = k2.*x(i,2519:2524);
% Flags 2525 to 2530: (translation)
for j=2525:2530
if x(i,j) there are continuous reactions
if flag_mult == 0 % if there are continuous reactions
% Calculate the maximum step size:
% Calculation of all continuous reaction rates
r_con=zeros(2531,1);
for j=2525:2531
rd(j)=a(j)*(1-flag(j));
end
r_con(2525)=rd(2525);
r_con(2526)=rd(2526);
r_con(2527)=rd(2527);
r_con(2528)=rd(2528);
r_con(2529)=rd(2529);
r_con(2530)=rd(2530);
r_con(2531)=-rd(2525)*444-rd(2526)*274-rd(2527)*279-rd(2528)*557-rd(2529)*240-rd(2530)*2127; % Amino acid depletion (host resource)
% All current states to calculate the relative rate of change ?????(earliest switch occurs at checksize (100))?????
rr = r_con(2525:2531); % contiunous reaction rates
r_index = find(rr); % which reactions are continuous
r_new = rr(r_index); % r_new reacion rates are all non-zero
% last molecule level of continuous states:
w = x(i,2525:2531)';
w_new=w(r_index);
% time to change 100% in one state:
t_norm=w_new./r_new;
% maximum step: minimum of max_step and tolerance percentage multiplied with the smallest time
maxstep=min([max_step;tol*min(abs(t_norm))]);
% maxstep=min(abs(x_norm));
end; %End of calculating the maximum step size.
check_step = 0;
% All stochastic reaction rates:
r=a.*flag;
% Sum of it:
a0 =sum(r);
if (a0 == 0) %If there are no stochastic reactions
step = maxstep;
check_step = 1;
else
%Find timestep:
step = (1/a0)*log(1/r1);
%Find reaction:
m=sum(cumsum(r)<=r2*a0)+1; %Calculating which reaction occurs
if step > maxstep
check_step = 1;
%If the step is too big, only let it go a little bit.
%This is fine because the timesteps are exponentially distributed (and therefore memoryless).
step = maxstep;
m=0;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%END OF STOCHASTIC CALCULATION%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save current time for check of delayed reactions:
current_time=T(i)+step;
event_times=[et1(1) et2(1) et3(1) et4(1) et5(1) et6(1) et7(1) et8(1) et9(1) et10(1) et11(1) et12(1) et13(1) et14(1) et15(1)]; % All first event times
et=min(event_times); % first event time of all first event times
n=find(event_times==et); % which event comes first
n=n(1);
if et>current_time % if stored events are later than the next reaction:
%New states:
if m==0 % If no stochastic reaction (time step too big for continuous update) do nothing
else % else
z=v1(m,:); % Take the m-th row of the stoichiometric matrix v
end
% If reactions are delayed, store the completion of the reaction in seperate event vectors and move the pointer of the next free
% storage postition (makes storing and accesing the next event faster)
% m==2519...2525: transcription reactions
% m==1 and m==1260: replication reactions
if m==2519
et1(pos(1))=current_time+tau(1);
et9(pos(9))=current_time+tau(1);
pos(1)=pos(1)+1;
pos(9)=pos(9)+1;
end;
if m==2520
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et10(pos(10))=current_time+tau(2);
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(10)=pos(10)+1;
end;
if m==2521
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et11(pos(11))=current_time+tau(3);
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(11)=pos(11)+1;
end;
if m==2522
et4(pos(4))=current_time+tau(4);
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et12(pos(12))=current_time+tau(4);
pos(4)=pos(4)+1;
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(12)=pos(12)+1;
end;
if m==2523
et5(pos(5))=current_time+tau(5);
et4(pos(4))=current_time+tau(4);
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et13(pos(13))=current_time+tau(5);
pos(5)=pos(5)+1;
pos(4)=pos(4)+1;
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(13)=pos(13)+1;
end;
if m==2524
et6(pos(6))=current_time+tau(6);
et5(pos(5))=current_time+tau(5);
et4(pos(4))=current_time+tau(4);
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et14(pos(14))=current_time+tau(6);
pos(6)=pos(6)+1;
pos(5)=pos(5)+1;
pos(4)=pos(4)+1;
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(14)=pos(14)+1;
end;
if m==1
et7(pos(7))=current_time+tau(7);
et15(pos(15))=current_time+tau(7);
pos(7)=pos(7)+1;
pos(15)=pos(15)+1;
end;
if m==1260
et8(pos(8))=current_time+tau(7);
et15(pos(15))=current_time+tau(7);
pos(8)=pos(8)+1;
pos(15)=pos(15)+1;
end;
else %if there are stored events that are earlier than the next reaction
step = et - T(i); % update to the next event time
z = v2(n,:); % add delayed reaction
% delete the first event of the event vector and move the pointer of the last event one spot to the left
switch n
case{1}
et1 = [et1(2:length(et1));inf];
pos(1)=pos(1)-1;
case{2}
et2 = [et2(2:length(et2));inf];
pos(2)=pos(2)-1;
case{3}
et3 = [et3(2:length(et3));inf];
pos(3)=pos(3)-1;
case{4}
et4 = [et4(2:length(et4));inf];
pos(4)=pos(4)-1;
case{5}
et5 = [et5(2:length(et5));inf];
pos(5)=pos(5)-1;
case{6}
et6 = [et6(2:length(et6));inf];
pos(6)=pos(6)-1;
case{7}
et7 = [et7(2:length(et7));inf];
pos(7)=pos(7)-1;
case{8}
et8 = [et8(2:length(et8));inf];
pos(8)=pos(8)-1;
case{9}
et9 = [et9(2:length(et9));inf];
pos(9)=pos(9)-1;
case{10}
et10 = [et10(2:length(et10));inf];
pos(10)=pos(10)-1;
case{11}
et11 = [et11(2:length(et11));inf];
pos(11)=pos(11)-1;
case{12}
et12 = [et12(2:length(et12));inf];
pos(12)=pos(12)-1;
case{13}
et13 = [et13(2:length(et13));inf];
pos(13)=pos(13)-1;
case{14}
et14 = [et14(2:length(et14));inf];
pos(14)=pos(14)-1;
case{15}
et15 = [et15(2:length(et15));inf];
pos(15)=pos(15)-1;
end;
end;
% Continuous Langevin reaction rate calculations
if flag_mult==0 % if there are continuous reactions
rn = randn([6,1]); % normally distributed random number
for j=2525:2530
%% avoid some nuisance generation of complex values;
%% jbr, 12/21/2008
D = real(rd(j)*step);
D = max(0, D);
sd(j)=rd(j)*step+sqrt(D)*rn(j-2524); % first order Langevin equation
end
% set Langevin reaction rates
s_con(2525)=sd(2525);
s_con(2526)=sd(2526);
s_con(2527)=sd(2527);
s_con(2528)=sd(2528);
s_con(2529)=sd(2529);
s_con(2530)=sd(2530);
s_con(2531)=-sd(2525)*444-sd(2526)*274-sd(2527)*279-sd(2528)*557-sd(2529)*240-sd(2530)*2127;
r_con=s_con;
else
r_con=zeros(2531,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%delta x (stochastic and Langevin part together)
z=z+r_con';
%update all states
x(i+1,:) = x(i,:) + z;
% update all partially encapsidated genomes (for later calculations)
g(i+1,1)=sum(x(i,2:1258));
g(i+1,2)=sum(x(i,1261:2517));
%update time
T(i+1) = T(i) + step;
% if (mod(i,2000)==0)
% fprintf(1,'i=%3d T=%3.2f ',i,T(i));
% end;
% jump to next implementation part when tend is reached or all amino
% acids are depleted (which actually never happens)
if T(i+1)>tend || x(i+1,2531)<=0
break
end
end
%% Second loop
% Initialization:
start=i; % first iteration of the second part of the simulation
maxi2=120000; % maximum number of iterations in the second part of the simulation
tend=tfinal; % maximum simulated time, which is always smaller than t(start+120000)
if (sum(x(i,1260:2517))+sum(x(i,1:1258)))==0 % Only when there is no naked strand --> there might not me a quasi-steady state yet
warning_message=str('Possibly, steady-state not yet reached.')
alpha(i)=1;
tau_encap=inf;
else % Calculation of the mean alpha of the quasi-steady state distribution of N
alpha(i) = k2 * x(i,2519) / (k3*(sum(x(i,1260:2517))+sum(x(i,1:1258))));
tau_encap=1258/(k3*alpha(i)); % Average of the total encapsidation time (approximation): 1258 * r_encap
end
index1=find(x(i,2:1258)); % Find all partially encapsidated (-)RNA strands
for j=1:length(index1) % Set the remaining time for the full encapsidation of the partially encapsidated (-)RNA strands
et16(j)=T(i)+(1258-index1(length(index2)+1-j))/1258*tau_encap;
pos(16)=pos(16)+1; % Increase the pointer of the number of delayed (-)RNA encapsidations by 1
end
index2=find(x(i,1261:2517)); % Find all partially encapsidated (+)RNA strands
for j=1:length(index2) % Set the remaining time for the full encapsidation of the partially encapsidated (+)RNA strands
et17(j)=T(i)+(1258-index2(length(index2)+1-j))/1258*tau_encap;
pos(17)=pos(17)+1; % Increase the pointer of the number of delayed (+)RNA encapsidations by 1
end
tau_encap_old=tau_encap; % Set the old encapsidation time tau_encap_old to tau_encap (which will be used to update the end time of encapsidation)
% We reduce the number of states from more than 2500 to 19, because we do
% not model all intermediate states of the encapsidation (faster implementation)
x(:,1)=x(:,1); % (-)RNA (naked)
x(:,2)=g(:,1); % all partially encapsidated (-)RNA genomes
x(:,3)=x(:,1259); % fully encapsidated (-)RNA genome
x(:,4)=x(:,1260); % (+)RNA (naked)
x(:,5)=g(:,2); % all partially encapsidated (-)RNA genomes
x(:,6)=x(:,2518); % fully encapsidated (-)RNA genome
x(:,7:19)=x(:,2519:2531); % all other states: mRNAs and proteins and also the host resource amino acids in x(:,19)
x=[x(:,1:19);zeros(maxi2-i,19)]; % Reserve memory for all states and iterations
x_tau = zeros(7,19); % All past states at t-tau (tau = delay for all seven different delay times)
T=[T;zeros(maxi2-i,1)]; % Time vector
alpha=[alpha;zeros(maxi2-i,1)]; % Mean of N distribution vector
genomes=[genomes;zeros(maxi2-i,1)]; % genomes vector
a=zeros(19,1); % Calculation of all reaction rates:
ap=zeros(19,1); % Calculation of past reaction rates:
rd=zeros(19,1); % all continuous reaction rates
r_con=zeros(19,1); % all continuous reaction rates that do not take part in the stochastic time step / reaction calculation
sd = zeros(19,1); % all Langevin reaction rates
s_con=zeros(19,1); % all Langevin reaction rates that do not take part in the stochastic time step / reaction calculation
% v1: stoichiometric matrix of the non-delayed reactions
% v2: stoichiometric matrix of the delayed reactions
v1=zeros(19,19);
v2=zeros(17,19);
v1(1,18)=-1;
v1(2,1)=-1;v1(2,2)=1;
v1(4,18)=-1;
v1(5,4)=-1;v1(5,5)=1;
v1(7:12,18)=-1;
for j=13:18
v1(j,j)=1;
end
v1(13,19)=-444;
v1(14,19)=-274;
v1(15,19)=-279;
v1(16,19)=-557;
v1(17,19)=-240;
v1(18,19)=-2127;
for j=7:12
v2(j-6,j)=1;
end
v2(7,1)=1;v2(7,18)=1;
v2(8,4)=1;v2(8,18)=1;
v2(9:15,18)=1;
v2(16,2)=-1;v2(16,3)=1;
v2(17,5)=-1;v2(17,6)=1;
%checks for fast and slow reactions, Default to slow.
flag = ones(19,1); % flags for non-delayed reactions, 1 for stochastic, 0 for Langevin implementation
flagp= ones(19,1); % flags for delayed reactions, 1 for stochastic, 0 for Langevin implementation
%% Second Main Loop
for i=start:maxi2
genomes(i)=sum(x(i,1:2)); % all (-)RNA genomes that are not fully encapsidated (still produce mRNAs and not used for replication)
genomes2=(x(i,1)+x(i,2)+x(i,4)+x(i,5)); % all genomes that are involved in reactions with the N protein (only used for Poisson distribution of N)
% If there are no genomes that react with N, the genome level is set to N, such that we do not get an error message for dividing by 0
% It will not have an effect, as genomes2 is only used to calculate the mean of the distribution and the time delay, but not for the reaction
% rate, thus it will not have any effect on the stochastic reaction and time step
if genomes2==0
genomes2=1;
end
% Mean of the Poisson distribution of N
alpha(i) = k2 * x(i,7) / (k3*genomes2);
% Average time of encapsidation (1258 * delta_t) which is our way to simplify all encapsidation reactions to one reaction
tau_encap=1258/(k3*alpha(i));
% Update all encapsidation completion times with the newly calculated tau_encap
et16(1:pos(16))=T(i)+(et16(1:pos(16))-T(i))./tau_encap_old.*tau_encap; % encapsidation completion times for the (-)RNA strands
et17(1:pos(17))=T(i)+(et17(1:pos(17))-T(i))./tau_encap_old.*tau_encap; % encapsidation completion times for the (+)RNA strands
tau_encap_old=tau_encap; % set the old encapsidation time to the new one
% When the mean of the poisson distribution is big, N will be drawn from a normal distribution
if alpha(i)>=20
%% avoid some nuisance generation of complex values;
%% jbr, 12/21/2008
D = real(alpha(i));
D = max(0, D);
x(i,13) = round(alpha(i) + sqrt(D) * randn);
else
% Draw N from Poisson distribution
r3 = rand;
for j=1:length(facts)
%% avoid some nuisance generation of complex values;
%% jbr, 12/21/2008
D = real(alpha(i));
D = max(0, D);
p(j)=exp(-alpha(i))*1/facts(j)*D^(j-1); % Probability of N being 0,1,2,....
r3=r3-p(j);
if r3<0
x(i,13)=j-1;
break;
end;
end
end
% Stop simulation when tend (4 hours) is reached
if T(i)>tend
break
end
% Calculation of all reaction rates:
a(1) = 50*k1(1)*x(i,6)*x(i,18);
a(2) = k3*x(i,13)*x(i,1);
a(3) = 0;
a(4) = k1(1)*x(i,3)*x(i,18);
a(5) = k3*x(i,13)*x(i,4);
a(6) = 0;
a(7:12) = k1(1:6).*genomes(i).*x(i,18);
% If level of N is greater than the checksize (100), N will be created by a continuous Langevin reaction and not be drawn from the Poisson distribution
if x(i,13) there are continuous reactions
if flag_mult==0 % if there are continuous reactions
% Calculate the maximum step size:
% all continuous reaction rates
r_con=zeros(19,1); % r_con is used to calculate the next maximum continuous step (set to 1% change in the continuous states)
% this has to be done before the Langevin calculations, because we need the step size for the Langevin reaction rate calculations
% rd will later be used for the Langevin reaction rate calculation (only when the reactions are continous)
for j=13:19
rd(j)=a(j)*(1-flag(j));
end
r_con(13)=rd(13);
r_con(14)=rd(14);
r_con(15)=rd(15);
r_con(16)=rd(16);
r_con(17)=rd(17);
r_con(18)=rd(18);
r_con(19)=-rd(14)*274-rd(15)*279-rd(16)*557-rd(17)*240-rd(18)*2127; % Amino acid depletion (host resource)
% all current states to calculate the relative rate of change (I take the max of (100,w) because we don't want to divide through 0)
rr = r_con(13:19); % contiunous reaction rates
r_index = find(rr); % which reactions are continuous
r_new = rr(r_index); % r_new reacion rates are all non-zero
% last molecule level of continuous states:
w = x(i,13:19)';
w_new=w(r_index);
% time to change 100% in one state:
t_norm=w_new./r_new;
% maximum step: minimum of max_step and tolerance percentage multiplied with the smallest time
maxstep=min([max_step;tol*min(abs(t_norm))]);
end; %End of calculating the maximum step size.
check_step = 0;
% All stochastic reaction rates:
r=a.*flag;
% Sum of it:
a0 =sum(r);
if (a0 == 0) %If there are no stochastic reactions
step = maxstep;
check_step = 1;
else
%Find timestep:
step = (1/a0)*log(1/r1);
%Find reaction:
m=sum(cumsum(r)<=r2*a0)+1; %Calculating which reaction occurs
if step > maxstep
check_step = 1;
%If the step is too big, only let it go a little bit.
%This is fine because the timesteps are exponentially distributed (and so memoryless).
step = maxstep;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%END OF STOCHASTIC CALCULATION%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save current time for check of delayed reactions:
current_time=T(i)+step;
event_times=[et1(1) et2(1) et3(1) et4(1) et5(1) et6(1) et7(1) et8(1) et9(1) et10(1) et11(1) et12(1) et13(1) et14(1) et15(1) et16(1) et17(1)]; % All first event times
et=min(event_times); % first event time of all first event times
n=find(event_times==et); % which event comes first
n=n(1);
if et>current_time % if stored events are later than the next reaction:
%New states:
if m==0 % If no stochastic reaction (time step too big for continuous update) do nothing
else
z=v1(m,:); % else take the m-th row of the stoichiometric matrix v
end
% If reactions are delayed, store the completion of the reaction in seperate event vectors and move the pointer of the next free
% storage postition (makes storing and accesing the next event faster)
% m==7...12: transcription reactions
% m==1 and m==4: replication reactions
% m==2 and m==5: encapsidation reactions
if m==7
et1(pos(1))=current_time+tau(1);
et9(pos(9))=current_time+tau(1);
pos(1)=pos(1)+1;
pos(9)=pos(9)+1;
end;
if m==8
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et10(pos(10))=current_time+tau(2);
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(10)=pos(10)+1;
end;
if m==9
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et11(pos(11))=current_time+tau(3);
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(11)=pos(11)+1;
end;
if m==10
et4(pos(4))=current_time+tau(4);
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et12(pos(12))=current_time+tau(4);
pos(4)=pos(4)+1;
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(12)=pos(12)+1;
end;
if m==11
et5(pos(5))=current_time+tau(5);
et4(pos(4))=current_time+tau(4);
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et13(pos(13))=current_time+tau(5);
pos(5)=pos(5)+1;
pos(4)=pos(4)+1;
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(13)=pos(13)+1;
end;
if m==12
et6(pos(6))=current_time+tau(6);
et5(pos(5))=current_time+tau(5);
et4(pos(4))=current_time+tau(4);
et3(pos(3))=current_time+tau(3);
et2(pos(2))=current_time+tau(2);
et1(pos(1))=current_time+tau(1);
et14(pos(14))=current_time+tau(6);
pos(6)=pos(6)+1;
pos(5)=pos(5)+1;
pos(4)=pos(4)+1;
pos(3)=pos(3)+1;
pos(2)=pos(2)+1;
pos(1)=pos(1)+1;
pos(14)=pos(14)+1;
end;
if m==1
et7(pos(7))=current_time+tau(7);
et15(pos(15))=current_time+tau(7);
pos(7)=pos(7)+1;
pos(15)=pos(15)+1;
end;
if m==4
et8(pos(8))=current_time+tau(7);
et15(pos(15))=current_time+tau(7);
pos(8)=pos(8)+1;
pos(15)=pos(15)+1;
end;
if m==2
et16(pos(16))=current_time+tau_encap;
pos(16)=pos(16)+1;
end;
if m==5
et17(pos(17))=current_time+tau_encap;
pos(17)=pos(17)+1;
end;
else %if there are stored events that are earlier than the next reaction
step = et - T(i); % update to the next event time
z = v2(n,:); % add delayed reaction
% delete the first event of the event vector and move the pointer of the last event one spot to the left
switch n
case{1}
et1 = [et1(2:length(et1));inf];
pos(1)=pos(1)-1;
case{2}
et2 = [et2(2:length(et2));inf];
pos(2)=pos(2)-1;
case{3}
et3 = [et3(2:length(et3));inf];
pos(3)=pos(3)-1;
case{4}
et4 = [et4(2:length(et4));inf];
pos(4)=pos(4)-1;
case{5}
et5 = [et5(2:length(et5));inf];
pos(5)=pos(5)-1;
case{6}
et6 = [et6(2:length(et6));inf];
pos(6)=pos(6)-1;
case{7}
et7 = [et7(2:length(et7));inf];
pos(7)=pos(7)-1;
case{8}
et8 = [et8(2:length(et8));inf];
pos(8)=pos(8)-1;
case{9}
et9 = [et9(2:length(et9));inf];
pos(9)=pos(9)-1;
case{10}
et10 = [et10(2:length(et10));inf];
pos(10)=pos(10)-1;
case{11}
et11 = [et11(2:length(et11));inf];
pos(11)=pos(11)-1;
case{12}
et12 = [et12(2:length(et12));inf];
pos(12)=pos(12)-1;
case{13}
et13 = [et13(2:length(et13));inf];
pos(13)=pos(13)-1;
case{14}
et14 = [et14(2:length(et14));inf];
pos(14)=pos(14)-1;
case{15}
et15 = [et15(2:length(et15));inf];
pos(15)=pos(15)-1;
case{16}
et16 = [et16(2:length(et16));inf];
pos(16)=pos(16)-1;
case{17}
et17 = [et17(2:length(et17));inf];
pos(17)=pos(17)-1;
end;
end;
% Continuous Langevin reaction rate calculations
if flag_mult==0 % If there are continous reactions
rn = randn([5,1]); % normally distributed random number
for j=14:18
%% avoid some nuisance generation of complex values;
%% jbr, 12/21/2008
D = real(rd(j)*step);
D = max(0, D);
sd(j)=rd(j)*step+sqrt(D)*rn(j-13); % First order Langevin reaction rate
end
s_con(14)=sd(14);
s_con(15)=sd(15);
s_con(16)=sd(16);
s_con(17)=sd(17);
s_con(18)=sd(18);
s_con(19)=-sd(14)*274-sd(15)*279-sd(16)*557-sd(17)*240-sd(18)*2127;
r_con=s_con;
else
r_con=zeros(19,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%delta x (stochastic and Langevin part together)
z=z+r_con';
%update species
x(i+1,:) = x(i,:) + z;
%update time
T(i+1) = T(i) + step;
% if (mod(i,5000)==0)
% fprintf(1,'i=%3d T=%3.2f ',i,T(i));
% end;
end %main loop
message = ['End of run ',num2str(u),' with ',num2str(i),' iterations']
index_old=1;
index=zeros(tfinal+1,1);
for j=1:length(time_vec)
index(j)=searcht(index_old,T,time_vec(j),i,0)-1;
index_old=index(j);
end
% Save all states of all runs in new vectors
for c=1:18
eval(['x',num2str(c),'(:,u)=x(index,c);'])
end
end
save xmean.debug
g1 = x1+x2+x3; % all (-)RNA genomes
g2 = x4+x5+x6; % all (+)RNA genomes
x_mean1 = [mean(g1,2), mean(g2,2)]; % mean of genomes
x_mean2 = [mean(x7,2), mean(x8,2), mean(x9,2), mean(x10,2), mean(x11,2), mean(x12,2)]; % mean of mRNAs
x_mean3 = [mean(x13,2), mean(x14,2), mean(x15,2) mean(x16,2), mean(x17,2), mean(x18,2)]; % mean of proteins
f1=0:1400:140000; % bins for the GFP distribution
GFP35=x17(12601,:); % GFP level at 3.5 hours
GFP4=x17(14401,:); % GFP level at 4 hours
[zz,yy]=hist(GFP35,f1); % create history of GFP level at 3.5 hours
table4=[yy' zz']; % distibution table of GFP at 3.5 hours
[zz,yy]=hist(GFP4,f1); % create history of GFP level at 4 hours
table5=[yy' zz']; % distibution table of GFP at 4 hours
f2=0:30:3000; % bins for the RNA distribution
RNA25=g1(9001,:); % (-)RNA level at 2.5 hours
RNA3=g1(10801,:); % (-)RNA level at 3 hours
RNA4=g1(14401,:); % (-)RNA level at 4 hours
[zz,yy]=hist(RNA25,f2); % create history of RNA level at 2.5 hours
table6=[yy' zz']; % distibution table of RNA at 2.5 hours
[zz,yy]=hist(RNA3,f2); % create history of RNA level at 3 hours
table7=[yy' zz']; % distibution table of RNA at 3 hours
[zz,yy]=hist(RNA4,f2); % create history of RNA level at 4 hours
table8=[yy' zz']; % distibution table of RNA at 4 hours
table1 = [time_vec/3600, x_mean1];
table2 = [time_vec/3600, x_mean2];
table3 = [time_vec/3600, x_mean3];
table9 = [x12(5400,:)' g1(14401,:)']; % create table for LmRNA, (-)RNA dependency
save genomes.dat table1
save mRNAs.dat table2
save proteins.dat table3
save GFPdist35.dat table4
save GFPdist4.dat table5
save RNAdist25.dat table6
save RNAdist3.dat table7
save RNAdist4.dat table8
save LmRNA.dat table9
Main Script