The Coalescence of Intrahost HIV Lineages Under Symmetric CTL Attack

Abstract

Cytotoxic T lymphocytes (CTLs) are immune system cells that are thought to play an important role in controlling HIV infection. We develop a stochastic ODE model of HIV–CTL interaction that extends current deterministic ODE models. Based on this stochastic model, we consider the effect of CTL attack on intrahost HIV lineages assuming that CTLs attack several epitopes with equal strength. In this setting, we introduce a limiting version of our stochastic ODE under which we show that the coalescence of HIV lineages can be described through Poisson–Dirichlet distributions. Through numerical experiments, we show that our results under the limiting stochastic ODE accurately reflect HIV lineages under CTL attack when the HIV population size is on the low end of its hypothesized range. Current techniques of HIV lineage construction depend on the Kingman coalescent. Our results give an explicit connection between CTL attack and HIV lineages.

Appendix

1.1 A.1 Proof of Lemma 5.2

In this section, we provide the technical details that support Conclusions 3 and 4 of Lemma 5.2. From Conclusions 1 and 2, we can reduce (12) to the following:

(59)

The O(μ) terms in the last two equations directly above can be ignored because $T_{i}^{h} - T_{i} = O(|\log(\delta)|)$ and μ|log(δ)|→0. Dropping these terms, we note the following relation:

$$ \frac{d (\log(e_i) - \log(e_{i-1}) )}{dt} = \Delta k. $$

(60)

Integrating the above equation and using our assumptions on e _i−1(T _i ) and e _i (T _i ), we find

$$ \frac{e_i(t)}{e_{i-1}(t)} = O_p\bigl( \delta \exp\bigl[\Delta k(t - T_i)\bigr]\bigr). $$

(61)

If we can show that e _i is bounded, then as t grows, (61) implies that e _i−1 collapses. To see that e _i is bounded, set

$$ z(t) = \frac{h(t)}{g} + \frac{e_i(t)}{\gamma} + \frac{e_{i-1}(t)}{\gamma}. $$

(62)

Then by straightforward differentiation,

(63)

Since z(t) is nonnegative, we find that z(t) must be bounded. In turn, h, e _i , and e _i−1 must be bounded. Returning to (61) and setting t≥T _i +2/Δk|log(δ)|, we find, since e _i is bounded,

$$ e_{i-1}(t) = O_p(\delta), $$

(64)

and this gives Conclusion 4.

Once t>2/Δk|log(δ)|, we can further reduce (59) to

(65)

Consider then (65). Ignoring the O _p (δ) term for a moment, the system is not dependent on δ. Since we have shown h,e _i to be bounded, application of the Poincaré–Bendixon theorem shows that the system converges to its nontrivial equilibrium, $h = \frac{k_{i}}{\gamma}$ and $e_{i} = \frac{1-h}{h}$. Now consider the O _p (δ) term. Given some fixed distance ϵ>0, if we run the system from t=2/Δk|log(δ)| to t=3/Δk|log(δ)|, we are guaranteed by choosing δ sufficiently small to be within ϵ of the equilibrium. In turn, taking ϵ small, we can linearize (65) about its equilibrium.

Straightforward computation shows that both eigenvalues of the linearized system have negative real part bounded above by

$$ \rho = -\bigl(g\gamma + O_p(\delta) \bigr) \min\biggl(1, \frac{4k^2}{g\gamma}\biggl(1 - \frac{k}{\gamma}\biggr) \biggr). $$

(66)

Running the system from t=3/Δk|log(δ)| to $t = (3 + \frac{2}{|\rho |})|\log(\delta)|$ forces (65) to within O _p (δ) of the equilibrium. This gives Conclusion 3.

1.2 A.2 Proof of Propositions 5.4 and 5.5

We first prove Proposition 5.4. Consider the Laplace transform of I(T)/ω,

$$ E\biggl[\exp\biggl[-\lambda \frac{I(T)}{\omega}\biggr]\biggr] = \exp[-\psi], $$

(67)

where

$$ \psi = \int_{-\infty}^{T} dt \mu \exp[r_1 t] \biggl(1 - E[\exp\biggl[-\frac{\lambda}{\omega} B(T - t)\biggr] \biggr), $$

(68)

and where B(t) is a continuous-time branching process with birth and death rates b and d run for time t and given B(0)=1.

The formula for E[exp[−λB(t)/ω]] is well known (Athreya and Ney 1972). Plugging into this formula and shifting t leads to

(69)

Making the substitution s=(b/r ₂)exp[r ₂ t]/ω gives

$$ \psi = A_0 \int_{\frac{b}{r_2 \omega}}^\infty ds \frac{1}{s^{1+\alpha}} \biggl(1 - \frac{1}{1 + \lambda s}\biggr) \biggl(1 + O\biggl(\frac{1}{\omega}\biggr)\biggr), $$

(70)

where α=r ₁/r ₂ and

$$ A_0 = \biggl(\frac{b}{r_2}\biggr)^\alpha \frac{1}{b}. $$

(71)

The expression 1/(1+λs) in (70) can be recognized as the Laplace transform of an exponential. Substituting the density of an exponential and taking the ω→∞ limit leads to

$$ \psi \to A_0 \int_{0}^{\infty} ds \frac{1}{s^{1+\alpha}} \int_0^\infty dz \bigl(1 - \exp[-\lambda s z]\bigr) \exp[-z]. $$

(72)

Applying the Fubini theorem, we arrive at

$$ \psi \to A_1 \int_{0}^{\infty} \frac{\alpha}{\varGamma(1 - \alpha)} \,ds \frac{1}{s^{1+\alpha}} \bigl(1 - \exp[-\lambda s]\bigr), $$

(73)

where Γ is the gamma function, and

$$ A_1 = \biggl(\frac{b}{r_2}\biggr)^\alpha \frac{1}{b} \varGamma(1 + \alpha) \varGamma(1 - \alpha) \frac{1}{\alpha}. $$

(74)

All of the above gives

$$ E\bigl[-\lambda I(T_\text {sample})/\omega\bigr] \to \exp\bigl[-A_1| \lambda|^{\frac{r_1}{r_2}}\bigr]. $$

(75)

The Laplace transform above is that of $A_{1}^{r_{2}/r_{1}} \mathcal {S}(r_{1}/r_{2}, 1, 1)$ (Nolan 2011; Bertoin 1996), which demonstrates Proposition 5.4.

Proposition 5.5 follows almost directly from the arguments that gave (73). By the uniqueness of the Levy measure, $I(T_{\text {sample}})/\omega$ converges to a nonhomogeneous Poisson process with Levy measures proportional to $s^{-1-r_{1}/r_{2}}$. In other words, the distribution of the jumps that form I(T)/ω converge to those of a stable process of index r ₁/r ₂. Each such jump corresponds to the descendants of an immigrant. Partitioning according to immigrants is then given by PD(r ₁/r ₂,0) by existing results (Pitman 2002).

1.3 A.3 Proofs of Conclusion 2 of Lemma 5.3 and Lemma 5.6

We need to frame Propositions 5.4 and 5.5 in the context of spawning period dynamics under the SPL. In the period [T _i ,T _i+1], we consider “immigrants” representing v _i+1→v _i+2 mutations. We have the following facts that relate to the assumptions of Proposition 5.4:

• 1. The number of v _i+1 variants at t is $\mathbb {E}\delta \exp[-(\Delta k+ O(\delta))(T_{i+1} - t)]$.

• 2. v _i+1→v _i+2 mutations occur at rate $\mu \mathbb {E}\delta \exp[-(\Delta k+ O(\delta))(T_{i+1} - t)]$.

• 3. v _i+2 variants produced descendants according to a continuous-time binary branching process with birth and death rates k _i +O(δ) and k _i+2+O(δ), respectively.

The O(δ) terms in the rates above will not affect the SPL limits and can be ignored. Essentially, under the SPL, the effects of O(δ) perturbations in the rates will disappear in the integral (69). If we set T=T _i+1 and $\omega = (\mu \mathbb {E}k_{i} \delta)^{2}$, then we may apply Proposition 5.4 to the spawning period except that “immigrants” arrive only after time T _i . This does not impact the ψ limit in the proof of Proposition 5.4. Indeed, (70) is altered to

$$ \psi = A_0 \int _{\frac{b}{r_2 \omega}}^{\frac{b\exp[r_2(T_{i+1}-T_i)}{r_2 \omega}} ds \frac{1}{s^{1+\alpha}} \biggl(1 - \frac{1}{1 + \lambda s}\biggr) \biggl(1 + O\biggl( \frac{1}{\omega}\biggr)\biggr), $$

(76)

where r ₂=2Δk, r ₁=Δk, and α=1/2. Under the SPL, ω→∞, so the lower bound of integration in (76) goes to 0. For the upper bound, note that

$$ \exp\bigl[r_1(T_{i+1} - T_i)\bigr] = \frac{e_{i+1}(T_{i+1})}{e_{i+1}(T_i)} = \frac{\delta}{O(\mu^2 \mathbb {E}\delta^2)}, $$

(77)

giving

(78)

So the upper bound of integration in (76) goes to ∞, and Proposition 5.5 applies. Lemma 5.6 now applies by the same arguments as Proposition 5.5.

Plugging in the appropriate values for r ₁,r ₂,b and considering the k _i factor in ω, we find

$$ \frac{e_{i+2}(T_{i+1})}{\mu^2 \mathbb {E}\delta^2} \to \biggl(\frac{k_i}{\Delta k} \frac{\pi^2}{2} \biggr) \mathcal {S}\biggl(\frac{1}{2}, 1, 1\biggr). $$

(79)

This gives Conclusion 2 of Lemma 5.3.

1.4 A.4 Justification of Deterministic Approximations

To justify our deterministic approximations, we describe the steps needed to make the results of Sect. 5.1 rigorous. Associated with the time interval [T _i ,T _i+1] we define $h^{\text {det}}$ and $e_{j}^{\text {det}}$ for j≤i+2 as follows:

(80)

For j≤i,

$$ \frac{\text{d}e_j^\text {det}}{\mathrm{d}t} = e_j^\text {det}h^\text {det}- k_j e_j^\text {det}+\mu \gamma e_{j-1}^\text {det}h^\text {det}, $$

(81)

and

(82)

(83)

System (80)–(82) reflects the assumptions we make in the proofs of Lemmas 5.2 and 5.3. Namely, e _j dynamics for j≥i are taken as deterministic, e _i+1 dynamics are taken as deterministic only after $T_{i}^{h}$, and e _i+2 dynamics are analyzed stochastically on all of [T _i ,T _i+1].

Let $\Delta h = h - h^{\text {det}}$ and $\Delta e_{j} = e_{j} - e_{j}^{\text {det}}$ for j≤i+2. We assume that Δh(T _i )=Δe _j (T _i )=0 for all j. Set

(84)

The following lemma shows that system (80)–(82) behaves very similarly to the fully stochastic system (12).

Lemma A.1

$$ P\biggl(\sup_{t \in [T_i, T_{i+1}]} \big|v(t)\big| > \frac{1}{\mathbb {E}^{\frac{1}{8}}}\biggr) \to 0 $$

(85)

and

$$ \sup_{t \in [T_i, T_{i+1}]} \bigg|\frac{e_{i+1}(t)}{e_{i+1}^{\mathrm{det}}(t)} - 1\bigg| \to 0. $$

(86)

The results of Sect. 5.1 can be made rigorous by using Lemma A.1 to make straightforward modifications of the proofs. For example, to demonstrate Conclusion 2 of Lemma 5.3, we need to modify the arguments made in Sect. A.2. More precisely, instead of considering (68), we consider

$$ \psi = \int _{T_i}^{T_{i+1}} dt \mu \mathbb {E}e_{i+1}(T_{i+1} - t) \biggl(1 - E\biggl[\exp\biggl[-\frac{\lambda}{\mu^2 \mathbb {E}\delta^2} B(T_{i+1} - t)\biggr]\biggr] \biggr), $$

(87)

which can be decomposed into

$$ \psi = I_1 + I_2, $$

(88)

where

(89)

The analysis of I ₁ proceeds exactly as in Appendix A.2 with the modifications of Appendix A.3. I ₂→0 given (86) and the convergence result for I ₁.

We have left the task of proving Lemma A.1. Our approach is a specific implementation of the general approach outlined in Darling and Norris (2008). See Leviyang (2011) for another example of these types of arguments in the context of viral dynamics.

Proof

Define the process M _j (t) by

$$ dM_j = de_j - \bigl(e_j(\gamma h - k_j) + \mu \gamma h e_{j-1} \bigr)\,dt. $$

(90)

M _j (t) is a martingale on the filtration defined by h,e _j ,e _j−1.

We first focus on (85). Using stopping times, we force Δh,Δe _j for j≤i to be small. We then linearize the equations governing Δh,Δe _j and show that the stopping times may be removed while still insuring that Δh,Δe _j are small.

For the subinterval $[T_{i}, T_{i}^{h}]$, define the following stopping times:

(91)

In the definition of T _C , κ is any positive constant less than 1/4. The factor $(\mu \mathbb {E})^{\kappa}$ ensures that $P(T_{C} < T_{i}^{h}) \to 0$ while maintaining sufficient control of e _i+1. Set T=T _A ∧T _B ∧T _C . The equations for Δh and Δe _j on [T _i ,T∧t] with j≤i are

(92)

Set

(93)

Then (92) can be integrated:

(94)

where Σ is a matrix defined through (92) that depends on $h^{\text {det}}$ and $e_{j}^{\text {det}}$. Note that the first integral to the right of the equality directly above is a martingale.

Now, consider $t \wedge T \wedge T_{i}^{h}$. Since, by definition, $T_{i}^{h} - T_{i} = O(|\log(\delta)|)$, we have

$$ \bigg\|\exp\biggl[\int_s^{t \wedge T \wedge T_i^h} ds' \varSigma\bigl(s'\bigr)\biggr]\bigg\|_\infty = O\biggl( \frac{1}{\delta^c}\biggr) $$

(95)

for some constant c that depends only on the parameters of (12). Then using the Doob inequality and noting that the quadratic variation of M(t) is $O(\frac{1}{\mathbb {E}})$, we can arrive at

(96)

where the c _i are constants independent of the SPL. A Chebyshev bound and the limit $(\mu^{2} \mathbb {E})^{2} (\mu \mathbb {E})^{\kappa}/\delta^{c} \to 0$ then give

$$ P\biggl(\sup_{s \in [T_i, T_i^h \wedge T]} \big|v(s)\big| > \frac{1}{\mathbb {E}^{\frac{1}{4}}}\biggr) \to 0. $$

(97)

The above limit shows that the stopping times T _A ,T _B may be removed from T outside of a set with collapsing probability under the SPL. The arguments given in the proof of Lemma 5.2 can be used to show $P(T_{C} < T_{i}^{h}) \to 0$, and so we arrive at

$$ P\biggl(\sup_{s \in [T_i, T_i^h]} \big|v(s)\big| > \frac{1}{\mathbb {E}^{\frac{1}{4}}}\biggr) \to 0. $$

(98)

Now we consider the time interval $[T_{i}^{h}, T_{i+1}]$. Define

(99)

Set $T' = T_{A}' \wedge T_{B}' \wedge T_{C}'$. Note that $T_{C}'$ is defined with respect to 2δ.

On $[T_{i}^{h}, T']$ the arguments of Sect. 5.1 that involve deterministic dynamics give $e_{j}^{\text {det}}= O(\delta)$ for j<i and $\gamma h^{\text {det}}- k_{i} = O(\delta)$. For j<i, we have the bound

$$ d\Delta e_j = \bigl(\gamma h^\text {det}- k_j \bigr) \Delta e_j + O(\delta)\Delta h + dM_j + O\biggl( \frac{1}{(\mathbb {E}^{\frac {1}{8}})^2} + \frac{\mu}{ \mathbb {E}^{\frac{1}{8}}}\biggr). $$

(100)

Applying the same arguments as above and noting that $T' - T_{i}^{h} = O(1/\delta)$ gives for j<i,

$$ P\biggl(\sup_{s \in [T_i^h, T']} \big|\Delta e_j(s)\big| > \frac{1}{\mathbb {E}^{\frac{1}{8}}} \biggr) \to 0. $$

(101)

We need to control Δe _i ,Δh. It is easy to check that the 2×2 matrix that gives the linearized dynamics of Δe _i ,Δh on $[T_{i}^{h}, T']$ is negative definite. Then similar arguments as used above show

$$ P\biggl(\sup_{s \in [T_i^h, T']} \big|\Delta e_i(s)\big| > \frac{1}{\mathbb {E}^{\frac{1}{8}}}\biggr) \to 0 $$

(102)

and similarly for Δh. Equation (102) allows us to remove $T_{A}', T_{B}'$ from the definition of T′. As a result, we have demonstrated (85), except that t is restricted to $[T_{i}, T_{C}']$ instead of [T _i ,T _i+1].

We now turn to e _i+1 and the proof of (86). Clearly, Δe _i+1(t)=0 for $t \in [T_{i}, T_{i}^{h}]$, and we need only consider the interval $[T_{i}^{h}, T_{i+1}]$. From the results above and the analysis of Sect. 5.1, for $t \in [T_{i}^{h}, T_{C}']$, we have

$$ de_{i+1}^\text {det}= e_{i+1}^\text {det}\bigl(\Delta k+ O(\delta)\bigr). $$

(103)

This gives

$$ d \biggl(\frac{e_{i+1}}{e_{i+1}^\text {det}} \biggr) = \frac{e_{i+1}^\text {det}}{e_{i+1}^\text {det}} \biggl(O\biggl(\frac{1}{(\mathbb {E}^{\frac {1}{8}})^2} + \frac{\mu}{ \mathbb {E}^{\frac{1}{8}}} \biggr) \biggr) + \frac{dM_{i+1}}{e_{i+1}^\text {det}}. $$

(104)

Integrating (104) gives for $t \in [T_{i}^{h}, T_{C}']$,

(105)

Taking the second moment, noting that $T_{C}' - T_{i}^{h} = O(1/\delta)$, and using the Doob inequality give

(106)

for some constant c. Since $\mu^{2} \mathbb {E}^{2} \delta^{c} \to \infty$ by the Chebyshev inequality, we have proved (86) except that $t \in [T_{i}, T_{C}']$.

The lemma is proved except that we have considered $[T_{i}, T_{C}']$ rather than [T _i ,T _i+1]. But we have

$$ \frac{e_{i+1}(T_C')}{e_{i+1}^\text {det}(T_C')} \to 1, $$

(107)

and so e _i+1(T _C )/2δ→1. But this implies $P(T_{i+1} < T_{C}') \to 1$, and we are done. □