Birth-death model

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

Birth model

We assume a population in which each individuals gives birth to new individual with a constant per-capita birth rate. Death of individuals is ignored.

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

t=.

deq={n'[t]==rateBirth n[t],n[0]==n0}

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

para={rateBirth->0.1,n0->10}

deqN=deq/.para

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

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

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

This ODE is too easy and Mathematica would not be necessary. The solution is that the population size n(t) exponentiall increases with time t starting from an initial value n0. In log scale, it linearly increases with time t.

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 stochastically.

The rule is that time to the next event (an individual joins the population by birth) follows the exponential distirbution with the parameter rateBirth*n given in the right hand side of the ODE.

rateBirth=0.1;
popSize=10;

Table[-Log[Random[]]/(rateBirth 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. The stochastic process can be implemented as follows.

rateBirth = 0.1;  (* Set the per-capita birth rate *)
tEnd=30;        (* Set the end of time *)

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

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateBirth popSize);
t+=wt;
popSize++;
AppendTo[seq,{t,popSize}]
]

seq

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

Note that the population size increases by one at each time a birth occurs. Note that inter-event time becomes smaller and smaller as the population size increases. 

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.

rateBirth = 0.1;  (* Set the per-capit birth rate *)
tEnd=30;        (* 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,
wt=-Log[Random[]]/(rateBirth popSize);

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

t+=wt;
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.

rateBirth = 0.1;  (* Set the per-capita birth rate *)
tEnd=50;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* 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=1;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateBirth popSize);

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

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

seqSims

seqSimsEq

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

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

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

Take out the time sequence

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]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeAverageLog,{1,t},t]

gFit=Plot[fit,{t,0,30}]

Show[gSimAverageLog,gFit]

gSol=Plot[Exp[0.1 t],{t,0,30}]

Show[gSol,gSimAverage]

How does the averaged population size change with time? Increase the number of realizations and check the result. 

It can be shown that the averaged populatoin size is given as n0 exp(rateBirth t).

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

Death model

We assume a population in which each individual dies with a constant per-capita death rate. Birth of individuals is ignored and hence extinction (population size becomes zero) is inevitable.

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

t=.

deq={n'[t]==-rateDeath n[t],n[0]==n0}

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

para={rateDeath->0.1,n0->100}

deqN=deq/.para

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

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

This ODE is too easy and Mathematica would not be necessary. The solution is that the population size n(t) linearly increases with time t starting from an initial value n0.

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 stochastically by death.

The rule is that time to the next event (an individual dies) follows the exponential distirbution with the total death rate given in the right hand side of ODE.

rateDeath=0.1;
popSize=10;

Table[-Log[Random[]]/(rateDeath 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 decremented by one. The stochastic process can be implemented as follows.

Note that division by zero (wt = infinitly large) occurs when popSize becomes zero.

rateDeath = 0.1;  (* Set the per-capita death rate *)
tEnd=100;        (* Set the end of time *)

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

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateDeath popSize);
t+=wt;
popSize--;
AppendTo[seq,{t,popSize}]
]

seq

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

Note that the population size decreases by one at each time an individual dies. 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.

rateDeath = 0.1;  (* Set the per-capita death rate *)
tEnd=30;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

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

AppendTo[seq,{t,popSize}];

While[t<tEnd,
(* We check if popSize is positive or zero *)
If[popSize>0,
wt=-Log[Random[]]/(rateDeath popSize),
wt=999;
];

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

t+=wt;
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.

rateDeath = 0.1;  (* Set the per-capita death rate *)
tEnd=30;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* 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=100;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,

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

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

t+=wt;
popSize--;
AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSims

seqSimsEq

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

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

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

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]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeAverageLog,{1,t},t]

gFit=Plot[fit,{t,0,30}]

Show[gSimAverageLog,gFit]

gSol=Plot[100Exp[-0.1 t],{t,0,30}]

Show[gSol,gSimAverage]

How does the averaged population size change with time? 

It can be shown that the averaged populatoin size is given as n0 exp(-rateDeath t).

It should be noted that in this death model, population size eventually becomes zero when the last individual dies. Then the inter-event time to the next death is infinite because of division by zero. This means that once the population goes extinct (population size becomes zero), the population remains empty forever.

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

Birth-death model

We assume a population in which each individual gives birth and dies with constant per-capita birth and death rate, respectively.
Extinction is possible even if average populatoin size exponentially increases.

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

rateBirth=.
rateDeath=.

deq={n'[t]==(rateBirth-rateDeath)n[t],n[0]==n0}

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

para={rateBirth->0.2,rateDeath->0.1,n0->10}

deqN=deq/.para

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

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

This ODE is too easy and Mathematica would not be necessary. The solution is that the population size n(t) linearly increases with time t starting from an initial value n0.

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

And 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.

rateBirth = 0.2; (* set the per-capita birth rate *)
rateDeath=0.1; (* set the per-capita death rate *)
popSize = 10;

Table[-Log[Random[]]/((rateBirth+rateDeath)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 rateBirth/(rateBirth + rateDeath) 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.

rateBirth = 0.2; (* set the per-capita birth rate *)
rateDeath=0.1; (* set the per-capita death rate *)
tEnd=30;        (* Set the end of time *)

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

AppendTo[seq,{t,popSize}];

While[t<tEnd,

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

t+=wt;
If[Random[]<rateBirth/(rateBirth+rateDeath),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.

rateBirth = 0.2; (* set the per-capita birth rate *)
rateDeath=0.1; (* set the per-capita death rate *)
tEnd=30;        (* 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,

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

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

t+=wt;
If[Random[]<rateBirth/(rateBirth+rateDeath),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.

rateBirth = 0.2; (* set the per-capita birth rate *)
rateDeath=0.1; (* set the per-capita death rate *)
tEnd=30;        (* 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,

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

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

t+=wt;

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

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

seqSims

seqSimsEq

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

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

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

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

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]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeAverageLog,{1,t},t]

gFit=Plot[fit,{t,0,30}]

Show[gSimAverageLog,gFit]

gSol=Plot[10Exp[0.1 t],{t,0,30}]

Show[gSol,gSimAverage]

How does the averaged population size change with time? 

It can be shown that the averaged populatoin size is given as n0 exp( (rateBirth – rateDeath) t).

It should be noted that in this birth-death model, population size can become zero when the last individual dies. Then the inter-event time to the next event (either birth or death) is infinite because of division by zero. This means that once the population goes extinct (population size becomes zero), the population remains empty forever.

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

Extinction in the birth-death model: pop size n(t) is an integer

If the initial population size is small in size, the population can easily go extinct. Run the simulation with a small number of initial population size.

When population goes extinct, inter-event time to the next event (birth or death) is infinitely large and no events can happen (division by zero).

rateBirth = 0.2; (* set the per-capita birth rate *)
rateDeath=0.1; (* set the per-capita death rate *)
tEnd=30;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* 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=1; (* Initial popSize *)
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,

(* We check popSize before we calculate wt *)
If[popSize>0,
wt=-Log[Random[]]/((rateBirth+rateDeath)popSize),
wt=999;
];

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

t+=wt;

If[Random[]<rateBirth/(rateBirth+rateDeath),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"}]

Take out the sequences of the population size

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

Transpose[seqPopSize]

popSizeLast=Last[Transpose[seqPopSize]]

Count[popSizeLast,0]

numRealizations

It can be shown that the probability of population extinction P0 is analytically given as

	P0 = (rateDeath/rateBirth) ^ n0	when rateBirth > rateDeath
	P0 = 1					when rateBirth < rateDeath

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

This is called “demographic stchasticity” that a small-sized population has a higher risk of extinction even if per-capita brith rate is larger than the per-capita death rate.