For the deterministic model, including a death rate just changed $\lambda$ to $\lambda-\mu$. Rather than adding a death rate of $\mu=0.05$, we would have obtained the same deterministic model if we had simply changed $\lambda$ from $0.1$ to $0.05$ and didn't add any deaths.
The stochastic model, however, treats deaths differently from simply decreasing births. Before the introduction of deaths, in a small time window of width $\Delta t$, the population size $N$ could only jump up by one, from the number $n$ to the number $n+1$, or stay at the same population size of $n$. With the introduction of deaths, there is now the possibility that the population decreases by one to $n-1$.
The probability that a bacterium divides in time interval $\Delta t$ and gives birth to a new cell is unchanged from before the introduction of the death rate. Using the value $\lambda=0.1$, this birth probability is
$P(N(t+\Delta t) = n+1 \,|\, N(t) = n) = $
.
The probability that a bacterium dies in the time interval $\Delta t$ is identical in form to the birth probability. In this case $N$ decreases from the number $n$ to the number $n-1$, and the rate depends on $\mu$ rather than $\lambda$. Using the value $\mu=0.05$, this death probability is
$P(N(t+\Delta t) = n-1 \,|\, N(t) = n) = $
.
Since we view the time interval as being small, we allow at most one birth or one death in the interval. It turns out, we also don't have to consider the possibility that we have both one birth and one death in the interval. (Just like the event of two births or two deaths, the probability of the event of one birth and one death is vanishingly small as $\Delta t$ gets small.) Therefore, the probability that $N$ stays at the value $n$ during the whole interval is just one minus the sum of the above two probabilities:
$P(N(t+\Delta t) = n \,|\, N(t) = n) = $
Again, we don't allow $N$ to change by more than one during the interval, so the probability of any other change in $N$ is
$P(N(t+\Delta t) = m \,|\, N(t) = n) = $
for any $m$ not equal to $n$, $n+1$, or $n-1$.
The stochastic model can be written more compactly as
$$P(N(t+\Delta t) = m \,|\, N(t)=n) = \begin{cases}
_ & \text{if $m=n+1$}\\
_ &\text{if $m=n$}\\
_ & \text{if $m=n-1$}\\
_ & \text{otherwise.}
\end{cases}$$
In terms of the parameters $\lambda$ and $\mu$, the probabilistic model is:
$P(N(t+\Delta t) = n+1 \,|\, N(t) = n) = $
$P(N(t+\Delta t) = n \,|\, N(t) = n) = $
$P(N(t+\Delta t) = n-1 \,|\, N(t) = n) = $
$P(N(t+\Delta t) = m \,|\, N(t) = n) = $
for any $m$ not equal to $n$, $n+1$, or $n-1$.
The stochastic model can be written more compactly as
$$P(N(t+\Delta t) = m \,|\, N(t)=n) = \begin{cases}
_ & \text{if $m=n+1$}\\
_ &\text{if $m=n$}\\
_ & \text{if $m=n-1$}\\
_ & \text{otherwise.}
\end{cases}$$
We refer to such a stochastic model that includes births and deaths as a birth-death model.
Now that we have two possible events to track, births and deaths, simulating the stochastic system is slightly more complicated. The simulation is simplified by first ignoring any difference between births and deaths so that the model appears similar to the birth-only model. If we don't care about distinguishing births from deaths but just want to know when the next event occurs (where an event is either a birth or death), then we just need to know the combined rate of these events. If the population size is $N$, the total rate of births is
and the total rate of deaths is
. The combined rate of birth or death events in their sum:
.
The first step in the stochastic simulation is to determine when the next event occurs. If we let $T$ be the waiting time until the next event, then $T$ will be an
random variable with mean that is equal to
.
Only after we draw a random waiting time $T$, do we worry about the identity of the event. The likelihood of the event being a birth or death depends on the relative magnitudes of the birth and death rates. If these rates were equal, then the event would be a birth with probability
and the event would be a death with probability
. However, since in our model the per capita birth rates is $\lambda=0.1$ per minute and per capita death rates is $\mu=0.05$ per minute, what's the likelihood the event was a birth?
What's the likelihood that the event was a death?
(Since one of the two must have occurred, this probabilities must sum to 1. Simply rescale the values of $\lambda$ and $\mu$ so that they sum to 1.) If rounding to a decimal, include at least five significant digits in your response.
In the simulation, we want to be able to determine these probabilities even as you change $\lambda$ and $\mu$, so we want expressions in terms of the parameters $\lambda$ and $\mu$. To rescale the quantities $\lambda$ and $\mu$ so that they sum to 1, you should divide by
. Therefore, in terms of the parameters $\mu$ and $\lambda$, what is the probability that the event was a birth?
What is the probability that the event was a death?
.
To summarize, we will simulate each event of the birth-death process with two steps. First, we'll draw an _ random variable at the combined rate of $_$, giving the waiting time $T$ until the next event. Second, we'll flip a coin that will come out Heads with probability $_$. If we get Heads, the event is a birth and we increase the bacteria count $N$ by one ; otherwise the event is a death, and we decrease $N$ by one.
We have one more wrinkle to address given that we now have deaths. What if we have one bacterium left and the next event is a death? At that point, how bacteria are there in the colony?
Given that population size, what will be the overall rate of births?
What will be the overall rate of deaths?
. We refer to this event as an extinction event. Once an extinction occurs, we stop the simulation, as nothing else will occur.
Can the phenomenon of extinction occur in the deterministic model? If the initial condition $N(0)$ is larger than zero, will the solution of the deterministic equation ever reach zero?
These extinctions are a new feature of the stochastic model.
To simulate the birth-death process, you can start with the incomplete birth-death process R script. The script needs the information you determined above: the overall event rate, the probability that an event is a birth, and the solution to the deterministic solution. See comments at beginning of file to instruction on how to edit the file to make it run. Fill in those values and also set the parameters such as lambda and mu to the values corresponding to this bacteria growth problem. Then, you can source the file to simulate the birth-death process and plot the results.
Set tmax to about 40 and n_samples to about 10. With the original parameter values $\lambda=0.1$ and $\mu=0.05$, the deterministic model shows exponential
. With the initial condition $N(0)=2$, do the stochastic model solutions always exhibit exponential _?
Even though the growth rate is larger than the death rate, you will often observe
, i.e., events where the population completely dies out (although with only 10 samples, you may not see them every simulation).