Extension of an ODE model to a stochastic model

Any ODE model can be readily extended as a "stochastic" model in which population size is integer-valued (non-zero) and it changes stochastically with a certain rule (increase by immigration and birth, and decrease by emigration and death, etc.).

Try to extend these ODE models to stochastic versions and see how they behave!

Constant birth and death + Constant immigration and emigration

t=.

deq={n'[t]==b n[t]-d n[t]+flow,n[0]==n0}

sol=DSolve[deq,n[t],t]

para={b->0.1,d->0.15,flow->1}

solN=sol/.para/.n0->0.1

Plot[n[t]/.solN,{t,0,100},PlotRange->All]

How can the ODE model be extended as a stochastic pop dynamics?

rateBirth = 0.1; (* set the per-capita birth rate *)
rateDeath=0.15; (* set the per-capita death rate *)
rateInflow=1;

tEnd=100;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=1;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,

rateTotal=(rateBirth+rateDeath)popSize+rateInflow;

If[popSize>0,
wt=-Log[Random[]]/rateTotal,
wt=999;
];

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

probIncrease=(rateBirth popSize+rateInflow)/rateTotal;

t+=wt;
If[Random[]<probIncrease,popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];
AppendTo[seq,{t,popSize}]
]

ListPlot[seqEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

rateBirth = 0.1; (* set the per-capita birth rate *)
rateDeath=0.15; (* set the per-capita death rate *)
rateInflow=1;

tEnd=100;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=10;  (* Number of realizations *)
seqSimsEq={};

Do[

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=1;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,

rateTotal=(rateBirth+rateDeath)popSize+rateInflow;

If[popSize>0,
wt=-Log[Random[]]/rateTotal,
wt=999;
];

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

probIncrease=(rateBirth popSize+rateInflow)/rateTotal;

t+=wt;
If[Random[]<probIncrease,popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];
AppendTo[seq,{t,popSize}]
];

AppendTo[seqSimsEq,seqEq],{numRealizations}

]

gSims=ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]];

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq];

Display the population sizes at t = tEnd

Last[Transpose[seqPopSize]]

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]];

seqTPopSizeAverage=Transpose[{seqT,seqPopSizeAverage}];
seqTPopSizeAverageLog=Transpose[{seqT,Log[seqPopSizeAverage]}];

gSimAverage=ListPlot[seqTPopSizeAverage,Joined->True,PlotStyle->{Thickness[0.01],RGBColor[1,0,0]}]

Show[gSims,gSimAverage]

Logistic growth + Constant immigration and emigration

t=.
b=.
b1=.
d=.
d1=.

deq={n'[t]==(b-b1 n[t])n[t]-(d+d1 n[t])n[t]+flow,n[0]==n0}

para={b->0.2,b1->0,d->0.1,d1->0.001,flow->10}

deqN=deq/.para/.n0->1

sol=NDSolve[deqN,n[t],{t,0,100}]

Plot[n[t]/.sol,{t,0,100},PlotRange->All]

How can the ODE model be extended as a stochastic pop dynamics?

b = 0.2; (* set the per-capita birth rate *)
d=0.1; (* set the per-capita death rate *)
b1=0.0;
d1=0.001;
rateInflow=10;

Print["r=",b-d];
Print["k=",(b-d)/(b1+d1)];

tEnd=200;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=10;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,

rateBirthTotal=(b-b1 popSize)popSize;
rateDeathTotal=(d+d1 popSize)popSize;

rateTotal=rateBirthTotal+rateDeathTotal+rateInflow;

If[popSize>0,
wt=-Log[Random[]]/rateTotal,
wt=999;
];

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
If[Random[]<(rateBirthTotal+rateInflow)/rateTotal,popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];
AppendTo[seq,{t,popSize}]
]

ListPlot[seqEq,Joined->True,PlotRange->All,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

b = 0.2; (* set the per-capita birth rate *)
d=0.1; (* set the per-capita death rate *)
b1=0.0;
d1=0.001;
rateInflow=10;

Print["r=",b-d];
Print["k=",(b-d)/(b1+d1)];

tEnd=200;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=10;  (* Number of realizations *)
seqSimsEq={}; 

Do[

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=10;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,

rateBirthTotal=(b-b1 popSize)popSize;
rateDeathTotal=(d+d1 popSize)popSize;

rateTotal=rateBirthTotal+rateDeathTotal+rateInflow;

If[popSize>0,
wt=-Log[Random[]]/rateTotal,
wt=999;
];

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
If[Random[]<(rateBirthTotal+rateInflow)/rateTotal,popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];
AppendTo[seq,{t,popSize}]
]

AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

gSims=ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"},PlotRange->All]

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]];

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq];

Display the population sizes at t = tEnd

Last[Transpose[seqPopSize]]

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]];

seqTPopSizeAverage=Transpose[{seqT,seqPopSizeAverage}];
seqTPopSizeAverageLog=Transpose[{seqT,Log[seqPopSizeAverage]}];

gSimAverage=ListPlot[seqTPopSizeAverage,Joined->True,PlotStyle->{Thickness[0.01],RGBColor[1,0,0]},PlotRange->All]

Show[gSims,gSimAverage]

Epidemic SIS model

We assume two popualtions, Susceptible and Infectious. Let S and I denote the population size of each population, respectively.

Assumption is that S and I contact randamly by the law of mass action, and infection occurs with a constant infection rate when S meets I.

I recoveres with a constant recovery rate to be S again. 

Because the total population size remains the same (s(t) + i(t): const), this dynamics is essenntially the same as the logistic model.

t=.

deq={s'[t]==-beta s[t]i[t]+gamma i[t],i'[t]==beta s[t]i[t]-gamma i[t],s[0]==s0,i[0]==i0}

para={beta->0.1,gamma->1}

deqN=deq/.para/.{s0->99,i0->1}

sol=NDSolve[deqN,{s[t],i[t]},{t,0,10}]

Plot[Evaluate[{s[t],i[t]}/.sol],{t,0,10},AxesLabel->{"t","s(t),i(t)"}]

How can the ODE model be extended as a stochastic pop dynamics?

Epidemic SIR model

We assume three popualtions, Susceptible, Infectious and Removed with lifetime immunity acqured after recovery from I. Let S, I and R denote the population size of each population, respectively.

Assumption is that S and I contact randamly by the law of mass action, and infection occurs with a constant infection rate when S meets I.

I recoveres with a constant recovery rate to be R.

t=.

deq={s'[t]==-beta s[t]i[t],i'[t]==beta s[t]i[t]-gamma i[t],r'[t]==gamma i[t],s[0]==s0,i[0]==i0,r[0]==r0}

para={beta->0.1,gamma->1}

deqN=deq/.para/.{s0->99,i0->1,r0->0}

sol=NDSolve[deqN,{s[t],i[t],r[t]},{t,0,10}]

Plot[Evaluate[{s[t],i[t],r[t]}/.sol],{t,0,10},AxesLabel->{"t","s(t),i(t),r(t)"}]

How can the ODE model be extended as a stochastic pop dynamics?

Lotka Volterra competition model

Two populations compete with each other.

t=.
b1=.
b11=.
b21=.
b22=.
d1=.
d11=.
d2=.
d21=.

deq={n1'[t]==(b1-b11 n1[t]-b12 n2[t])n1[t]-(d1+d11 n1[t]+d12 n2[t])n1[t],
n2'[t]==(b2-b21 n1[t]-b22 n2[t])n2[t]-(d2+d21 n1[t]+d22 n2[t])n2[t],
n1[0]==n10,n2[0]==n20}

para={b1->0.2,b11->0,b12->0.0001,d1->0.1,d11->0.001,d12->0,
b2->0.2,b21->0,b22->0.001,d2->0.1,d21->0.0002,d22->0}

deqN=deq/.para/.{n10->1,n20->2}

tEnd=200;
sol=NDSolve[deqN,{n1[t],n2[t]},{t,0,tEnd}]

Plot[Evaluate[{n1[t],n2[t]}/.sol],{t,0,tEnd},AxesLabel->{"t","n1(t), n2(t)"}]

ParametricPlot[Evaluate[{n1[t],n2[t]}/.sol],{t,0,tEnd},AxesLabel->{"n1(t)", "n2(t)"},PlotRange->All]

How can the ODE model be extended as a stochastic pop dynamics?