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 \alphas and deaths by \deltas. 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 mth generation we call d_m. For d_m=1, or extinction in the mth 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.


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.


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s