Birth Death Process with Delay

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 = λ #[];
a2 = μ #[];
a0 = Total@{a1, a2};
reaction = 1/a0 Accumulate@{a1, a2};
pos = First@FirstPosition[reaction, _?(RandomReal[] < # &)];
Δt =
RandomVariate@ExponentialDistribution[a1 + a2];

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

]

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

ListLinePlot[sim, Epilog -> {Red, PointSize[Medium], Point[sim]},
Frame -> True, PlotTheme -> "Detailed",
FrameLabel -> {"Time", "Population"}, ImageSize -> Large,
InterpolationOrder -> 0] 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 = λ #[];
a2 = μ #[];
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, {#[] + t, #[] - 1},
pos == 1,

Which[

t + td < t +
tw, {#[] + t + td, #[] + 1},

t + td >= t + tw,

Which[

pos2 == 2, {#[] +t + tw, #[] - 1},

pos2 == 1, {#[] + t + tw, #[]}

]

]

]

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

ListLinePlot[sim, Epilog -> {Red, PointSize[Medium], Point[sim]},
Frame -> True, PlotTheme -> "Detailed",
FrameLabel -> {"Time", "Population"}, ImageSize -> Large,
InterpolationOrder -> 0] 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 #[]];

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

) &, {0, initialPop}, 10]];
ListLinePlot[sim, Epilog -> {Red, PointSize[Medium], Point[sim]},
Frame -> True, PlotTheme -> "Detailed",
FrameLabel -> {"Time", "Population"}, ImageSize -> Large,
InterpolationOrder -> 0] 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 []'s in your definition of a1 and a2 should be []'s. – Chris K Feb 3 '18 at 3:00
• Yes they should be [], I don't know how that happened. I'll fixed. – OkkesDulgerci Feb 3 '18 at 6:02

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}]] 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]] • reaction = First@FirstPosition[reactionVec, _?(RandomReal[] < # &)]; this does not work in 2 species but I used alternative one FirstPosition[reaction - RandomReal[], _?Positive][] and worked perfectly. mathematica.stackexchange.com/questions/162666/… – OkkesDulgerci Apr 28 '18 at 7:32

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 =
tDelayY =

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"][]},
PlotLegends -> {"Y"}, Frame -> True, PlotTheme -> "Detailed",
FrameLabel -> {"Time", "Population"}, ImageSize -> Large,
InterpolationOrder -> 0],
ListLinePlot[simX, PlotLegends -> {"X"},
PlotStyle -> {ColorData[97, "ColorList"][]},
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}]] 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') • reaction = First@FirstPosition[reactionVec, _?(RandomReal[] < # &)]; this does not work in 2 species but I used alternative one FirstPosition[reaction - RandomReal[], _?Positive][] and worked perfectly. mathematica.stackexchange.com/questions/162666/… – OkkesDulgerci Apr 28 '18 at 7:45
• Or one can use rand = RandomReal[]; reaction = First@FirstPosition[reactionVec, _?( rand < # &)]; – OkkesDulgerci Apr 28 '18 at 7:58