Trying to model the dynamics of the non-regulated construct and the auto-regulated construct using the Gillespie algorithm, do you know why it doesn't reach 200 and stay?

t_max= 1200; % in s

k=0.01; % in s^-1

g_nr=2; % in molecules/s

t_nr=0; % in s

x_nr=5; % in molecules

step_counter=1;

% Simulation of the non-regulated construct

t=0;

prod_rate=g_nr;

deg_rate=k*x_nr;

while t<t_max

wait_time=-log(rand)/((prod_rate)+(deg_rate*x_nr(step_counter)));

prob_prod=(prod_rate)/((prod_rate)+(deg_rate*x_nr(step_counter))); % Propensity of production

t=t+wait_time; % Update current time

step_counter=step_counter+1; % Update the number of steps (reactions) associated with the experiment.

t_nr(step_counter)=t; % Add the current time to the time log

if rand < prob_prod % Defines whether production takes place based on Monte Carlo method.

x_nr(step_counter)=x_nr(step_counter-1)+1; % Implements production

else

x_nr(step_counter)=x_nr(step_counter-1)-1; % Implements degradation

end

if t>t_max

t_nr(end)=[];

x_nr(end)=[];

end

end

close all

figure,

plot(t_nr,'r','LineWidth',2)

hold on

plot(x_nr,'b','LineWidth',2)

legend('t_nr','x_nr')

x_nr should oscillate around a value of 200 within ~500 s after the start of the experiment, but my code only goes up to approximately 90. The production and degradation rate as well as the probability of production are all given as equations that I'm sure I got right, so I'm not sure what's causing it to only go up to 90 molecules.

Answer 1

[Jeffrey Clark]

@Jonathan Louie, I believe your deg_rate=k*x_nr should just be deg_rate=k.