Plotting the dynamics of the non-regulated construct and the auto-regulated construct utilizing the Gillespie algorithm and perturbations

I'm trying to model the dynamics of the non-regulated construct and the auto-regulated construct utilizing the Gillespie algorithm. There's 2 perturbations implemented at 400 and 800 seconds. I was wondering if that's the correct way to code adding 150 molecules at 400 seconds only one time, and subtracting 150 molecules at 800 seconds one time?

Also my graph extends until 5000 on the x-axis, and I want it to end at 1200, does anyone know what I need to change in my code? Thanks!

clearvars

% Setup

t_max= 1200; % in s

k=0.01; % in s^-1

g_nr=2; % in molecules/s

t_nr=0; % in s

x_nr(1)=5; % in molecules

step_counter=1;

% Perturbations

perturbation_magnitude=150; % in molecules

perturbation_time_400=400;% in s

perturbation_time_800=800;% in s

is_perturbation_400=0;

is_perturbation_800=0;

% 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

if t>perturbation_time_400

if is_perturbation_400==0;

x_nr=[x_nr perturbation_magnitude];

is_perturbation_400=1;

if t>perturbation_time_800

if is_perturbation_800==0;

x_nr=x_nr-perturbation_magnitude;

is_perturbation_800=1;

end

end

end

end

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

% autoregulated construct

t_max= 1200; % in s

k=0.01; % in s^-1

step_counter=1;

t=0;

g_nr=2;

x_nr(1)=5;

g_ar=4;

n_hill=8;

x_ar(1)=5;

t_ar=0;

theta=(g_nr)/k;

prod_ratea=g_ar*((theta^n_hill)/((theta^n_hill)+(x_ar^n_hill)));

deg_rate=k*x_ar;

while t<t_max

wait_time=-log(rand)/((prod_ratea)+(deg_rate*x_ar(step_counter)));

prob_prod=prod_ratea/((prod_ratea)+(deg_rate*x_ar(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_ar(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_ar(step_counter)=x_ar(step_counter-1)+1; % Implements production

else

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

if t>perturbation_time_400

if is_perturbation_400==0;

x_ar=[x_ar perturbation_magnitude];

is_perturbation_400=1;

if t>perturbation_time_800

if is_perturbation_800==0;

x_ar=x_ar-perturbation_magnitude;

is_perturbation_800=1;

end

end

end

end

end

if t>t_max

t_ar(end)=[];

x_ar(end)=[];

end

end

figure,

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

hold on

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

legend('t_ar','x_ar')

figure,

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

hold on

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

legend('x_ar','x_nr')

Answer 1

[sai charan sampara]

Hello Jonathan,

I understand that you want to add 150 molecules exactly at 400 seconds and remove 150 molecules exactly at 800 seconds.

This can be done as follows. You are storing all the values of the time steps “t” in an array “t_nr”. This array can be used to access previous time step value. Use an if conditional to check if the current time step is greater than 400 and the previous time step from the array is less than 400. If this is satisfied, then it indicates that the time 400 seconds was just passed. This means that 150 molecules should be added at this instant. So, at this instance check whether the “is_perturbation” flag is still zero. Once this is satisfied add 150 to the current value of “x_nr” that is “x_nr(step_count)”. To prevent adding multiple times once the molecules are added the flag is set to 1.

if  t>perturbation_time_400 && t_nr(step_counter-1)<400 && is_perturbation_400==0
    x_nr(step_counter)=x_nr(step_counter-1)+150;
    is_perturbation_400=1;
end

Similar method can be applied at 800 seconds as well.

if  t>perturbation_time_800 && t_nr(step_counter-1)<800 && is_perturbation_800==0
    x_nr(step_counter)=x_nr(step_counter-1)-150;
    is_perturbation_800=1;
end

One issue with this is that once 150 molecules are added the “x_nr” values become big so thereby reducing the “wait_time” values to very small numbers this can increase the number of times the loop is run significantly. This will increase run time.

In the plot the values on the x axis correspond to the indexes of the plotted arrays. Since no x variable was mentioned in the “plot” the index values are being used and the 5000 is because the size of “x_nr” is nearly 5000. To change this, use the plot function with “t_nr” as the x variable. Using “hold on” makes them plot on the same graph but as y variables and hence the unwanted behaviour. To fix this plot the variables against “t_nr” and “t_ar” instead of along them. The three plots can be changed as follows to obtain the correct graphs.

For first figure:

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

For the next two figures:

figure,
plot(t_ar,x_ar,'r','LineWidth',2) 

figure,
plot(t_ar,x_ar,'r','LineWidth',2) 
hold on
plot(t_nr,x_nr,'c','LineWidth',2)
legend('x_ar','x_nr')

You can refer to the below documentation to learn more about plot:

https://www.mathworks.com/help/matlab/ref/plot.html

Hope this helps.

Thanks,

Charan.