Extinction in an HIV model

Via Professor Jessica Conway’s really nice paper and complete notes in her supplementary file, I’m redoing a calculation for my own model with my own notation so that I can see how it works.

First let me write down the full set of ordinary differential equations in a mean field picture as described in Intra-host HIV model with latency.

$\dot{S}=\alpha_S-\delta_S S - \beta S V\\ \dot{L}=\theta_L L + (1-\tau)\beta S V\\ \dot{A}=\xi L -\delta_A A + \tau\beta S V\\ \dot{V}=\pi A - \gamma' V$

We are thus considering latent and actively infected cells that become active with probability $\tau$ . Birth rates are given by $\alpha$ s and deaths by $\delta$ s. Activation is given by $\xi$ and the remaining variables are typical but expressed using greek letters.

Note, we can be clear on where the viruses all go if we use $\gamma'=\gamma+\beta S$ then, using the quasi equilibrium $V^*(t)=\pi A/\gamma'$ . If the latent pool is disregarded, the basic reproduction number $R_0$ is calculated when the population is fully susceptible, i.e. $S^*=\alpha_S/\delta_S$ , so we solve the equation

$\dot{A}= \delta_A(\frac{\tau\beta S^* \pi}{\delta_A\gamma'}-1)$

which gives the solution, defining $R_0=\frac{\tau\beta S^* \pi}{\delta_A\gamma'}$

$A= A_0e^{\delta_A(R_0-1)t}$

which goes to zero precisely when $R_0<1$ . The result then is when the basic reproduction number is less than one, the virus will always die out with probability equal to one.

In the situation where the virus does die out, calculating the generations until the die out is useful. We formulate the branching process stochastic analogue to the ODEs above considering only activated and latent cells. We assume the original number of active cells is then $A_0 = 1$ . The active cell can then make virions which go on to infect other cells. Thus we have the possible transitions:

$A \xrightarrow{\pi} A+V\\ V \xrightarrow{\tau\beta S^*} A\\ A \xrightarrow{\delta_A} 0\\ V \xrightarrow{\gamma} 0$

assuming that the number of susceptible cells is at equilibrium $S^*=\alpha_S/\delta_S$ .

The probability of $A$ cells arising from $N$ viruses assuming a single virus infects a single cell is then

$p_A = {N \choose A} \left( \frac{\tau\beta S^*}{\gamma + \tau\beta S^*}\right)^A\left(\frac{\gamma}{\gamma+\tau\beta S^*}\right)^{N-A}$ .

Probability generating functions

With the probability, we can then write the probability generating function (pgf, as described in my post: Adding beans to a pot) for the probability of active cells as

$G(s) = \sum_{A=0}^N p_A s^A$

putting in the probability $p_A$ , we can group the terms with exponent $A$ , and use the binomial theorem $(x+a)^n = \sum_{k=0}^\infty {n \choose k} x^k a^{n-k}$ . This leads to

$G(s) = \sum_{A=0}^N {N \choose A} \left( \frac{\tau\beta S^*s}{\gamma + \tau\beta S^*}\right)^A\left(\frac{\gamma}{\gamma+\tau\beta S^*}\right)^{N-A} = \left( \frac{\tau\beta S^*s + \gamma}{\gamma + \tau\beta S^*}\right)^N$

so that defining $r=1+ \frac{\gamma}{\tau\beta S^*}$ we have

$G(s) = \frac{1}{r^N} (s + r^2)^N$

We might take $N$ to also be a random variable because active cells can produce varying numbers of virions, $N$ . The probability distribution $p_N$ arises from the sequence of events that leads to $N$ virions before the cell dies. Thus the probability distribution is the product of $N$ birth events having normalized probability $\frac{\pi}{\pi+\delta_A}$ and a single event with probability $\frac{\delta_A}{\pi+\delta_A}$ :

$p_N = \left( \frac{\pi}{\pi+\delta_A}\right)^N\left(\frac{\delta_A}{\pi+\delta_A}\right)$

We can then take the average generating function with respect to the probability of $N$ , in the discrete sense, this is a sum over all $N$ multiplied by $p_N$ so that we have

$\sum_{N=0}^\infty p_N G(s) = \sum_{N=0}^\infty \left( \frac{\pi}{\pi+\delta_A}\right)^N\left(\frac{\delta_A}{\pi+\delta_A}\right) \left( \frac{\tau\beta S^*s + \gamma}{\gamma + \tau\beta S^*}\right)^N$

so we can group terms with exponent $N$ , and then use that the sum of the geometric series $\sum_{n=0}^\infty x^n = \frac{1}{1-x}$ for $latex|x| \leq 1$. This leads to the average generating function with respect to $N$ , i.e. $\langle \cdot \rangle_N$ :

$\langle G(s) \rangle_N = \left(\frac{\delta_A}{\pi+\delta_A}\right) \left[1-\left( \frac{\tau\beta S^*s + \gamma}{\gamma + \tau\beta S^*} \frac{\pi}{\pi+\delta_A} \right)\right]^{-1}$

or pulling the terms into the denominator

$\langle G(s) \rangle_N = \left[\frac{\pi+\delta_A}{\delta_A}-\frac{\tau\beta S^*s + \gamma}{\gamma + \tau\beta S^*} \left(\frac{\pi}{\delta_A} \right)\right]^{-1}$

Now factoring

$\langle G(s) \rangle_N = \frac{\gamma+\tau\beta S^*}{\frac{\pi(\gamma + \tau\beta S^*)}{\delta_A} + \gamma + \tau\beta S^* -\frac{\pi\tau\beta S^*}{\delta_A}s - \frac{\gamma \pi}{\delta_A} }$

we can cancel the first and last terms in the denominator

$\langle G(s) \rangle_N = \frac{\gamma+\tau\beta S^*}{\frac{\pi\tau\beta S^*}{\delta_A} + \gamma + \tau\beta S^* -\frac{\pi\tau\beta S^*}{\delta_A}s }$

factoring again we have

$\langle G(s) \rangle_N = \frac{1}{\frac{\pi\tau\beta S^*}{\delta_A( \gamma + \tau\beta S^*)}(1-s) +1 }$

Conway defines $R_c = \frac{\pi\tau\beta S^*}{\delta_A( \gamma + \tau\beta S^*)}$ so that the average generating function is finally

$\langle G(s) \rangle_N = \frac{1}{R_c(1-s) +1 }$

Extinction process

We now have the expression for the generating function

$\langle G(s) \rangle_N = \frac{1}{R_c(1-s) +1 }$

If we know an extinction will happen, that is, $R_0$ is less than 1, we have a subcritical branching process. We can consider the probability of an infected cell producing $0,1,2,\ldots N$ offspring as $p_0, p_1, p_2,\ldots p_N$ , where we assume a maximum amount does exist.

The probability of ultimate death (extinction) in the $m$ th generation we call $d_m$ . For $d_m=1$ , or extinction in the $m$ th generation, requires everything to die after the $(m-1)$ th generation. The probability of $j$ cells dying in generation $m-1$ is $(d_{m-1})^j$ . For all birth events to die out then we must balance each birth probability with the equivalent death probability, i.e.:

$d_m = p_0 + p_1(d_{m-1}) + p_2(d_{m-1})^2 + \ldots = \sum_{k=0}^N p_k (d_{m-1})^k$

for which the rhs is precisely the pgf we defined before. That means we have the condition

$d_m = G(d_{m-1})$

Now let’s try to solve this relationship for the definitions of $G(s)$ given above, identifying $s=d_{m-1}$ . Conway gives us $d_k=\frac{1-R^{k+1}}{1-R^{k+2}}$ , where I’ve dropped the subscript. Let’s check that with the average pgf:

First $d_{k-1} = \frac{1-R^{k}}{1-R^{k+1}}$

$\frac{1-R^{k+1}}{1-R^{k+2}} = \frac{1}{R(1-\frac{1-R^{k}}{1-R^{k+1}}) +1}$

cross multiplying

$(1-R^{k+1}) (R(1-\frac{1-R^{k}}{1-R^{k+1}}) +1) = 1-R^{k+2}$

totally writing all terms out

$R-R^{k+2}- R+R^{k+1} + 1-R^{k+1} = 1-R^{k+2}$

we can see immediately that the equality is true. Therefore, we have derived the probability of dying out in the k-th generation.

$d_k=\frac{1-R^{k+1}}{1-R^{k+2}}$

Note this result only works for $R<1$ otherwise we could have negative probabilities. In this “subcritical” regime, we can be interpreted the expression as the probability of not extinction in the next generation relative to probability of not-extinction in some generation.

Extinction in an HIV model

Probability generating functions

Extinction process

Published by dbrvs

Leave a Reply Cancel reply

Probability generating functions

Extinction process

Share this:

Like this:

Related

Published by dbrvs

Leave a Reply Cancel reply