I would like to simulate birth process with time delay. Here is a reference ref on page 3. Note there is two species whereas in my case I have one species. So only equation (1) and (2) should be considered. My code does not give desired result. Any suggestion.

Here is the code for Birth Death Process without delay(propensity rate depends on current population).

$\emptyset\xrightarrow[]{\text{$\lambda$ X}} X\quad \quad$ birth

$X\xrightarrow[]{\text{$\mu$ X}} \emptyset \quad \quad$ death

SeedRandom@2;
With[{λ = 4, μ = 1, initialPop = 10}, sim = NestList[(
      a1 = λ #[[2]];
      a2 = μ #[[2]];
      a0 = Total@{a1, a2};
      reaction = 1/a0 Accumulate@{a1, a2};
      pos = First@FirstPosition[reaction, _?(RandomReal[] < # &)];
      Δt = 
       RandomVariate@ExponentialDistribution[a1 + a2];

      Which[ 
       pos == 1, {#[[1]] + Δt, #[[2]] + 1},
       pos == 2, {#[[1]] + Δt, #[[2]] - 1}

       ]

      ) &, {0, initialPop}, 20]];
sim;

ListLinePlot[sim, Epilog -> {Red, PointSize[Medium], Point[sim]}, 
 Frame -> True, PlotTheme -> "Detailed", 
 FrameLabel -> {"Time", "Population"}, ImageSize -> Large, 
 InterpolationOrder -> 0]

enter image description here

Now Assume there is delay on birth reaction.

$\emptyset\xrightarrow[]{\text{$\lambda$ X}} X\quad \quad$ birth has a delay

$X\xrightarrow[]{\text{$\mu$ X}} \emptyset \quad \quad$ death has no delay

  1. Sample reaction time $t_1\sim Exp(\lambda X+\mu X)$
  2. Choose a reaction.
  3. If it is a death set $X=X-1$ and new time is $t_1$.
  4. If it is a birth sample $t_d\sim Gamma(4,2)$
  5. Put $t_1+t_d$ into queue.
  6. Sample new reaction time $t_w\sim Exp(\lambda X+\mu X)$
  7. If $t_1+t_d<t_1+t_w$ set $X=X+1$. New time $t_1+t_d$
  8. If $t_1+t_d>=t_1+t_w$
  9. Choose a reaction.
  10. If it is a death set $X=X-1$.
  11. If it is a birth set $X=X$. No birth. New time $t_1+t_w$
  12. Repeat

    SeedRandom@2;
    With[{λ = 30, μ = 1, initialPop = 10}, sim = NestList[(
      a1 = λ #[[2]];
      a2 = μ #[[2]];
      a0 = Total@{a1, a2};
      reaction = 1/a0 Accumulate@{a1, a2};
      pos = First@FirstPosition[reaction, _?(RandomReal[] < # &)];
      pos2 = First@FirstPosition[reaction, _?(RandomReal[] < # &)];
      t = RandomVariate@ExponentialDistribution[a1 + a2];
      td = RandomVariate@GammaDistribution[4, 2];
      tw = RandomVariate@ExponentialDistribution[a1 + a2];
    
      Which[ 
       pos == 2, {#[[1]] + t, #[[2]] - 1},
       pos == 1,
    
       Which[
    
        t + td < t + 
          tw, {#[[1]] + t + td, #[[2]] + 1},
    
        t + td >= t + tw,
    
        Which[
    
         pos2 == 2, {#[[1]] +t + tw, #[[2]] - 1},
    
         pos2 == 1, {#[[1]] + t + tw, #[[2]]}
    
         ]
    
        ]
    
       ]
    
      ) &, {0, initialPop}, 100]];
    sim;
    
    ListLinePlot[sim, Epilog -> {Red, PointSize[Medium], Point[sim]}, 
     Frame -> True, PlotTheme -> "Detailed", 
     FrameLabel -> {"Time", "Population"}, ImageSize -> Large, 
     InterpolationOrder -> 0]
    

enter image description here

EDIT: Let's look at Pure Birth process(there is no death).

Here is the code without delay.

SeedRandom@2
With[{A = 5, initialPop = 10}, sim = NestList[(

      Δt = 
       RandomVariate@ExponentialDistribution[A #[[2]]];

      {#[[1]] + Δt, #[[2]] + 1}

      ) &, {0, initialPop}, 10]];
ListLinePlot[sim, Epilog -> {Red, PointSize[Medium], Point[sim]}, 
 Frame -> True, PlotTheme -> "Detailed", 
 FrameLabel -> {"Time", "Population"}, ImageSize -> Large, 
 InterpolationOrder -> 0]

enter image description here

Here is the algorithm for delayed Pure Birth process

  1. Sample reaction time $t_1\sim Exp(\lambda X)$
  2. Sample delay time $t_{d_1}\sim Gamma(4,2)$
  3. Put $s_1=t_1+t_{d_1}$ into stack. Stack={$s_1$}
  4. Sample new reaction time $t_2\sim Exp(\lambda X)$
  5. If $t_1+t_2<s_1$, sample new delay time $t_{d_2}\sim Gamma(4,2)$ and Put $s_2=t_1+t_2+t_{d_2}$ into stack. Stack={$s_1$,$s_2$}. Order stack min to max. Min will be birth time if there is a birth in the future.
  6. If $t_1+t_2>=s_1$, let X=X+1. Remove $s_1$ from stack. set time $s_1$
  7. Repeat

I don't know how to code this. Any suggestion.

  • I believe the [[3]]'s in your definition of a1 and a2 should be [[2]]'s. – Chris K Feb 3 at 3:00
  • Yes they should be [[2]], I don't know how that happened. I'll fixed. – Okkes Dulgerci Feb 3 at 6:02
up vote 2 down vote accepted

I figure it out how to simulate birth-death process with one species. Black curve is DDE (delayed DE).

 A = 30;(*Birth rate*)
    β = 0.05;(*Death rate*)
    {a, b} = {72, 1/12};(*delay parameter, mean=6,var=0.5*)

   sim = Block[{X, x}, Table[

        X = x = 0;(*Initial population*)

        currentTime = 0;(*Initial time*)
        stackTime = {};
        Prepend[Join @@ Last@Reap@Do[

             a1 = A   ;
             a2 = β X;
             a0 = Total@{a1, a2};


     reactionVec = 1/a0 Accumulate@{a1, a2};

                reaction = 
  First@FirstPosition[reactionVec - RandomReal[], _?Positive];

                 tDelay = RandomVariate@GammaDistribution[a, b];
                 tWait = RandomVariate@ExponentialDistribution[a0];
                 currentTime = currentTime + tWait;

                 stackTime = Sort@stackTime;

                 minStack = Min[stackTime];

             Which[

              currentTime < minStack,


              Which[


               reaction == 1, {Sow@{currentTime, X}, 
                AppendTo[stackTime, currentTime + tDelay]},

               reaction == 2, {Sow@{currentTime, X -= 1}}],


              minStack < currentTime, 

              {Sow@{minStack, X += 1}, currentTime = minStack, 
               stackTime = Rest@stackTime} ]

             , 10000], {0, x}], 10]];
    With[{τ = 6},
      sol = NDSolveValue[{x'[t] == A UnitStep[t - τ] - β x[t], 
         x[t /; t <= 0] == 0}, x[t], {t, 0, 120}]];

    fig = Show[
      ListLinePlot[sim, Frame -> True, PlotTheme -> "Detailed", 
       FrameLabel -> {"Time", "Population"}, ImageSize -> Large], 
      Plot[sol, {t, 0, 120}, PlotStyle -> {Thick, Black}]]

enter image description here

Here we can extract population in unit interval.

  sim2 = Table[First /@ SplitBy[sim[[i]], Last], {i, Length@sim}];

iFun = Table[
   Interpolation[sim2[[i]], InterpolationOrder -> 1], {i, 
    Length@sim2}];
\[Tau] = 1;
range = Table[Range[0, sim2[[i, -1, 1]], \[Tau]], {i, Length@sim2}];
data = Table[
   Transpose@{range[[i]], Floor@iFun[[i]]@range[[i]]}, {i, 
    Length@sim2}];

fig2 = Show[
  ListLinePlot[data, Frame -> True, PlotTheme -> "Detailed", 
   FrameLabel -> {"Time", "Population"}, ImageSize -> Large], 
  Plot[sol, {t, 0, 120}, PlotStyle -> {Thick, Black}, 
   PlotRange -> All]]

enter image description here

I would like to extend my simulation for two species. Here is my attempt. It doesn't work as expected. I also have Matlab code that I don't understand. Any suggestions.

λ1 = 30;(*Birth rate of X*)
μ1 = 0.05;(*Death rate of X*)
λ2 = 30;(*Birth rate of Y*)
μ2 = 0.05;(*Death rate of Y*)
{α1, β1} = {72, 1/12}; (*mean=6,var=0.5 of X*)
{α2, β2} = {72, 1/12};
θ = 10;
n = 2;



 sim = Block[{X, x, Y, y}, Table[

X = x = 0;(*Initial population of X*)
Y = y = 0;

currentTime = 0;(*Initial time*)
  stackTimeX = {};
 stackTimeY = {};

Prepend[Join @@ Last@Reap@Do[

     a1 = λ1   ;
     a2 = μ1 X;

     a3 = λ2 X^n/(θ^n + X^n)   ;
     a4 = μ2 Y;
     a0 = Total@{a1, a2, a3, a4};
     reactionVec = 1/a0 Accumulate@{a1, a2, a3, a4};



      reaction = 
  First@FirstPosition[reactionVec - RandomReal[], _?Positive];



     nextReaction = RandomVariate@ExponentialDistribution[a0];

     currentTime = currentTime + nextReaction;

     tDelayX = 
      RandomVariate@GammaDistribution[α1, β1];
     tDelayY = 
      RandomVariate@GammaDistribution[α2, β2];

     stackTimeX = Sort@stackTimeX;
     stackTimeY = Sort@stackTimeY;

     minStack = Min[Min@stackTimeX, Min@stackTimeY];

     Which[

      currentTime < minStack,

      Which[


       reaction == 1, {Sow@{currentTime, X, Y}, 
        AppendTo[stackTimeX, currentTime + tDelayX]},

       reaction == 2, Sow@{currentTime, X -= 1, Y},


       reaction == 3, {Sow@{currentTime, X, Y}, 
        AppendTo[stackTimeY, currentTime + tDelayY]},

       reaction == 4, Sow@{currentTime, X, Y -= 1}

       ],

      minStack < currentTime,



      posStack = 
       First@Ordering[{Min@stackTimeX, Min@stackTimeY}, 1];

      Which[

       posStack == 1,

       {Sow@{minStack, X += 1, Y}, currentTime = minStack, 
        stackTimeX = Rest@stackTimeX;},

       posStack == 2,

       {Sow@{minStack, X, Y += 1}, currentTime = minStack, 
        stackTimeY = Rest@stackTimeY;}



       ]], 10000], {0, x, y}], 1]];




 {simX, simY} = {sim[[All, All, {1, 2}]], sim[[All, All, {1, 3}]]};

plot1 = Show[
   ListLinePlot[simY, PlotStyle -> {ColorData[97, "ColorList"][[2]]}, 
    PlotLegends -> {"Y"}, Frame -> True, PlotTheme -> "Detailed", 
    FrameLabel -> {"Time", "Population"}, ImageSize -> Large, 
    InterpolationOrder -> 0], 
   ListLinePlot[simX, PlotLegends -> {"X"}, 
    PlotStyle -> {ColorData[97, "ColorList"][[1]]}, 
    InterpolationOrder -> 0]];

With[{τ1 = τ2 = 6},
  {solX, solY} = 
   NDSolveValue[{x'[
       t] == λ1 UnitStep[t - τ1] - μ1 x[t], 
     y'[t] == λ2 (x[t - τ2])^
        n/(θ^n + (x[t - τ2])^n) - μ2  y[t], 
     x[t /; t <= 0] == 0, y[t /; t <= 0] == 0}, {x, y}, {t, 0, 70}]];

Show[plot1, 
 Plot[{solX[t], solY[t]}, {t, 0, 70}, 
  PlotStyle -> {Thick, Black, Red}]]

enter image description here

Here is the Matlab code and output.

    function dydt = Delayed_2species(t,y,Z,par)

lambda1 = par.lambda1;
lambda2 = par.lambda2;
mu1     = par.mu1;
mu2     = par.mu2;
n       = par.n;
theta   = par.theta;
tau_1   = par.tau_1;

ylag = Z(:,1);
dydt=[lambda1*(t>=tau_1) - mu1*y(1)
      lambda2*ylag(1)^n/(theta^n+ylag(1)^n) - mu2*y(2)];

clearvars
clc

%% simulation of the delay differential equation
% using Matlab's dde23

% define parameters
lambda1 = 30;
lambda2 = 30;
mu1     = 0.05;
mu2     = 0.05;
tau_1   = 6;
tau_2   = 6;
n_Hill  = 2;
theta   = 10;

pars.lambda1 = lambda1;
pars.lambda2 = lambda2;
pars.mu1     = mu1;
pars.mu2     = mu2;
pars.tau_1   = tau_1;
pars.n       = n_Hill;
pars.theta   = theta;

lags    = tau_2;
history = [0;0];
t_final = 100;
tspan   = [0 t_final];

sol     = dde23(@Delayed_2species,lags,history,tspan,[],pars);
t_sol   = sol.x;
species = sol.y; % in row format unlike ode45

figure; hold on
h_X=plot(t_sol,species(1,:),'LineWidth',2);
h_Y=plot(t_sol,species(2,:),'LineWidth',2);
legend([h_X h_Y],'X','Y')

%% Simulate using delayed Gillespie algorithm
% define parameters
lambda1 = 30;
lambda2 = 30;
mu1     = 0.05;
mu2     = 0.05;
n_Hill  = 2;
theta   = 10;
alpha1  = 72;
beta1   = 1/12;
alpha2  = 72;
beta2   = 1/12;


%transition matrix for Gillespie algorithm (column j represents the change 
%in the state vector due to reaction j)
S=[1 -1 0  0
    0  0 1 -1];


nt=16000; %simulation time (number of steps)

%initial conditions
t        = 0;
t_vec    = zeros(1,nt);
t_vec(1) = t;
X        = zeros(2,nt);

queue1=[]; % queue for the first delayed reaction
queue2=[]; % queue for the second delayed reaction
for n=1:nt-1
    %birth-death process of two species

    % propensities
    a_v=[lambda1; mu1*X(1,n); ...
        lambda2*X(1,n)^n_Hill/(theta^n_Hill+X(1,n)^n_Hill); mu2*X(2,n)]; 
    a_v_cumsum=cumsum(a_v);
    a_0=a_v_cumsum(end);

    t_wait=(1/a_0)*log(1/rand); %time to next event
    t_prop=t+t_wait;
    % check if a reaction is coming from either of the queues
    queue_check = isempty(queue1)+isempty(queue2);
    switch queue_check
        case 2 %both queues are empty
            t=t_prop;
            t_vec(n+1)=t;
            u=find(a_v_cumsum>=rand*a_0,1,'first');
            if u==1
                delay_1=gamrnd(alpha1,beta1);
                queue1=t+delay_1;
                X(:,n+1)=X(:,n);
            elseif u==2
                X(:,n+1)=X(:,n)+S(:,2); % update X vector
            elseif u==3
                delay_2=gamrnd(alpha2,beta2);
                queue2=t+delay_2;
                X(:,n+1)=X(:,n);
            else
                X(:,n+1)=X(:,n)+S(:,4); % update X vector
            end
        case 1 %one of the queues is empty and the other one is non-empty
            if isempty(queue1)
                if t_prop < queue2(1)
                    t=t_prop;
                    t_vec(n+1)=t;
                    u=find(a_v_cumsum>=rand*a_0,1,'first');
                    if u==1
                        delay_1=gamrnd(alpha1,beta1);
                        queue1=t+delay_1;
                        X(:,n+1)=X(:,n);
                    elseif u==2
                        X(:,n+1)=X(:,n)+S(:,2); % update X vector
                    elseif u==3
                        delay_2=gamrnd(alpha2,beta2);
                        queue2=sort([t+delay_2 queue2]);
                        X(:,n+1)=X(:,n);
                    else
                        X(:,n+1)=X(:,n)+S(:,4); % update X vector
                    end
                else
                    t=queue2(1);
                    t_vec(n+1)=t;
                    X(:,n+1)=X(:,n)+S(:,3);
                    queue2=queue2(2:end);
                end
            else
                if t_prop < queue1(1)
                    t=t_prop;
                    t_vec(n+1)=t;
                    u=find(a_v_cumsum>=rand*a_0,1,'first');
                    if u==1
                        delay_1=gamrnd(alpha1,beta1);
                        queue1=sort([t+delay_1 queue1]);
                        X(:,n+1)=X(:,n);
                    elseif u==2
                        X(:,n+1)=X(:,n)+S(:,2); % update X vector
                    elseif u==3
                        delay_2=gamrnd(alpha2,beta2);
                        queue2=t+delay_2;
                        X(:,n+1)=X(:,n);
                    else
                        X(:,n+1)=X(:,n)+S(:,4); % update X vector
                    end
                else
                    t=queue1(1);
                    t_vec(n+1)=t;
                    X(:,n+1)=X(:,n)+S(:,1);
                    queue1=queue1(2:end);
                end
            end

        case 0 %both queues are non-empty
            [firstInQueue,which] = min([queue1(1) queue2(1)]);
            if t_prop < firstInQueue
                t=t_prop;
                t_vec(n+1)=t;
                u=find(a_v_cumsum>=rand*a_0,1,'first');
                if u==1
                    delay_1=gamrnd(alpha1,beta1);
                    queue1=sort([t+delay_1 queue1]);
                    X(:,n+1)=X(:,n);
                elseif u==2
                    X(:,n+1)=X(:,n)+S(:,2); % update X vector
                elseif u==3
                    delay_2=gamrnd(alpha2,beta2);
                    queue2=sort([t+delay_2 queue2]);
                    X(:,n+1)=X(:,n);
                else
                    X(:,n+1)=X(:,n)+S(:,4); % update X vector
                end
            else
                t=firstInQueue;
                t_vec(n+1)=t;
                if which==1
                    X(:,n+1)=X(:,n)+S(:,1);
                    queue1=queue1(2:end);
                elseif which==2
                    X(:,n+1)=X(:,n)+S(:,3);
                    queue2=queue2(2:end);
                end
            end
    end
end

%
figure(1);
plot(t_vec,X(1,:),'LineWidth',2,'Color',h_X.Color); %plots X
plot(t_vec,X(2,:),'LineWidth',2,'Color',h_Y.Color); %plots Y
xlim([0 t_final])
legend([h_X h_Y],'X','Y')

enter image description here

  • reaction = First@FirstPosition[reactionVec, _?(RandomReal[] < # &)]; this does not work in 2 species but I used alternative one FirstPosition[reaction - RandomReal[], _?Positive][[1]] and worked perfectly. mathematica.stackexchange.com/questions/162666/… – Okkes Dulgerci Apr 28 at 7:45
  • Or one can use rand = RandomReal[]; reaction = First@FirstPosition[reactionVec, _?( rand < # &)]; – Okkes Dulgerci Apr 28 at 7:58

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.