Figure 2a:

Stochastic simulation run including all encapsidation (chain)

Code for Figure 2a

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=1;                 % number of simulation runs

% 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)

% 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=100000;         % maximum number of iterations in the first part of the simulation
tend=14400;         % maximum simulated time

% 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
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)

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
        sd(j)=rd(j)*step+sqrt(rd(j)*step)*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 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


message = ['End of run ',num2str(u),' with ',num2str(i),' iterations']

end

%% store results at times given in tplot
%% use 1 second time interval for plotting

%tplot = linspace(T(1), T(end), 1)';
%% jbr, removing this linspace error introduced in changest 32
%% 12/25/2008
tplot = (T(1):1:T(end))';
[tcon, idx] = sort([T; tplot]);
test = (idx > length(T));
howmanyoldguys = cumsum(1-test);
trow  = howmanyoldguys(test);

x1 = x(trow,1);                    % (-)RNA_0 genome
x2 = sum(x(trow,1:1259),2);        % all (-)RNA genomes
x3 = x(trow,1259);                 % (-)RNA_{1258}

x4 = x(trow,1260);                 % (+)RNA genome
x5 = sum(x(trow,1260:2518),2);     % all (+)RNA genomes 
x6 = x(trow,2518);                 % (+)RNA_{1258}

x7 = x(trow,2519);                 % N mRNA
x8 = x(trow,2525);                 % N proteins


table = [tplot/3600, x1, x2, x3, x4, x5, x6 x7 x8] ;
save fullrun.dat table