Derivation of extinction coefficient
Stochastic extinction is most relevant during the initial stage of infection. At this stage the number of target cells is small. We can consider that the target cells are constant or equal to the initial number of target cells (T ≈ T0). As a result, the states become decoupled making the stochastic events independent of each other. In addition, each event produces progeny over a lifetime which is also independent of the lifetime of all other events. More details on how to derive a branching process from a CTMC can be found in [43]. Under these conditions, the CTMC model becomes a multi-type branching process where the reduced state vectors now represent \(\vec {m}\)=\((n_{E_{1}}\), \(n_{I_{1}}, n_{V_{1}}\), \(n_{E_{2}}, n_{I_{2}}, n_{V_{2}})\), where \(n_{E_{1}}\) and \(n_{E_{2}}\) are the numbers of eclipse cells, \(n_{I_{1}}\) and \(n_{I_{2}}\) are the infected cells, and \(n_{V_{1}}\) and \(n_{V_{2}}\) are the virions of both viruses. Including the assumption of a constant number of target cells, the reduced model is
$$\begin{array}{*{20}l} V_{1}\xrightarrow {\beta_{1} T} E_{1} && V_{2}\xrightarrow {\beta_{2} T} E_{2}\\ E_{1}\xrightarrow {k_{1}} I_{1} && E_{2}\xrightarrow {k_{2}} I_{2}\\ I_{1}\xrightarrow {p_{1}} V_{1} && I_{2}\xrightarrow {p_{2}} V_{2}\\ I_{1}\xrightarrow {\delta_{1}} \emptyset && I_{2}\xrightarrow {\delta_{2}} \emptyset \\ V_{1}\xrightarrow {c_{1}} \emptyset && V_{2}\xrightarrow {c_{2}} \emptyset. \\ \end{array} $$
Thus the continuous-time Markov chain becomes a multi-type branching process that describes the dynamics of a population of individuals having birth and death independently according to the specified (in this case exponential) probability mass function. If a time-homogeneous CTMC is a branching process, the only absorbing state is \(\vec {0}\). For this model we defined the absorbing state as \(\vec {0}\) and the probability to reach this state from, say, \(\vec {m}\), is \(\xi (\vec {m})\). This probability is referred to as the extinction probability. Biologically, the extinction probability is defined as the probability that the two types of viruses and all their infected cells are completely eliminated from the host. Once a transition occurs, the current state \(\vec {m}\) is incremented by one of the transition vectors given below.
$$\begin{array}{*{20}l} d\vec{m_{1}}=(0, +1, 0, 0, 0,0)\ & \text{for}\ V_{1}\xrightarrow {\beta_{1} T} E_{1}\\ d\vec{m_{2}}=(0, -1, +1, 0, 0, 0)\ & \text{for}\ E_{1}\xrightarrow {k_{1}} I_{1} \\ d\vec{m_{3}}=(+1, 0, 0, 0, 0, 0) \ &\text{for} \ I_{1}\xrightarrow {p_{1}} V_{1} \\ d\vec{m_{4}}=(0, 0, -1, 0, 0, 0) \ & \text{for}\ I_{1}\xrightarrow {\delta_{1}} \phi \\ d\vec{m_{5}}=(-1, 0, 0, 0, 0, 0)\ & \text{for}\ V_{1}\xrightarrow {c_{1}} \phi \\ d\vec{m_{6}}=(0, 0, 0, 0, +1, 0)\ & \text{for}\ V_{2}\xrightarrow {\beta_{2} T} E_{2} \\ d\vec{m_{7}}=(0, 0, 0,0, -1, +1)\ & \text{for}\ E_{2}\xrightarrow {k_{2}} I_{2} \\ d\vec{m_{8}}=(0, 0, 0, +1, 0, 0)\ &\text{for} \ I_{2}\xrightarrow {p_{2}} V_{2} \\ d\vec{m_{9}}=(0, 0, 0, 0, 0, -1)\ & \text{for}\ I_{2}\xrightarrow {\delta_{2}} \phi \\ d\vec{m_{10}}=(0, 0, 0, -1, 0, 0)\ &\text{for} \ V_{2}\xrightarrow {c_{2}} \phi. \end{array} $$
If the rate of the ith reaction is defined as ai such that a1=β1TV1,a2=β2TV2, a3=k1E1,a4=k2E2, a5=δ1I1,a6=δ2I2, a7=p1I1,a8=p2I2, a9=c1V1,a10=c2V2, then the probability that the ith reaction is the next reaction is given by
$$\begin{array}{*{20}l} P_{i}(\vec{m})&=\frac{a_{i}(\vec{m})}{Z(\vec{m})}\\ \text{where}\ Z(\vec{m})&=\sum_{i}^{n_{max}} a_{i}(\vec{m}), \end{array} $$
and nmax is the number of transitions involved in the model and is equal to 10. The time of the next reaction is a random variable with distribution \(Z(\vec {m}) \exp (-Z(\vec {m})t)\) with mean \(\frac {1}{Z(\vec {m})}\) (according to the Gillespie algorithm). The probability that a simultaneous exposure to both viruses eventually evolves to extinction, or reaches the absorbing state, (0,0,0,0,0,0), from state \(\vec {m}\) or the extinction coefficient, \(\xi (\vec {m})\), is
$$\begin{array}{*{20}l} \xi(\vec{m})&=\sum_{i}P_{i}(\vec{m})\xi(\vec{m}+d\vec{m_{i}}), \vec{m}\neq \vec{0},\\ \xi(\vec{0})&=1\ \text{when} \ \vec{m}= \vec{0}\text{.}\notag \end{array} $$
(1)
Substituting the expressions for \(P_{i}(\vec {m})\) and \(\xi (\vec {m}+d\vec {m_{i}})\) in Eq. (1), the extinction coefficient becomes:
$$\begin{array}{*{20}l} \xi(\vec{m})&=\frac{\beta_{1}T{V_{1}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}+1}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}\\&\quad+\frac{k_{1}{E_{1}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}-1}_{E_{1}}\rho^{n_{I_{1}}+1}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}\notag \\&+\frac{p_{1}{I_{1}}}{Z}\rho^{n_{V_{1}}+1}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}\\&\quad+\frac{\delta_{1}{I_{1}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}-1}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}\notag \\&+\frac{c_{1}{V_{1}}}{Z}\rho^{n_{V_{1}}-1}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}\\&\quad+ \frac{\beta_{2}T{V_{2}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}+1}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}\\&+ \frac{k_{2}n_{E_{2}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}-1}_{E_{2}}\rho^{n_{I_{2}}+1}_{I_{2}}\\&\quad+ \frac{p_{2}{I_{2}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}+1}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}} \\&+ \frac{\delta_{2}{I_{2}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}-1}_{I_{2}}\\&\quad+ \frac{c_{2}{V_{2}}}{Z}\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}-1}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}. \end{array} $$
(2)
Although the general solution of this expression is intractable, the CTMC assumption of independent events means that the functional equation of \(\xi (\vec {m})\) can be reduced to an algebraic equation. Thus the extinction probability from a given state is the product of the extinction probabilities from each of the constituents of that state [44], so we can write
$$\begin{array}{*{20}l} \xi(\vec{m})&=\xi(n_{E_{1}},n_{I_{1}}, n_{V_{1}}, n_{E_{2}}, n_{I_{2}}, n_{V_{2}})\\&\quad=\rho^{n_{V_{1}}}_{V_{1}}\rho^{n_{E_{1}}}_{E_{1}}\rho^{n_{I_{1}}}_{I_{1}}\rho^{n_{V_{2}}}_{V_{2}}\rho^{n_{E_{2}}}_{E_{2}}\rho^{n_{I_{2}}}_{I_{2}}, \end{array} $$
(3)
where \(\rho ^{n_{V_{1}}}_{V_{1}}\) is the probability that the virus, V1, initiates a process with \(n_{V_{1}}\) number of virus particles that results in extinction. In a similar manner, \(\rho ^{n_{E_{1}}}_{E_{1}}, \rho ^{n_{I_{1}}}_{I_{1}}\) and others are probabilities for eclipse cell, E1, or infected cell, I1 and so on. Eq. (3) is recognizable as the fixed-point equation \(\vec {\varepsilon }= P(\vec {\varepsilon })\), where \(\vec {\varepsilon }= [\varepsilon _{1}, \ldots, \varepsilon _{J}]\) and \(P(\vec {\varepsilon })\) is the probability generating function of the progeny distributions. Now substituting Eq. (3) in Eq. (2), we get
$$\begin{array}{*{20}l} \rho^{n_{V_{1}}}_{V_{1}}&=\frac{\beta_{1}T}{\beta_{1}T+c_{1}}\rho^{n_{V_{1}}}_{V_{1}}\rho_{E_{1}}+\frac{c_{1}}{\beta_{1}T+c_{1}}\rho^{n_{V_{1}}-1}_{V_{1}}\\ \text{or,}\ \rho_{V_{1}}&=\frac{\beta_{1}T}{\beta_{1}T+c_{1}}\rho_{V_{1}}\rho_{E_{1}}+\frac{c_{1}}{\beta_{1}T+c_{1}}\text{.}\\ \rho^{n_{I_{1}}}_{I_{1}}&=\frac{p_{1}}{p_{1}+\delta_{1}}\rho_{V_{1}}\rho^{n_{I_{1}}}_{I_{1}}+\frac{\delta_{1}}{p_{1}+\delta_{1}}\rho^{n_{I_{1}}-1}_{I_{1}}\\ \text{or,}\ \rho_{I_{1}}&=\frac{p_{1}}{p_{1}+\delta_{1}}\rho_{V_{1}}\rho_{I_{1}}+\frac{\delta_{1}}{p_{1}+\delta_{1}},\\ \rho^{n_{E_{1}}}_{E_{1}}&=\rho^{n_{E_{1}}-1}_{E_{1}}\rho_{I_{1}},\ or \ \rho_{E_{1}}=\rho_{I_{1}}\ \text{and}\\ \rho^{n_{V_{2}}}_{V_{2}}&=\frac{\beta_{2}T}{\beta_{2}T+c_{2}}\rho^{n_{V_{2}}}_{V_{2}}\rho_{E_{2}}+\frac{c_{2}}{\beta_{2}T+c_{2}}\rho^{n_{V_{2}}-1}_{V_{2}}\\ \text{or,}\ \rho_{V_{2}}&=\frac{\beta_{2}T}{\beta_{2}T+c_{2}}\rho_{V_{2}}\rho_{E_{2}}+\frac{c_{2}}{\beta_{2}T+c_{2}},\\ \rho^{n_{I_{2}}}_{I_{2}}&=\frac{p_{2}}{p_{2}+\delta_{2}}\rho_{V_{2}}\rho^{n_{I_{2}}}_{I_{2}}+\frac{\delta_{2}}{p_{2}+\delta_{2}}\rho^{n_{I_{2}}-1}_{I_{2}}\\ \text{or,}\ \rho^{n_{I_{2}}}_{I_{2}}&=\frac{p_{2}}{p_{2}+\delta_{2}}\rho_{V_{2}}\rho_{I_{2}}+\frac{\delta_{2}}{p_{2}+\delta_{2}},\\ \rho^{n_{E_{2}}}_{E_{2}}&=\rho^{n_{E_{2}}-1}_{E_{2}}\rho_{I_{2}}\ or \ \rho_{E_{2}}=\rho_{I_{2}}, \end{array} $$
where \(\rho _{V_{i}}\), \(\rho _{I_{i}}\) and \(\rho _{E_{i}}\) are the extinction probabilities when the processes are initiated with a single virus particle or eclipse cell or infectious cell. Solving for each probability, we get \(\rho _{V_{i}}=1\) and \(\rho _{V_{i}}= \frac {c_{i}(p_{i}+\delta _{i})}{p_{i}(c_{i}+\beta _{i}T)}\), \(\rho _{I_{i}}=1\) and \(\rho _{I_{i}}=\frac {\delta _{i}(c_{i}+\beta _{i}T)}{\beta _{i}T(p_{i}+\delta _{i})}\), and \(\rho _{E_{i}}=\rho _{I_{i}}\) where i=1,2. Since probability has to be less than or equal to 1, we can write the solutions of the extinction probabilities as:
$$\begin{array}{*{20}l} \rho_{V_{1}}&=\text{min}\left(\frac{c_{1}(p_{1}+\delta_{1})}{p_{1}(c_{1}+\beta_{1}T)}, 1\right),\\ \rho{I_{1}}&=\text{min}\left(\frac{\delta_{1}(c_{1}+\beta_{1}T)}{\beta_{1}T(p_{1}+\delta_{1})}, 1\right), \\ \rho_{E_{1}}&=\rho{I_{1}},\\ \rho_{V_{2}}&=\text{min}\left(\frac{c_{2}(p_{2}+\delta_{2})}{p_{2}(c_{2}+\beta_{2}T)}, 1\right),\\ \rho{I_{2}}&=\text{min}\left(\frac{\delta_{2}(c_{2}+\beta_{2}T)}{\beta_{2}T(p_{2}+\delta_{2})}, 1\right),\\ \rho_{E_{2}}&=\rho{I_{2}}\text{.} \end{array} $$
Probability of virus extinction Since the extinction of each event is independent, we can write for the probability that both viruses go extinct if the simultaneous infection is initiated with a single virus of each type by the expression \(\rho _{V_{1}}\rho _{V_{2}}\),
$$\rho_{V_{1}}\rho_{V_{2}}=\frac{c_{1}(p_{1}+\delta_{1})}{p_{1}(c_{1}+\beta_{1}T_{0})}\frac{c_{2}(p_{2}+\delta_{2})}{p_{2}(c_{2}+\beta_{2}T_{0})}\text{.} $$
Stochastic dynamics of identical viruses
While the probability of virus extinction is an important feature of stochastic models, we are also interested in understanding if stochasticity affects the predicted dynamics of coinfections that survive. Previously in our ODE model [27], we found that the virus with the higher growth rate always out-competes the slower growing virus. While ODEs can give us the average behaviors of the coinfection process, in real systems the biological processes are stochastic. The randomness associated with births and deaths during the initial infection process may lead to virus extinction even in an exponentially growing virus population [45]. Yan et al. [34] reported that the invasion of viral infection is dependent on the initial viral dose and growth rate of each virus. Here, we are interested in knowing how the coinfection dynamics change with change in growth rates of each virus. First, we will observe the dynamics of coinfection with identical viruses.
Keeping all the initial conditions and transition rates for both viruses equal, we examine the time course of coinfection by plotting number of virus over time. 1000 sample stochastic trajectories of viral load curve for coinfection with identical viruses are shown in Fig. 1. We find that both viruses have peaks above the threshold of detection (100 virions) 88% of the time and 12% of the time one of the viruses experiences extinction. Among 120 (12%) extinctions, virus 1 and virus 2 experience extinction 49 and 65 times out of 1000 simulations, respectively. In another words, there is a 4.9% chance that what begins as a coinfection will result in a single virus infection with virus 2 or 6.5% chance with virus 1.
The ODE model predicts that when all parameters are equal, both viruses will have the same time course, splitting available target cells equally. In the stochastic model, we find that despite having identical growth rates, one virus out-competes the other virus in particular realizations of the model. Virus 1 has a higher peak viral titer 513 times within 1000 simulations, while virus 2 has the higher peak viral titer 487 times. So while a particular realization of the model will have a clear dominant virus, on average, the viruses are equivalent, in agreement with the ODE model. Additional file 1 includes additional figures examining the distributions when viruses differ. To characterize the viral time course, we calculate peak viral load, time of peak for each virus, as well as duration of coinfection (Fig. 2). The median time of peak for virus 1 is found to be 2.374 ±0.64 days and for virus 2, it is 2.375 ±0.65 days. The median of peak viral load for virus 1 and 2 are (4.0 ±2.6) ×107 and (4.1 ±2.6) ×107, respectively. From the distributions (Fig. 2), we see that even if the viruses behave differently for a particular model realization, on average, they tend to behave identically. Finally, the distribution of coinfection duration is given in Fig. 2 where the median coinfection duration is found to be 5.730 ±0.059 days. Despite fluctuations in the time course of each virus, the coinfection duration does not vary much.
Stochastic dynamics for different viruses
Since growth rate determines which virus is the stronger competitor [27], we investigate how differences in growth rate between the two viruses changes stochastic infections. Unfortunately, growth rate is not a parameter in the model, so we need to determine which model parameter(s) to change to systematically vary the growth rate. We use the expression for growth rate derived by Smith et al. [46] and determine how growth rate depends on different model parameters (Fig. 3). We find that growth rate varies approximately linearly with virus production rate, p, over a large range of p (p>1), so we will systematically alter p for one virus to alter its growth rate.
For ease of interpretation, we define the relative viral production rate \(r = \frac {p_{1}}{p_{2}}\). We first examine how competition between the viruses changes as the relative growth rates change. Here variation is introduced for virus 1 keeping virus 2 fixed for a range, r=1×10−1×102. We count the number of times, out of 1000 simulated infections, a particular virus has a higher viral titer peak than the other virus. The results are shown in Fig. 4. When the viruses have identical growth rates, there is a 50% chance that a particular virus will have the higher peak titer, as seen in the previous section. The probability of having a higher peak viral load increases rapidly as the production rate of a virus increases, reaching 90% with a less than 2-fold change in viral production. Note that the probability of having the higher peak viral titer never quite reaches 100%, even when there are large differences in growth rate. This indicates that early stochastic events can significantly alter the time course of the infection.
In Fig. 5, we compare coinfection dynamics for the ODE and CTMC models, looking at peak viral load, time of viral peak, and coinfection duration. ODEs predict that if the growth rate of one virus is higher than the other it will always have a higher peak viral load (Fig. 5 (top left)). For the CTMC model, the transition from one virus dominating to the other dominating is not as sharp. Unlike the predictions of ODEs, the CTMC allows for the slower growing virus to dominate infection dynamics. In fact, the median peak viral loads for virus 1 and virus 2 cross closer to a relative viral production rate of 101 rather than 100 as seen in the ODE model. Stochastic variability in the peak viral load (as indicated by the shaded area) for both viruses overlaps for a wide range of relative viral production, indicating that the viruses can have similar peak viral loads.
The time of viral peak also shows some differences between the ODE and CTMC models. For the ODE model, the time of viral peak is similar for both viruses when the relative viral production rate is greater than 100, although the time of peak decreases as the relative viral production rate increases. This is because the viral production rate of virus 1 is increased over its baseline value causing an earlier time of peak. This drives the earlier time of peak of virus 2, which is the weaker competitor in this case. The decline in time of viral peak is not as sharp in the CTMC model since stochasticity can temper the effect of the increased production rate of virus 2 by allowing virus 1 to still have an opportunity to infect some cells.
Finally, we compare the predicted duration of coinfection variation for ODE and stochastic models (Fig. 5 (bottom row)). Viruses do not coexist for more than about a week in either model. The longest coinfection durations are seen, for both models, when the two viruses have the same growth rates. This is because the faster growing virus out-competes the slower growing virus leading to short infections for the slower growing virus.
One feature of viral infections that cannot be captured by ODE models is extinction of the infection. Therefore, we simulate the probability of virus extinction, defined as the fraction of times when one virus does not grow above the virus detection limit (detection limit is equal to 100 virus particles), when the coinfection is initiated with a single virus of each type (Fig. 6). Note that this is slightly different than the definition for the probability of extinction calculated in “Derivation of extinction coefficient” section which requires that virus, along with infectious and eclipse cells, all go to zero. The probability of both viruses growing to detectable levels is highest for viruses with similar relative production rates. When relative viral production rates are very different (about 10–100 fold difference), there is a high probability that one virus becomes extinct. When the viruses have very different production rates, the virus with a higher production rate will out-compete the one with a low production rate driving it to extinction. However, since one virus (in this case virus 1) experiences decreased production rate from the base value but initiates infection with the same amount of virus, the probability of extinction reaches close to 100% in a quicker manner for lower relative production rate than that of the higher relative rates.