3.3: Simulation of Population Growth
As we have seen, stochastic modeling is significantly more complicated than deterministic modeling. As the modeling becomes more sophisticated, a numerical simulation becomes necessary. Here, for illustration, we show how to simulate individual realizations of population growth.
A naive approach would make use of the birth rate \(b\) directly. During the short time interval \(\Delta t\) , each individual has probability \(b \Delta t\) of giving birth. We can decide if an individual gives birth by generating a random deviate (a pseudo-random number between zero and one): if the random deviate is less than \(b \Delta t\) , then the individual gives birth; if larger than \(b \Delta t\) , then the individual does not. With \(N\) individuals at time \(t\) , we then simply compute \(N\) random deviates. Counting the number of random deviates less than \(b \Delta t\) allows us to update the population size to the time \(t+\Delta t\) . For accuracy, \(\Delta t\) must be small, making this a computationally slow method.
There is, however, a much more efficient way to simulate population growth. Define a random variable \(\tau=\tau(N)\) to be the time it takes for a population to grow from size \(N\) to size \(N+1\) because of a single birth. The random variable \(\tau\) is called the interevent time and represents the elapsed time between births. A simulation from population size \(N_{0}\) to size \(N_{f}\) would then simply require computing \(N_{f}-N_{0}\) different random values of \(\tau\) , a relatively easy and quick computation if we know the probability density function (pdf) of \(\tau\) .
Accordingly, we define \(P(\tau)\) to be the pdf of \(\tau\) for a population of size \(N\) . The cumulative distribution function (cdf), \(F(\tau)\) , defined as the probability that the interevent time is less that \(\tau\) is given by
\[F(\tau)=\int_{0}^{\tau} P(\tau) d \tau \nonumber \]
where \(P(\tau)=F^{\prime}(\tau)\) . The complementary cumulative distribution function (ccdf), \(G(\tau)\) , defined as the probability that the interevent time is greater than \(\tau\) is given by \(G(\tau)=1-F(\tau)\)
Now the probability that the interevent time is greater than \(\tau+\Delta \tau\) , with \(\Delta \tau\) small, is given by the probability that it is greater than \(\tau\) times the probability that there are no births in the time interval \(\Delta \tau\) . Therefore, \(G(\tau+\Delta \tau)\) satisfies
\[G(\tau+\Delta \tau)=G(\tau)(1-b N \Delta \tau) \nonumber \]
Differencing \(G\) and taking the limit \(\Delta \tau \rightarrow 0\) yields the differential equation
\[\frac{d G}{d \tau}=-b N G \nonumber \]
which can be integrated using the initial condition \(G(0)=1\) to obtain
\[G(\tau)=e^{-b N \tau} \nonumber \]
From \(G(\tau)\) , we can find
\[F(\tau)=1-e^{-b N \tau}, \quad P(\tau)=b N e^{-b N \tau} \nonumber \]
The pdf, \(P(\tau)\) , has the form of an exponential distribution with parameter \(b N\) .
Here, we make use of a well-known result from probability theory that enables us to compute \(\tau\) using random deviates. With \(y\) a random deviate, \(\tau\) can be computed from \(\tau=F^{-1}(y)\) , where \(F^{-1}\) is the inverse function of \(F .\) The correct formula for the exponential distribution is
\[\tau=-\frac{\ln (1-y)}{b N} \nonumber \]
To simulate a population growing from \(N_{0}\) to \(N_{f}\) , we compute \(N_{f}-N_{0}\) random deviates \(y\) , and then compute the corresponding interevent times using \((3.3.6)\) , taking care to adjust the population size \(N\) as the population grows.
Below, we illustrate a simple MATLAB function that simulates one realization of population growth from initial size \(N_{0}\) to final size \(N_{f}\) , with birth rate \(b\) .
The function population_growth_simulation.m can be driven by a MATLAB script to compute realizations of population growth. For instance, the following script computes 25 realizations for a population growth from 10 to 100 with \(b=1\) and plots all the realizations: