Logistic model

The birth-death model assumes a constant per-capita birth and death rate. But per-capita birth and death rate are not likely constant, they would change depending on population size, e.g., shortage of food, etc.

Here we assume that per-capita birth rate linearly decreases with populatio size n as

	b - b’ n

and per-capita death rate linearly increases with popualtoin size

	d + d’ n
	
where b’ and d’ are parameters controling “density dependency”.

We first introduce a deterministic dynamics given as ODE and then extend it to a stochastic dynamics.

Deterministic model as ODE: pop size n(t) is a real

The logistic model using parameters, b, b1, d, d1, can be arranged in a simpler form that only has two parameters, r and k. r is called intrinsic growth rate given as

	r = b - d

and k as carrying capacity

	k = (b - d)/(b1 + d1)

t=.

deq={n'[t]==r(1-n[t]/k)n[t],n[0]==n0}

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

para={r->1,k->1,n0->0.01}

deqN=deq/.para

solN=NDSolve[deqN,n[t],{t,0,20}]

Plot[Evaluate[n[t]/.solN],{t,0,20},AxesLabel->{"t","n(t)"},PlotRange->All]

LogPlot[Evaluate[n[t]/.solN],{t,0,20},AxesLabel->{"t","n(t)"},PlotRange->All]

This ODE is too easy and Mathematica would not be necessary.

Stochastic model: pop size n(t) is an integer

We now extend the deterministic ODE model as a stochastic model where population size is integer and population size (state variable) jumps from n to n + 1 by birth and jumps from n to n – 1 by death stochastically.

The rule is that time to the next event, an individual gives brith or dies, follows the exponential distirbution with the total rate rateBirth n(t) + rateDeath n(t) where rateBirth and rateDeath is a linear function of n, respectively.

After the inter-event time, we have to choose which event, a birth or a death, has occurred. Here we use the rule that a birth (n --> n+1) and a death (n --> n – 1) occurs with the probability proportional to the total birth rate  and the total death rate, respectively.

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

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

popSize = 10;

Table[-Log[Random[]]/((b-b1 popSize)popSize+(d+d1 popSize)popSize),{10}]

This sequence gives a series of inter-event times (waiting times to the next event). After each inter-event time, population size is incremented by one with the probability (b - b1 popSize)/(b - b1 popSize + d + d1 popSize) and is decreased by one with the rest of the probability. Note that the population size in the numerator and the denominator cancells each other. The stochastic process can be implemented as follows.

Note again that popSize can be zero and division by zero is possible.

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

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

tEnd=200;        (* Set the end of time *)

t=0;  (* Initial condition *)
popSize=10;
seq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,

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

If[popSize>0,
wt=-Log[Random[]]/(rateBirthTotal+rateDeathTotal),
wt=999;
];

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

seq;

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

Because this is a stochastic process, each realization results in different dynamics. 
To view the stochastic dynamics, it is often convenient to observe the populatoin size at a fixed time interval, say DT = 1.0. To do this, we introduce a time tEq that always advances by a constant DT and record the population size at each tEq.

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

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;

If[popSize>0,
wt=-Log[Random[]]/(rateBirthTotal+rateDeathTotal),
wt=999;
];

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

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

seq;

seqEq

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

If we run the process many times with the same conditions, we may have a hint to understand the stochastic process.Let' s run the simulation multiple times.

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

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

seqSims={};  (* Sequence that records the population size in each realization *)
seqSimsEq={}; (* Sequence that records the population size measured at tEq in each realization *)

Do[
t=0;
tEq=0;
popSize=10;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,

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

If[popSize>0,
wt=-Log[Random[]]/(rateBirthTotal+rateDeathTotal),
wt=999;
];

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

t+=wt;

If[Random[]<rateBirthTotal/(rateBirthTotal+rateDeathTotal),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];

AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSimsEq;

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

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

Take out the time sequence

seqSimsEq[[1]]

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

Take out the sequences of the population size

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

Transpose[seqPopSize];

Calculate the average (or mean) of the populatoin size measured at tEq = 0, 1, 2, 3, ...

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]]

Combine the time sequence with the sequence of the averaged pouplation size

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

gSimAverage=ListPlot[seqTPopSizeAverage]

gSimAverageLog=ListPlot[seqTPopSizeAverageLog,PlotRange->All]

Comparison with ODE solution

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

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

deqN=deq/.{b->0.2,d->0.1,b1->0,d1->0.001}/.n0->10

tEnd=200

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

gSolODE=Plot[Evaluate[n[t]/.sol],{t,0,tEnd},PlotRange->All]

Show[gSolODE,gSimAverage]

How does the averaged population size change with time? The averaged population size converges to an equilibrium that is close to the carrying capacity k defined in the ODE logistic model.

However, it can be shown that the averaged populatoin size of the stocashtic logistic model does NOT match the solution of the ODE logistic model. In the ODE logistic model, population size converges to the carrying capacity k. But in stochastic logistic model, the averaged popualtion size converges to the value less than the carrying capacity. The difference depends on the density-dependency b1 and d1. 

Furthermore, it can be shown that the stochastic logistic model ALWAYS results in population extinction after sufficient time has passed. This is a bit counter-intuitive because we have observed that the averaged populatoin size converges to a positive equilibrium.

For analytical proof, see the documents at http://gi.ics.nara-wu.ac.jp/~takasu/lecture/graduate19.html

Time to extinction in the stochastic logistic model

Let' s run a single realization until the population is gone extinct.

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

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

tEnd=10000;        (* 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[popSize>0,  (* <---  We run the simulation while popSize is positive *)

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

If[popSize>0,
wt=-Log[Random[]]/(rateBirthTotal+rateDeathTotal),
wt=999;
];

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

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

t

Last[seq]

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

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

Let' s run the simulation many times and see the distribution of the times to extinction.

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

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

tEnd=10000;        (* Set the end of time *)

DT=1;     (* Constant time interval *)

numRealizations=100;  (* Number of realizations *)

seqTimeExtinct={};

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[popSize>0,  (* <---  We run the simulation while popSize is positive *)

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

If[popSize>0,
wt=-Log[Random[]]/(rateBirthTotal+rateDeathTotal),
wt=999;
];

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

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

PrintTemporary[t];
AppendTo[seqTimeExtinct,t];
,{i,numRealizations}
]

seqTimeExtinct

Histogram[seqTimeExtinct]

Mean[seqTimeExtinct]