Array indices must be positive integers or logical values.

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=5; % in molecules
time_log=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(t_nr)));
    prob_prod=prod_rate/(prod_rate+deg_rate*x_nr(t_nr)); % Propensity of production
    
    t=t+wait_time; % Update current time
    t_nr=t_nr+1; % Update the number of steps (reactions) associated with the experiment.
    time_log(t_nr)=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(t_nr)=x_nr(t_nr-1)+1; % Implements production
    else
        x_nr(t_nr)=x_nr(t_nr-1)-1; % Implements degradation
    end
end

if t>t_max
    t_nr(end)=[];
    x_nr(end)=[];
end
close all
plot(time_log,t_nr,'r','LineWidth',2)
hold on
plot(time_log,x_nr,'r','LineWidth',2)

I get the error Array indices must be positive integers or logical values in line 24:

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

Does anyone know why this is?

[DGM]

You're trying to index into x_nr

t_nr=0; % in s
x_nr=5; % in molecules
%  ...
x_nr(t_nr)

t_nr is zero. Zero is not a valid array index, as the message states.

x_nr is a scalar. it only has one element, so the only valid value for t_nr is 1

[Jonathan Louie]

I understand now, thank you so much.

I fixed the script, but x_nr is supposed to osscilate around 200 within 500 seconds but my graph osscilates around 30. Do you know why?

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

end

if t>t_max

t_nr(end)=[];

x_nr(end)=[];

end

close all

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

hold on

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

legend('t_nr','x_nr')

[Voss]

"x_nr is supposed to oscillate around 200 ..."

Try changing

deg_rate=k*x_nr; % [1/s]*[molecules] = [molecules/s]

to

deg_rate=k; % [1/s]

The units make more sense that way in this expression:

    prod_rate + (deg_rate * x_nr(step_counter))
% molecules/s + (     1/s * molecules         ) = molecules/s <- coherent units

vs:

    prod_rate + (   deg_rate * x_nr(step_counter))
% molecules/s + (molecules/s * molecules         ) <- incoherent units

And that change also makes x_nr oscillate around 200.

[ Dislcaimer: I have no idea what problem this code is supposed to solve; I'm only observing that setting deg_rate=k seems to have the desired effect (oscillation ~200) and seems to make more sense units-wise than deg_rate=k*x_nr (assuming the units given in the comments at the top are accurate). ]

[Jonathan Louie]

Thanks so much for the help. I'm trying to model the dynamics of the non-regulated construct and the auto-regulated construct. So the equation deg_rate=k*x_nr; needs to have the x_nr because x_nr stands for the number of repressor molecules in the system. Also both t_nr and x_nr are vectors containing time series. So I'm using the Gillespie algorithm and plotting the time vector, t_nr, against the non-regulated fluorescence, x_nr

[Voss]

In that case, check that the units of your variables are correct.

[Jonathan Louie]

I checked the units and I understand that the deg rate should be 1/s molecules/s. But the instructions given for my assignment is:

the degradation rate of both models will be given by: deg(x(t)) = kx(t) where x(t) stands for the number of repressor molecules in the system. k is given as 0.01s^-1, but what I'm unsure of is if the statement below is referring to the variable for number of repressor molecules.

Define a variable (x_nr) with the number of molecules of the non-regulated system and initialize it with value 5 molecules. As for tnr &, xnr will be a vector that contains the changes in number of molecules corresponding to the different time points of vector tnr.

Because that's why I plugged in x_nr for x(t) for the deg rate equation.

Also for reference, I'm trying to plot the dynamics of the non-regulated construct and the auto-regulated construct by running the Gillespie algorithm until we reach the maximum time, t_max. To do this, set up a while loop so that the state of the system’s dynamics – its transitions – is simulated from 0 s up to time t_max, then drawing a wait time dt from an Negative exponential distribution (Poisson process).

Thank you so much for any help I really appreciate it.

[Image Analyst]

And what about the help you got below, you know, in the official Answers section of this page? Did you see it? If not, scroll down.

Answer 1

[Image Analyst]

Since this is one of the most common FAQs, see the FAQ for a thorough discussion:

https://matlab.fandom.com/wiki/FAQ#.22Subscript_indices_must_either_be_real_positive_integers_or_logicals..22

Answer 2

[Torsten]

if rand < prob_prod % Defines whether production takes place based on Monte Carlo method.
    x_nr(t_nr)=x_nr(t_nr-1)+1; % Implements production
else
    x_nr(t_nr)=x_nr(t_nr-1)-1; % Implements degradation
end

The first time you enter this if-clause, t_nr = 1.

Thus you address x_nr(t_nr-1) = x_nr(0).

This throws an error in MATLAB since array indices start with index 1, not 0.