1 Introduction

A host viral infection can be established with just one introduced virion. However, there is a chance that deaths of virions and infected cells occur such that the infection fails and the virus population becomes extinct. As a result, not every exposure to virus leads to systemic infection; there is a chance element in the number of times an individual is exposed to a pathogen before systemic infection is established, as demonstrated by animal studies in diseases such as SIV (Keele et al. 2009). Factors which lead to viral extinction before systemic infection is established include the innate and adaptive immune responses, spatial heterogeneity in infection processes, and the stochastic nature of processes such as infection of cells by virions, production of virions from cells and death of virions and infected cells.

The extinction probability when a number of virions is introduced to the host is of interest, not only because it reflects important information about the reproduction characteristics of the virus and the effectiveness of the mode of transmission, but also because there may be ways to manipulate this extinction probability and defend the host against infection. Many studies of the host immune response and of pre/post-exposure drug prophylaxis have concentrated on their effect on the viral load (Beauchemin and Handel 2011; Dobrovolny et al. 2013; Smith and Ribeiro 2010) rather than the extinction probability, but the ability to induce stochastic extinction is another way to quantify the effectiveness of the immune response or intervention.

Studies of successive respiratory viral infections in animals have shown that morbidity and mortality are reduced when infections are separated by days or weeks (Seo and Webster 2001; Walzl et al. 2000; Bodewes et al. 2011; Laurie et al. 2010); in the case of different influenza A subtypes, the subsequent infection can be prevented altogether (Laurie et al. 2015). The reduction of viral shedding can be attributed to cross-reactive adaptive immunity and a heightened innate immune response. We hypothesise that subsequent infections are prevented altogether when the immune response is so effective that stochastic extinction of the virus becomes likely (Cao et al. 2015). The proportion of prevented subsequent infections depends on the interval between primary and subsequent exposures, so the effect is time-dependent. Studies have also shown that pre-exposure interferon treatment in animals can prevent infection with pandemic influenza (Steel et al. 2010), and that both pre- and post-exposure prophylaxis reduce the probability of HIV infection (Cardo et al. 1997; Baeten et al. 2012; Grant et al. 2010; Thigpen et al. 2012). Because infection is often initiated by only one or a few virions (Keele et al. 2008), stochasticity plays an important role in the effectiveness of prophylaxis. The short timeframe within which post-exposure prophylaxis must be administered suggests that the extinction probability is time-dependent.

In the absence of a time-dependent immune response, stochastic models for HIV incorporating target cells and virions have been studied using Monte Carlo simulations of the multi-type branching process (Heffernan and Wahl 2005; Kamina et al. 2001; Tan and Wu 1998), or by simulating solutions to stochastic differential equations where the infection and death processes are diffusion processes, represented by noise terms in the equations (Tuckwell and Corfec 1998). Kamina et al. (2001) has shown that branching process simulations estimate the extinction probability more accurately than diffusion approximations. However, estimation of the extinction probability using simulations can be computationally costly; this makes the procedure infeasible for certain situations. For example, when conducting sensitivity analyses, sampling procedures generate many different parameter sets for which the extinction probability must be simulated. In addition, when fitting models to experimental data using methods such as maximum likelihood estimation, the likelihood of the observed data given model parameters at each point in parameter space depends on the extinction probability at that point. Analytic results for the extinction probability have been obtained, but most of these have either only kept track of the number of infected cells (Merrill 2005) or the number of virions (Tuckwell et al. 2008).

Pearson et al. (2011) derived analytic expressions for the extinction probability when an arbitrary number of infected cells and/or virions is introduced to the host. This method provides a fast, accurate calculation of the extinction probability, but it is assumed that the rates of target cell infection by a single virion, death of a single virion, production of a virion by a single cell, and death of a single cell are constant throughout the early stage of infection, when the extinction probability is non-negligible. This is not necessarily the case if the immune response is changing during the early stage of infection, such as through vaccination, prophylaxis or previous infection. Conway et al. (2013) extended the work of Pearson et al. (2011) to cover the case when there is explicit time-dependence of the infection rate and viral production rate, and used the results to analyse the effect of pre- and post-exposure prophylaxis for HIV. In general, the time dependence of the model parameters may not be explicit; the parameters may be functions of other model compartments, which may be time-varying.

Both Pearson et al. (2011) and Conway et al. (2013) have assumed that infected cells are immediately able to produce virus, whereas biologically there is a delay, often termed the latent period, caused by intracellular events which must happen before the production of virions. In addition, both Pearson et al. (2011) and Conway et al. (2013) have assumed that the mean lifetime of an infected cell is exponentially-distributed, such that the most likely mean lifetime is zero. Since this is biologically unlikely, models with different latent and infectious periods have been developed. Adding these complexities into the model lead to more realistic viral dynamics (Jensen 1948; Grossman et al. 1998; Lloyd 2001), and makes model parameters more biologically plausible (Baccam et al. 2006; Beauchemin et al. 2008).

We apply the method presented by Pearson et al. (2011) to calculate the extinction probability at different stages of infection for a within-host viral model with a time-dependent immune response modelled by ordinary differential equations, with multiple latent and infectious stages. We take a different approach to Conway et al. (2013) by directly solving for the probability function for extinction from a given state and time, rather than using the probability generating function. Using a case where the viral infection is strongly affected by the time dependence of the immune response (Laurie et al. 2015; Cao et al. 2015),we show that the calculations agree with the results obtained using Monte Carlo simulations, and examine the change in extinction probability when the number of latent and infectious stages, and the mean duration of these stages, changes.

2 Derivation of the extinction probability

2.1 The TIV model with continuous production of virions

A common model for within-host viral infection is the TIV model (Perelson 2002). The transitions are

$$\begin{aligned} V+T&\overset{\beta }{\rightarrow } I \end{aligned}$$
(1a)
$$\begin{aligned} I&\overset{\delta }{\rightarrow } \varnothing \end{aligned}$$
(1b)
$$\begin{aligned} I&\overset{p}{\rightarrow } I + V \end{aligned}$$
(1c)
$$\begin{aligned} V&\overset{c}{\rightarrow } \varnothing \end{aligned}$$
(1d)

where T is the number of uninfected target cells, I is the number of infected cells, V is the number of infectious virions, and \({\rightarrow }\varnothing \) denotes transition to the empty set i.e. clearance. Homogeneous mixing between uninfected cells and virions is assumed. Uninfected cells become infected by virus at rate \(\beta \) per target cell per virion. Infected cells produce infectious virions at a rate p per cell and die at a rate \(\delta \) per cell, where \(1/\delta \) is the mean life span of a productively infected cell; they also produce noninfectious virions, but as these do not contribute to the infection dynamics, they are ignored. Free virus is cleared at rate c per virion, as well as being lost to entry into target cells. This model is considered ‘continuous’ because the infected cells continually produce virions during their lifetime, in contrast to the ‘burst’ model where all virions are released when a cell dies. In this study, we will focus on the continuous production stochastic model, but the same methods can be applied to the burst model. The mean-field kinetics of both models are given by the deterministic equations

$$\begin{aligned} \frac{dT}{dt}&= -\beta TV \end{aligned}$$
(2a)
$$\begin{aligned} \frac{dI}{dt}&= \beta TV - \delta I \end{aligned}$$
(2b)
$$\begin{aligned} \frac{dV}{dt}&= pI - cV - \beta TV . \end{aligned}$$
(2c)

Some versions of the above deterministic model have a scaling factor in front of the \(\beta TV\) term in Eqs. 2a and 2b; this is because the number of infectious virions in the host usually cannot be measured in absolute terms, but only in terms of units such as plaque forming units (PFU), 50 % egg infectious doses (EID\(_{50}\)) or 50 % tissue culture infectious doses (TCID\(_{50}\)). Often, even these cannot be measured directly, and one only knows the concentration of these units in samples drawn from the host. The scaling factor is the ratio between the number of virions required to infect a cell and the number of virions corresponding to one measurement unit; if we assume that V is the absolute number of virions, and each cell is infected by exactly one virion, the scaling factor is not required. We will consider the most commonly modelled case, where infection of a target cell requires only one virion. A detailed study of the feasibility of different binding processes is beyond the scope of this paper and is left for future work.

The basic reproduction number \(R_0\) is given by

$$\begin{aligned} R_0 = \frac{\beta p T_0}{\delta (c + \beta T_0)} = \gamma \alpha \end{aligned}$$
(3)

where \(\gamma = \beta T_0/(c+\beta T_0)\) is the probability that a virion infects a cell and \(\alpha = p/\delta \) is the mean number of virions produced by an infected cell in its lifetime. It is defined as the number of secondary infected cells due to a single infected cell in a population of uninfected cells of size \(T_0\), i.e. the number of cells infected by virions produced directly by the first infected cell.

2.2 The extinction probability for the TIV model

For a time-independent immune response (constant \(\beta \),\(\delta \),p,c), the system is most susceptible to stochastic extinction in early infection, when virion and infected cell numbers are small. At this stage, the change in target cell numbers due to infection is small, so T is approximately constant (\(T \approx T_0\)) and we can simplify Eq. 1 to

$$\begin{aligned} V&\overset{\beta T}{\rightarrow } I \end{aligned}$$
(4a)
$$\begin{aligned} I&\overset{\delta }{\rightarrow } \varnothing \end{aligned}$$
(4b)
$$\begin{aligned} I&\overset{p}{\rightarrow } I + V \end{aligned}$$
(4c)
$$\begin{aligned} V&\overset{c}{\rightarrow } \varnothing . \end{aligned}$$
(4d)

The continuous-time Markov chain becomes a multi-type branching process. The extinction probabilities when the system is initiated with a single virion or a single infected cell are given respectively by (Pearson et al. 2011)

$$\begin{aligned} \rho _V&= \mathrm {min}\left( 1-\frac{R_0 - 1}{\alpha },1\right) \end{aligned}$$
(5a)
$$\begin{aligned} \rho _I&= \mathrm {min}(1/R_0,1). \end{aligned}$$
(5b)

If the system is initiated with \(n_V\) virions and \(n_I\) infected cells, because the extinctions of the lineages of each virus and infected cell are independent, the extinction probability is given by

$$\begin{aligned} \varepsilon (n_V,n_I) = \rho _V^{n_V}\rho _I^{n_I}. \end{aligned}$$
(6)

In natural infection or inoculation experiments where the hosts are initially naive to infection, the initial condition is \(n_I = 0, n_V > 0\).

It is worth noting that for acute diseases such as influenza where the viral load peak is orders of magnitude higher than the initial value, many studies neglect the loss of virions due to entry into target cells to initiate infection (Lee et al. 2009; Miao et al. 2010; Baccam et al. 2006; Saenz et al. 2010; Pawelek et al. 2012). By neglecting virion loss due to cell entry, the transitions for the model become

$$\begin{aligned} V&\overset{\beta T}{\rightarrow } I + V \end{aligned}$$
(7a)
$$\begin{aligned} I&\overset{\delta }{\rightarrow } \varnothing \end{aligned}$$
(7b)
$$\begin{aligned} I&\overset{p}{\rightarrow } I + V \end{aligned}$$
(7c)
$$\begin{aligned} V&\overset{c}{\rightarrow } \varnothing . \end{aligned}$$
(7d)

This simplification assumes that \(\beta T \ll c\), which becomes more valid as the infection progresses and target cell depletion occurs. As the aforementioned studies are not focussed on the very early stages of infection, the assumption is used to simplify the form of the model. However, for very early infection where the target cells have not been depleted, the assumption is not necessarily valid. The extinction probabilities become

$$\begin{aligned} \rho _V&= \mathrm {min}\left( \frac{c(\delta +p)}{p(\beta T + c)},1\right) \end{aligned}$$
(8a)
$$\begin{aligned} \rho _I&= \mathrm {min}\left( \frac{\delta (\beta T+c)}{\beta T(\delta +p)},1\right) . \end{aligned}$$
(8b)

By rewriting Eq. 5 as

$$\begin{aligned} \rho _V&= \mathrm {min}\left( \frac{c(\delta +p) + \delta \beta T}{p(\beta T + c)},1\right) \end{aligned}$$
(9a)
$$\begin{aligned} \rho _I&= \mathrm {min}\left( \frac{\delta (\beta T+c)}{\beta Tp},1\right) , \end{aligned}$$
(9b)

one can see that neglecting loss of virus due to entry into cells underestimates the extinction probability, as virions have the opportunity to infect more than one cell.

2.3 Extension to more realistic models with time-independent parameters

The continuous production TIV model has been extended to include an exponentially-distributed latent stage where infected cells do not yet produce virions, reflecting the delay in virion production due to biological processes which must happen in the cell to induce virion production. For many diseases, such as influenza and HIV, the length and form of the latent phase affect parameter estimates (Herz et al. 1996; Mittler et al. 1998; Nelson et al. 2000), and the parameters in the model take on more realistic values when fit to data if a latent stage is included (Baccam et al. 2006; Beauchemin et al. 2008). The extinction probability in this case is a function of the initial number of virions, latently infected cells and actively infected cells, and the absorbing state is when all of these quantities are zero.

In the general case there are multiple latent and infectious stages, each with exponentially-distributed transition times, and the simplest model recovered as a special case with one latent and one infectious stage. These models use the method of stages to make the distribution of the total time spent in the latent and infectious period more realistic (Jensen 1948; Grossman et al. 1998; Lloyd 2001), such that the most likely transition time is no longer zero, as in the case of the exponential distribution. The resulting distributions for the total latent and infectious period distributions are Erlang distributions with means \(1/\delta _L\) and \(1/\delta _I\), and modes \((K-1)/(K\delta _L)\) and \((M-1)/(M\delta _I)\) respectively, where K and M are the number of latent and infectious stages respectively. This is analogous to having gamma-distributed latent and infectious periods in an epidemic model (Anderson and Watson 1980). Because the transition times for each stage are exponentially-distributed, the resulting model is still a Markov chain, and the procedure by Pearson et al. (2011) can be used to calculate the extinction probability. The extinction probability is then a function of the initial number of virions, latently infected cells in each stage, and actively infected cells in each stage, and the absorbing state is when these \(1+K+M\) quantities are all zero.

The transitions for the continuous production model with K latent stages and M infectious stages are given by

$$\begin{aligned}&V \overset{\beta T}{\rightarrow } L_1 \end{aligned}$$
(10a)
$$\begin{aligned}&L_k \overset{K\delta _L}{\rightarrow } L_{k+1}, \quad k = 1,2,\ldots ,K-1 \end{aligned}$$
(10b)
$$\begin{aligned}&L_K \overset{K\delta _L}{\rightarrow } I_1 \end{aligned}$$
(10c)
$$\begin{aligned}&I_m \overset{M\delta _I}{\rightarrow } I_{m+1}, \quad m = 1,2,\ldots ,M-1 \end{aligned}$$
(10d)
$$\begin{aligned}&I_M \overset{M\delta _I}{\rightarrow } \varnothing \end{aligned}$$
(10e)
$$\begin{aligned}&I_m \overset{p}{\rightarrow } I_m + V, \quad m = 1,2,\ldots ,M \end{aligned}$$
(10f)
$$\begin{aligned}&V \overset{c}{\rightarrow } \varnothing \end{aligned}$$
(10g)

for \(K > 0\) and

$$\begin{aligned}&V \overset{\beta T}{\rightarrow } I_1 \end{aligned}$$
(11a)
$$\begin{aligned}&I_m \overset{M\delta _I}{\rightarrow } I_{m+1}, \quad m = 1,2,\ldots ,M-1 \end{aligned}$$
(11b)
$$\begin{aligned}&I_M \overset{M\delta _I}{\rightarrow } \varnothing \end{aligned}$$
(11c)
$$\begin{aligned}&I_m \overset{p}{\rightarrow } I_m + V, \quad m = 1,2,\ldots ,M \end{aligned}$$
(11d)
$$\begin{aligned}&V \overset{c}{\rightarrow } \varnothing \end{aligned}$$
(11e)

for \(K = 0\).

The independence of virion/cell lineage extinction can be written in this case as

$$\begin{aligned} \varepsilon (n_V,n_{L1},\ldots ,n_{LK},n_{I1},\ldots ,n_{IM}) = \rho _V^{n_V}\prod _{k = 1}^K \rho _{Lk}^{n_{Lk}}\prod _{m = 1}^M \rho _{Im}^{n_{Im}}. \end{aligned}$$
(12)

Following the procedure of Pearson et al. (2011), we obtain the following set of simultaneous equations for the extinction probabilities:

$$\begin{aligned}&\rho _V = \frac{1}{\beta T+c}(\beta T \rho _{I1} + c) \end{aligned}$$
(13a)
$$\begin{aligned}&\rho _{L1} = \rho _{L2} =\ldots = \rho _{LK} = \rho _{I1} \end{aligned}$$
(13b)
$$\begin{aligned}&\left[ \begin{array}{llll} \chi &{} \eta &{} &{} \\ &{} \ddots &{} \ddots &{} \\ &{} &{} \ddots &{} \eta \\ &{} &{} &{} \chi \end{array}\right] \begin{bmatrix} \rho _{I1}\\ \rho _{I2}\\ \vdots \\ \rho _{IM} \end{bmatrix} = \begin{bmatrix} 0\\ \vdots \\ 0 \\ -\eta \end{bmatrix} \nonumber \\&\text {where} \qquad \eta = \frac{M\delta _I}{p+M\delta _I} \nonumber \\&\chi = \frac{p \rho _V}{p+M\delta _I} - 1. \end{aligned}$$
(13c)

The M-by-M matrix in Eq. 13c is a Toeplitz tridiagonal matrix (with zeros in the diagonal below the main diagonal); inverting it yields (Huang and McColl 1997)

$$\begin{aligned} \begin{bmatrix} \rho _{I1}\\ \rho _{I2}\\ \vdots \\ \rho _{IM} \end{bmatrix} = \begin{bmatrix} \mathrm {min}((-\eta /\chi )^M,1)\\ \mathrm {min}((-\eta /\chi )^{M-1},1)\\ \vdots \\ \mathrm {min}((-\eta /\chi ),1) \end{bmatrix}. \end{aligned}$$
(14)

Note that \(\eta /\chi < 0\). Equation 13a can be solved simultaneously with the first row of Eq. 14 to obtain a numerical result for \(\rho _V\) and \(\rho _{I1}\), which can then be substituted into Eqs. 13b and 14 to find the remaining extinction probabilities.

The extinction probability for a latent cell in any stage is equal to the extinction probability for an infectious cell in the first stage, because for each latent stage, there is only one possible transition. Thus, the latent cell transitions to an infectious cell in the first stage, with an unchanged number of virions, with probability 1, and the latency does not affect the extinction probability. The extinction probability of an infected cell in stage \(M-m\) is equal to the extinction probability in stage M, to the power of m. This is because the infected cell is guaranteed to pass through each stage, and the stages are identical. Consequently, the probability mass function for the number of virions produced by an infected cell in m stages is the same as the probability mass function for the number of virions produced by m infected cells in one of those stages, and the extinction probability is the same in both cases. The independence of extinction of each infected cell leads to the aforementioned result.

2.4 Extension to a time-dependent immune response

For a time-dependent immune response, we focus on the case where during early infection, when virion and infected cell counts are low and stochastic extinction is most likely, the contribution to the immune response by the low number of virions and infected cells is small; rather, the immune response is mainly induced by factors such as a pre-existing infection or prophylaxis. The pre-existing infection or prophylaxis may also change the number of target cells over time, such that the assumption in the previous subsections that T is constant is no longer valid. The unidirectional interaction between immune response and virus/infected cells allows a generalisation of the TIV model parameters to be explicit functions of time.

We use \(\varepsilon (\overrightarrow{m},t)\) to denote the extinction probability when the system starts from the state \(\overrightarrow{m}\) at time t. We need the probabilities \(p_i(\overrightarrow{m},t,t+\tau )\) that at time t, given the state \(\overrightarrow{m}\), the ith reaction (infection of a target cell, death of an infected cell, production of a virion, or death of a virion) is the next reaction, and it occurs at time \(t+\tau \). We must integrate over all possible times to next event \(\tau \) to find the extinction probability at time t.

This is summarised by the equation

$$\begin{aligned} \varepsilon (\overrightarrow{m},t) = \int _0^\infty \sum _{i=1}^{R}p_i(\overrightarrow{m},t,t+\tau )\varepsilon (\overrightarrow{m}+d\overrightarrow{m}_i,t+\tau ) d\tau , \quad \overrightarrow{m} \ne \overrightarrow{0}, \end{aligned}$$
(15)

where \(R = 4\) is the total number of possible reactions, and the boundary condition

$$\begin{aligned} \varepsilon (\overrightarrow{0},t) = 1. \end{aligned}$$
(16)

This is a time-dependent generalisation of Eqs. 8 and 9 in Pearson et al. (2011). Equations 6 and 12 still hold, but the probabilities are now time-dependent:

$$\begin{aligned} \varepsilon (n_V,n_{L1},\ldots ,n_{LK},n_{I1},\ldots ,n_{IM},t) = \rho _V(t)^{n_V}\prod _{k = 1}^K \rho _{Lk}(t)^{n_{Lk}}\prod _{m = 1}^M \rho _{Im}(t)^{n_{Im}}. \end{aligned}$$
(17)

The time-dependent event probabilities \(p_i(\overrightarrow{m},t,t+\tau )\) are given by (Lu et al. 2004)

$$\begin{aligned} p_i(\overrightarrow{m},t,t+\tau ) = r_i(t+\tau ) \exp \left[ -\sum _{s=1}^{R} \int _t^{t+\tau } r_s(\tau ') d \tau ' \right] . \end{aligned}$$
(18)

Combining Eqs. 18 and 15, we obtain (after some algebra, including differentiating on both sides to turn the integral equation into a differential equation)

$$\begin{aligned} \frac{d}{dt}\left[ \varepsilon (\overrightarrow{m},t)\right] =\sum _{i=1}^{R} r_i(t)\left[ \varepsilon (\overrightarrow{m},t) - \varepsilon (\overrightarrow{m}+d\overrightarrow{m}_i,t)\right] . \end{aligned}$$
(19)

For the TIV model, substituting \(\overrightarrow{m} = (n_V,n_I) = (1,0)\) and \(\overrightarrow{m} = (0,1)\) separately into Eq. 19 yields the two simultaneous equations

$$\begin{aligned} \frac{d\rho _V}{dt}&= \beta (t)T(t)\left[ \rho _V - \rho _I\right] + c(t)\left[ \rho _V - 1\right] \end{aligned}$$
(20a)
$$\begin{aligned} \frac{d\rho _I}{dt}&= \delta (t)\left[ \rho _I - 1\right] + p(t)\left[ \rho _I - \rho _V\rho _I\right] . \end{aligned}$$
(20b)

Similarly, for the model with multiple latent and/or infectious stages, substituting \(\overrightarrow{m} = (n_V,n_{L1},\ldots ,n_{LK},n_{I1},\ldots ,n_{IM}) = (1,0,\ldots ,0), (0,1,0,\ldots ,0),\ldots ,(0,\ldots ,0,1)\) in turn, we obtain the following set of simultaneous equations for the extinction probabilities:

$$\begin{aligned} \frac{d\rho _V}{dt}&= \beta (t)T(t)\left[ \rho _V - \rho _{L1}\right] + c(t)\left[ \rho _V - 1\right] \end{aligned}$$
(21a)
$$\begin{aligned} \frac{d\rho _{Lk}}{dt}&= K\delta _L(t)[\rho _{Lk} - \rho _{Lk+1}], \quad k = 1,\ldots ,K-1\end{aligned}$$
(21b)
$$\begin{aligned} \frac{d\rho _{LK}}{dt}&= K\delta _L(t)[\rho _{LK} - \rho _{I1}] \end{aligned}$$
(21c)
$$\begin{aligned} \frac{d\rho _{Im}}{dt}&= M\delta _I(t)[\rho _{Im} - \rho _{Im+1}] + p(t)[\rho _{Im} - \rho _V\rho _{Im}], \quad m = 1,\ldots ,M-1 \end{aligned}$$
(21d)
$$\begin{aligned} \frac{d\rho _{IM}}{dt}&= M\delta _I(t)\left[ \rho _{IM} - 1\right] + p(t)\left[ \rho _{IM} - \rho _V\rho _{IM}\right] . \end{aligned}$$
(21e)

Equations 20 and 21 are a set of coupled ordinary differential equations which can be solved given initial conditions; in other words, if we know the extinction probabilities for introduction of one virion, introduction of one latent cell at each stage and introduction of one infected cell at each stage at any particular time, we can solve for the extinction probabilities at any time. Often, the time-dependent immune response is such that the parameters approach a constant value as \( t \rightarrow \infty \), so \(\rho _V(t)\), \(\rho _{Lk}(t)\) and \(\rho _{Im}(t)\) approach their time-independent values as given in Eq. 13. We can thus set the values of \(\rho _V(t)\), \(\rho _{Lk}(t)\) and \(\rho _{Im}(t)\) at some time t much longer than the timescale of the infection, and solve backwards in time to find the values of \(\rho _V(t)\), \(\rho _{Lk}(t)\) and \(\rho _{Im}(t)\) when the immune response is changing most. For the time-dependent case, it is no longer true that \(\rho _{L1}(t) = \cdots = \rho _{LK}(t) = \rho _{I1}(t)\), or that the extinction probability is independent of the mean and distribution of the latent period, because changing the time spent in the latent period shifts the period for which the cell is infectious, and the immune response during that period will be different.

3 Dependence of the extinction probability of virions on the number and mean duration of latent and infectious stages

Having developed a method for calculating the extinction probability, we examine how the extinction probability of introduced virions depends on the form of the model chosen, i.e. the number of latent and infectious stages and the mean latent and infectious periods, for both the TIV model and a model with a time-dependent immune response. We focus on the extinction probability when virions rather than infected cells or latent cells are introduced, because this is the route of transmission for both natural infection and inoculation experiments.

3.1 TIV model

From Eq. 13, we can see that the extinction probability of virions is independent of the number and duration of latent stages, but is dependent on the number of infectious stages and the mean infectious period. We vary the number of infectious stages and the mean infectious period, and fix the remaining parameters, to show the effect of these parameters. The fixed parameter values are taken from Cao et al. (2015), which are chosen to be typical values for influenza infection (Table 1).

Table 1 Parameter values for the TIV model with multiple infectious stages, adapted from Cao et al. (2015)
Full size table
Fig. 1
figure 1

Extinction probability when introducing ten virions, varying the mean infectious period \(1/\delta _I\) and the number of infectious stages M, indicated in the legend. Parameters are given in Table 1

Full size image

Figure 1 shows the extinction probability when introducing ten virions, varying the mean infectious period \(1/\delta _I\) and the number of infectious stages M, indicated in the legend. The number of virions is chosen to be 10 to mimic the number of infectious virions introduced to the upper respiratory tract in an influenza infection; justification is given in Online Resource 1. The results for different numbers of virions can be obtained using Eq. 12. The simultaneous equations 13a and the first row of Eq. 14 are solved for \(\rho _{I1}\) numerically using Matlab 2014b’s fzero function, with a starting guess of \(\rho _{I1} = 0\); the result is substituted into Eq. 13a to solve for \(\rho _V\). The code for generating all the figures in this paper is provided at https://bitbucket.org/prism2/extinction_probability_ayan.

We see that when the number of infectious stages is held constant, the extinction probability decreases as the mean infectious period increases. This is because as the expected number of virions produced by an infected cell (\(p/\delta _I\)) increases, the basic reproduction number \(R_0\) increases. As the mean infectious period increases to infinity, the extinction probability of one virion approaches \(c/(\beta T+c)\), because in this limit once the virion infects a cell, the cell lives forever and so the conditional extinction probability is zero; thus, the extinction probability is equal to the probability that the virion dies before infecting a cell. On the other hand, as the mean infectious period decreases such that \(R_0\) approaches 1, the extinction probability approches 1. When the mean infectious period is held constant, the extinction probability decreases as the number of stages increases. If we can determine the parameters p, c, \(\beta \), \(T_0\), \(V_0\) and one of \(\delta _I\) and M accurately using data such as viral load measurements, then data on the extinction probability can help us determine the remaining parameter. The distribution of the latent period, on the other hand, cannot be determined experimentally using the extinction probability, but nevertheless affects other aspects of infection such as the viral load time course.

3.2 A model of re-infection with a time-dependent immune response

We will now apply the method of calculating the extinction probability with a time-dependent immune response to the scenario detailed by Cao et al. (2015), in which the host is sequentially infected with two different influenza viruses. In the experiment which motivated this model, it was observed that when ferrets were infected sequentially with different strains of influenza virus, infection with the second virus was prevented in some ferrets and not others, depending on the strains used and the interval between primary and secondary exposures (the inter-exposure interval) (Laurie et al. 2015). Moreover, for some combinations of strains and inter-exposure intervals, infection with the second virus was only seen in a subset of ferrets, suggesting that stochastic effects are important for this temporary immunity. For those ferrets where infection with the second virus was eventually established, the level of the second virus stayed low for up to days after infection; this supports the hypothesis that stochastic extinction, which is most likely for low virion numbers, was responsible for the prevented infections. Motivated by the results of this experiment, Cao et al. (2015) developed a mechanistic model for strain-dependent and strain-independent components of the immune response to influenza in ferrets, and showed through numerical simulation that stochastic extinction was likely for certain numbers of initial virions, and a chosen inter-exposure interval and set of parameter values. Here, we explore extinction in this model systematically, and evaluate how changing the inter-exposure interval, initial number of virions, mean latent and infectious periods, and number of latent and infectious stages changes the extinction probability. This will lead to a better understanding of the phenomenon of temporary immunity.

The early kinetics of the first virus strain, when extinction is most likely to occur, are not affected by the second strain (which is not yet present), so the infection parameters can be treated as time-independent and the extinction probability is given by Eqs. 6 and 5. However, given the large inoculum and naive state of the host, extinction is unlikely to occur, and indeed was never observed in the experiments. Our interest lies in the early kinetics of the second strain, where the time-dependent immune response initiated by the first strain affects the extinction probability. Figure 2 shows the ways by which the first virus induces a time-dependent immune response [reproduced from Cao et al. (2015)].

Fig. 2
figure 2

Induction of a time-dependent immune response by the first virus. Reproduced from Cao et al. (2015)

Full size image

The target cell pool, T, is shared by the two viral strains. The time-dependent innate immune response is mediated by type I interferon F, which is produced in response to infected cells. It acts in three ways to modify infection kinetics, through:

  1. 1.

    Introduction of a temporary resistant state R in target cells;

  2. 2.

    Inhibition of the production of virions from infected cells; and

  3. 3.

    Direct killing of infected cells.

The model also captures the role of antibodies (A), responsible for strain-specific viral clearance and induction of long-term sterilising immunity, mediated by stimulation of B cells (B). In addition, target cells are replenished at a rate proportional to the product of the number of target and dead cells. Cao et al. (2015) has explored the dynamics of all three mechanisms; to aid exposition, we focus on one mechanism, which we choose to be the second mechanism (inhibition of the production of virions from infected cells). However, the extinction probability can be calculated for each of the mechanisms using the same method.

The equations for the dynamics of virus 1, the target cells and the induced immune response (for the second mechanism) are

$$\begin{aligned} \frac{dV_1}{dt}&= \frac{p_{V1}}{1+s_1F}I_1 - \delta _{V1} V_1 - \kappa _{A1}V_1A_1 - \beta _1 V_1T \end{aligned}$$
(22a)
$$\begin{aligned} \frac{dT}{dt}&= g_{T}T\left( 1 - \frac{T+I_1}{T_{0}}\right) - \beta _1V_1T \end{aligned}$$
(22b)
$$\begin{aligned} \frac{dI_1}{dt}&= \beta _1 V_1T - \delta _{I1} I_1 \end{aligned}$$
(22c)
$$\begin{aligned} \frac{dF}{dt}&= p_{F1}{I_1} - \delta _{F} F \end{aligned}$$
(22d)
$$\begin{aligned} \frac{dB_{1}}{dt}&= m_{11}V_1(1 - B_1) - m_{21}B_1 \end{aligned}$$
(22e)
$$\begin{aligned} \frac{dA_1}{dt}&= m_{31}{B_1} - r_1 A_1 - \kappa _{A1}'V_1A_1 \end{aligned}$$
(22f)

The solution of Eq. 22 for the number of target cells and amount of type I interferon is shown in Fig. 3 for the parameters in Table 2, which are adapted from Cao et al. (2015) and will be used for the remainder of the study. The time dependence of the extinction probability of the second virus is then entirely due to these two functions.

Fig. 3
figure 3

The number of target cells, T(t), and the amount of type I inferferon, F(t)

Full size image
Table 2 Parameter values for the two-strain model
Full size table

The second virus and infected cell equations can be written in the form of the TIV model with time-dependent parameters which are functions of the two time courses shown in Fig. 3:

$$\begin{aligned} \frac{dV_2}{dt}&= \frac{p_{V2}}{1+s_2F}I_2 - (\delta _{V2} + \beta _2 T) {V_2} \nonumber \\&= p(t)I_2 - (\delta _{V2} + \beta _2T(t))V_2 \nonumber \\ \text {where} \qquad p(t)&\equiv \frac{p_{V2}}{1+s_2F(t)} \end{aligned}$$
(23a)
$$\begin{aligned} \frac{dI_2}{dt}&= \beta _2 V_2T(t) - \delta _{I2}I_2. \end{aligned}$$
(23b)

We substitute the time courses T(t) and p(t) into Eq. 20. \(\beta (t)\), \(\delta (t)\) and c(t) in Eq. 20 are constants, equal to \(\beta _2\), \(\delta _{I2}\) and \(\delta _{V2}\) in Eq. 23 respectively. We solve the coupled ODEs in Eq. 20 for the extinction probability when one virion of \(V_2\) is introduced at time t after the introduction of a fixed number of virions of \(V_1\), then use Eq. 6 to calculate the extinction probability when \(n_V\) virions of \(V_2\) are introduced. To solve the ODEs, we use our knowledge that as \(t \rightarrow \infty \), \(T(t) \rightarrow T_0\) and \(p(t) \rightarrow p_{V2}\), so \(\rho _V(t)\) and \(\rho _I(t)\) approach the values given in Eq. 5. Thus, we set the initial values

$$\begin{aligned} \rho _V(t_e)&= \mathrm {min}\left( 1-\frac{R_{02} - 1}{\alpha _2},1\right) \end{aligned}$$
(24a)
$$\begin{aligned} \rho _I(t_e)&= \mathrm {min}(1/R_{02},1), \end{aligned}$$
(24b)

where \(t_e\) is a predetermined time much further than the end of the primary infection such that \(T(t_e)/T_0\), \(p(t_e)/p_{V2}\) and \(\delta (t)/\delta _{I2}\) are all approximately equal to one. The subscript 2 indicates that \(R_0\) and \(\alpha \) are to be calculated using the parameters of the second virus.

The ODEs are solved using MATLAB R2014b’s ode15s ODE solver (The MathWorks, Natick, MA), with absolute tolerance of \(10^{-12}\) and relative tolerance of \(10^{-6}\). We set \(t_e = 50\), which is much longer than the time course of a single infection. For numerical stability of the ODE solver, we substitute \(t' = -t\), then integrate from \(t' = -t_e\) to \(t' = 0\). We validate the results against stochastic simulations conducted using Gillespie’s tau-leap algorithm (Gillespie 2001) conducted with timestep \(dt = 10^{-3}\).

Fig. 4
figure 4

Calculated versus simulated extinction probabilities when introducing ten virions of the second strain a number of days after primary inoculation (the inter-exposure interval), indicated on the x-axis. Parameters are as per Table 2. The simulated probabilities are obtained over 1000 simulations for each scenario. The grey area indicates the 95 % prediction interval for the proportion of extinction events out of 1000, assuming that the number of extinction events follows a binomial distribution where the event probability is the calculated extinction probability

Full size image

Figure 4 shows the calculated versus simulated extinction probabilities when introducing ten virions of the second strain a number of days after primary inoculation (the inter-exposure interval), as indicated on the x-axis, with the time-dependent parameters in Fig. 3. The simulated probabilities are obtained over 1000 simulations for each scenario. The grey area indicates the 95 % prediction interval for the proportion of extinction events out of 1000, assuming that the number of extinction events follows a binomial distribution where the event probability is the calculated extinction probability. We get excellent agreement between the calculated and simulated results.

Fig. 5
figure 5

Calculated versus simulated extinction probabilities when introducing different numbers of virions of the second strain (as indicated on the x-axis) with a 3.5-day inter-exposure interval. Parameters are as per Table 2. The simulated probabilities are obtained over 1000 simulations for each scenario. The grey area indicates the 95 % prediction interval for the proportion of extinction events out of 1000, assuming that the number of extinction events follows a binomial distribution where the event probability is the calculated extinction probability

Full size image

Figure 5 shows the calculated versus simulated extinction probabilities when introducing different numbers of virions of the second strain (as indicated on the x-axis) with a 3.5-day inter-exposure interval. We get excellent agreement between the calculated and simulated results over a wide range of initial virion numbers.

We can modify the model to include K latent stages and M infectious stages for the cells infected with the second virus. Although one may be biologically motivated to do the same for the first virus, for ease of comparison we leave our model for the first virus unchanged with no latent stages and one infectious stage for the first virus, such that T(t) and p(t) are the same as for the previous case. Equation 23 then becomes [with the same definitions of T(t) and p(t)]

$$\begin{aligned} \frac{dV_2}{dt}&= p(t)\sum _{m=1}^{M}I_{m2} - (\delta _{V2} + \beta _2T(t))V_2 \end{aligned}$$
(25a)
$$\begin{aligned} \frac{dL_{12}}{dt}&= \beta _2 V_2T(t) - K\delta _{L2}L_{12} \end{aligned}$$
(25b)
$$\begin{aligned} \frac{dL_{k2}}{dt}&= K\delta _{L2}L_{k-1,2} - K\delta _{L2}L_{k2}, \quad k = 2,\ldots ,K-1 \end{aligned}$$
(25c)
$$\begin{aligned} \frac{dL_{K2}}{dt}&= K\delta _{L2}L_{K-1,2} - K\delta _{L2}L_{K2} \end{aligned}$$
(25d)
$$\begin{aligned} \frac{dI_{12}}{dt}&= K\delta _{L2}L_{K2} - M\delta _{I2}I_{12} \end{aligned}$$
(25e)
$$\begin{aligned} \frac{dI_{m2}}{dt}&= M\delta _{I2}I_{m-1,2} - M\delta _{I2}I_{m2}, \quad m = 2,\ldots ,M-1 \end{aligned}$$
(25f)
$$\begin{aligned} \frac{dI_{M2}}{dt}&= M\delta _{I2}I_{M-1,2} - M\delta _{I2}I_{M2}. \end{aligned}$$
(25g)

The initial values at \(t = t_e\) are given by Eq. 13, where all parameter values are for the second virus. We will investigate how changing the mean latent and infectious periods of the second virus, \(1/\delta _{L2}\) and \(1/\delta _{I2}\), and the number of latent and infectious stages, K and M, changes the extinction probability. All other parameters will be fixed, and are listed in Table 2.

First, we assume that for cells infected with the second virus, there is one latent stage, and one infectious stage with \(\delta _{I2} = 3\).

Fig. 6
figure 6

Extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different values of \(1/\delta _{L2}\) as indicated by the legend. Parameters are as per Table 2

Full size image

Figure 6 shows the calculated extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different values of \(1/\delta _{L2}\) as indicated by the legend. Parameters are as per Table 2. We see that for this particular model and set of initial conditions, if the second virus is introduced in the early stage of infection, increasing \(1/\delta _L\) increases the extinction probability. This is because increasing the mean latent period increases the chance that a cell infected with the second virus soon after introduction becomes infectious when the innate immune response is most active, so the expected number of virions it produces decreases. However, if the second virus is introduced in the middle stage of the first infection (2–5 days), when the innate immune response is most active, then increasing \(1/\delta _{L2}\) decreases the extinction probability. This is because increasing \(1/\delta _{L2}\) increases the chance that a cell infected with the second virus soon after introduction becomes infectious when the innate immune response is much less active, so the expected number of virions it produces increases. Moreover, as the mean latent period approaches 0, the extinction probability curve approaches that of the model without a latent stage, as shown in Fig. 4. A comparison with simulation results for this figure and subsequent figures is given in Online Resource 1, and agreement is excellent.

Fig. 7
figure 7

Extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different values of \(1/\delta _{I2}\) as indicated by the legend. Parameters are as per Table 2, with \(\delta _{L2} = 3\) and one latent stage

Full size image

By contrast, Fig. 7 shows the calculated extinction probabilities when \(1/\delta _{I2}\) is changed. ten virions of the second strain are introduced with varying inter-exposure intervals, for different values of \(1/\delta _{I2}\) as indicated by the legend. Parameters are as per Table 2, with \(\delta _{L2} = 3\) and one latent stage. In contrast to increasing \(1/\delta _{L2}\), increasing \(1/\delta _{I2}\) always decreases the extinction probability. This is because increasing \(1/\delta _{I2}\) increases the basic reproduction number \(R_{02}\) by increasing the expected number of virions produced by an infected cell (\(p_2/\delta _{I2}\)) in the absence of the immune response. In the presence of the immune response, increasing \(1/\delta _{I2}\) still increases the expected number of virions produced by a cell which becomes infected at time t, \(\displaystyle \int _{t}^{\infty } p_2/(1+s_2F(t)) \exp (-\delta _{I2}(t'))\,dt'\).

Fig. 8
figure 8

Extinction probabilities when introducing ten virions of the second strain with a 3.5-day inter-exposure interval, for different numbers of \(1/\delta _{L2}\) and \(1/\delta _{I2}\) as indicated by the x- and y-axes, when there is one latent stage and one infectious stage. Parameters are as per Table 2

Full size image

Figure 8 shows the extinction probabilities when introducing ten virions of the second strain with a 3.5-day inter-exposure interval, for different values of \(1/\delta _{L2}\) and \(1/\delta _{I2}\) as indicated by the x- and y-axes, when there is one latent stage and one infectious stage. We see that as \(1/\delta _{I2}\) is increased, for the extinction probability to be conserved, \(1/\delta _{L2}\) must be decreased. \(1/\delta _L\) and \(1/\delta _I\) cannot both be determined by measurements of extinction probability at one time point alone, but when considered with additional data, such as viral load data, the degeneracy in parameter values is less than if each type of data were considered separately. Measurements of extinction probability for different inter-exposure intervals, then fitting parameters to be consistent with the data at all time points, is another way of decreasing the degeneracy.

We now look at how changing the number of latent and/or infectious stages changes the extinction probability. First, we consider a model with no latent stages and M infectious stages. We hold \(\delta _{I2}\) constant at 3.

Fig. 9
figure 9

Extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different numbers of infectious stages as indicated by the legend. Parameters are as per Table 2, and \(\delta _{I2}\) is held constant at 3

Full size image

Figure 9 shows the extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different numbers of infectious stages as indicated by the legend. We see that for short inter-exposure intervals (0–1 days) and long inter-exposure intervals (6 days or more), increasing the number of infectious stages decreases the extinction probability, but for medium inter-exposure intervals (2–5 days), increasing the number of infectious stages increases the extinction probability. For long inter-exposure intervals, the number of target cells at the time of introduction is close to \(T_0\) and the immune response F is close to zero, so the dependence of the extinction probability on the number of infectious stages is like that of the TIV model in Fig. 1. For short inter-exposure intervals, the initial conditions suggest that exponential growth of infected cells and virions is likely (as target cells have not been depleted and the immune response is not yet active), so if extinction has not occurred by the time the temporary depletion of target cells and the immune response occur, the number of virions and infected cells will have grown large enough to survive the immune response. Hence, the dependence of the extinction probability on the number of infectious stages is once again like that of the TIV model. For medium inter-exposure intervals, it appears that the dependence of the extinction probability on the number of infectious stages depends on the parameter functions. For example, if we change \(q_1\) to \(5 \times 10^{-6}\) in the model, thus changing T(t) and F(t) such that the immune response is stronger and target cells do not become as depleted, we see in Fig. 10 that increasing the number of infectious stages decreases the extinction probability for all inter-exposure intervals. Compared to changing the mean infectious period, changing the number of infectious stages only has a small effect on the extinction probability for our model and parameters; this is because the chosen mean latent and infectious periods, which lie within the biologically realistic parameter ranges for influenza, are short compared to the timeframe within which the immune response changes significantly, so the distributions of latent and infectious periods do not have a large effect on the immune response experienced by an infectious cell. However, since the biologically realistic parameter ranges for the mean latent and infectious period are large, changing their values within this range has a large effect on the immune response experienced by an infectious cell.

Fig. 10
figure 10

Extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different numbers of infectious stages as indicated by the legend. Parameters are as per Table 2 except \(q_1 = 5\times 10^{-6}\), and \(\delta _{I2}\) is held constant at 3

Full size image

Finally, we consider a model with K latent stages and one infectious stage. Both \(\delta _{L2}\) and \(\delta _{I2}\) are held constant at 3.

Fig. 11
figure 11

Extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different numbers of latent stages as indicated by the legend. Parameters are as per Table 2, and both \(\delta _{L2}\) and \(\delta _{I2}\) are held constant at 3

Full size image

Figure 11 shows the extinction probabilities when introducing ten virions of the second strain with varying inter-exposure intervals, for different numbers of latent stages as indicated by the legend. As the number of latent stages increases, the extinction probability increases slightly for short inter-exposure intervals, but decreases slightly for long inter-exposure intervals (the difference is too small to be seen in the figure). Increasing the number of latent stages decreases the variance in the time taken for an infected cell to become infectious. As a result, for short inter-exposure intervals, the first infected cells are more likely to become infectious when the immune response is strong, rather than before or after the time at which there is a strong immune response. For long inter-exposure intervals, the immune response exponentially decays, so decreasing the variance in the time the first infected cells become infectious decreases the mean immune response experienced by the infectious cells. Overall, for this particular model of the immune response and initial values, the difference in extinction probability for different numbers of latent and infectious stages is small, suggesting that the number of latent and infectious stages may be difficult to determine from experimental measurements of extinction probability. However, the number of latent and infectious stages chosen will affect the viral load time course given a set of parameters, so when the extinction probability data is combined with other data such as viral load measurements, the number of latent and infectious stages chosen will impact parameter estimates.

4 Discussion

In this study, we have calculated the extinction probability of a disease upon introduction of a number of virions and/or infected cells to the host. Previous studies have derived these expressions for the TIV model, both with time-independent parameters (Pearson et al. 2011) and time-dependent parameters (Conway et al. 2013) which reflect the effects of prophylaxis on viral reproduction. We have generalised these expressions to models where cells have multiple latent and/or infectious stages, such that the times spent in the latent and infectious stages are Erlang-distributed. We have also shown how mechanistic models of the immune response can be written in terms of time-dependent viral parameters, such that the extinction probability given the number of initial virions and/or infected cells can be calculated.

For a time-dependent immune response, we have shown that the extinction probability depends on the time of introduction and the initial number of virions; these results are verified by full stochastic simulations, and excellent agreement is obtained. We have also demonstrated how the extinction probability changes when the mean latent and infectious periods are changed, and when the number of latent and infectious stages are changed. This suggests that the extinction probability at different times of introduction can be used to help determine the parameters of the model. Because the mapping of extinction probability to parameter values is not one-to-one, additional information, such as the viral load time course, is also needed.

The methods used in this study are applicable to many situations where the viral dynamics cannot be described with a simple TIV model. One of these situations, which we have used as a case study, is when the host has prior immunity due to a previous infection. We have shown that given a time-dependent immune response induced by the first virus, changing the introduction time of the second virus changes the probability that infection is established. This helps explain experimental observations by Laurie et al. (2015) that immunity conferred by a primary infection against subsequent infection is temporary, and is only observed in a subset of hosts under identical experimental conditions.

Other within-host influenza studies have constructed different models incorporating the effects of the immune response (Beauchemin and Handel 2011; Smith and Ribeiro 2010); the viral loads predicted by these models have been compared qualitatively (Dobrovolny et al. 2013), and different models have been fitted to the same viral load data to determine the best-fitting model (Pawelek et al. 2012). However, the effectiveness of the immune responses in these models in preventing subsequent infection have not been compared systematically. Our approach to calculating the extinction probability given a mechanistic model provides the tools to examine how different modelled components of the immune response work together to prevent subsequent infection. The emergence of drug-resistant viral strains during prophylaxis is an increasingly important issue, and several models, both deterministic and stochastic, have been developed to study the conditions for the emergence and survival of these strains (Handel et al. 2007; Canini et al. 2014; Haeno and Iwasa 2007; Moghadas 2011), but only some of these models take the immune response into account. Extinction probability calculations enable the study of how both stochasticity and a time-dependent immune response affect the proliferation of drug-resistant strains.

The effect of prophylaxis on the viral load had been studied using various deterministic models (Baccam et al. 2006; Beauchemin et al. 2008; Canini et al. 2014; Dobrovolny et al. 2011, 2013); our methods can be used to evaluate how different pre- and post-exposure drug prophylaxis regimes not only reduce the severity of infection, but prevent infection altogether. Conway et al. (2013) analysed how pre- and post-exposure prophylaxis affects the probability of HIV infection becoming established. In their model, prophylaxis has an explicit time-dependent effect on the infectivity of virus; we have shown how these time-dependent effects can be generated by mechanistic models.

A limitation of the extinction probability calculation method used in our study is that it can only calculate the probability of extinction given that the time-dependent immune response is independent of viral and infected cell dynamics. For example, some models predict that the viral load oscillates before reaching a steady state, such as due to delayed cellular immune responses (Canabarro et al. 2004); one would expect the extinction probability to be largest at these troughs. In these cases, extinction is due to target cell depletion and/or immune responses which vary greatly with the stochastic time courses of the virions and infected cells. Because of this feedback loop, the parameters in Eq. 20 are no longer functions known a priori, and the method cannot be used. Also, there are experimental difficulties in obtaining the extinction probability because extinction can only be observed once in each replica of the experiment. This makes it difficult to determine the effect of inter-host differences on the obtained extinction probability across hosts.

Moreover, when estimating parameters using extinction probability data, the initial number of virions must either be known or be included as an estimated parameter. The number of virions transmitted may be fewer than the number which reach the infected body parts, and it is often impractical to measure the viral load immediately after transmission. As a result, the initial viral load is often an unknown parameter which is fitted to the data. Even if the initial viral load can be thus estimated accurately, viral load is usually measured in units such as PFU, TCID\(_{50}\) or EID\(_{50}\), and as a sample concentration rather than as the total number of virions in the host. There can be much uncertainty in the conversion factor between measurement units and the number of virions (Handel et al. 2007). If we could experimentally determine the conversion factor between measurement units such as TCID\(_{50}\) and virions, this would enable a better understanding of the impact of extinction on infection outcomes.

Our study has looked at the effect of gamma-distributed latent and infectious periods, but some models have used other distributions, such as a constant infectious/latent period, a normal distribution or a log-normal distribution (Holder and Beauchemin 2011). In this case, the inter-event times are no longer Markovian, and the extinction probability from a given state is no longer only dependent on the number of virions and infected cells, but also on their ages. Hence, the method used in our study cannot be directly applied in these cases. Extension to general distributions for latent and infectious periods is the subject of future work.

The methods used in our study can also be used in an epidemic context; for example, where two viruses are co-circulating in the population and infection with one confers immunity to the other, depletion of susceptibles by one virus affects the extinction probability when another virus is introduced into the population. This is analogous to the within-host two-strain model presented in Sect. 3.2. Exploration of the application of our study to epidemic contexts is left for future research.