13 - A within-host HIV model

Case study: HIV

Human immunodeficiency virus (HIV) infects immune cells called ‘CD4\(^+\) T-cells’, which form an important part of the human immune system. HIV enters these T-cells, begins to replicate and then releases new virus particles in to the bloodstream. Immediately after infection, the virus grows rapidly and produces flu-like symptoms in the patient. After a few months, these symptoms disappear and the virus concentration reduces to a lower, but steady, level. This ‘asymptomatic period’ can last for years, with the virus density staying roughly constant, and the concentration of T-cells very slowly dropping. Eventually, the T-cell density becomes so low that the patient’s immune system is no longer effective, a condition called aquired immunodeficiency syndrome (AIDS). The patient is then at risk from life-threatening opportunistic infections. In the mid-1990s, much work was being done to explore the dynamics of HIV and to produce potential treatments, including the development of mathematical models led by Perelson and Nelson.

A (pre-treatment) model

As with the epidemic models, let’s start by establishing our variables and drawing a schematic of what is happening in this system. The variables that we want to keep track of (before treatment, at least), are:

• Healthy T-cells (\(T\))

• Infected T-cells (\(T^*\))

• Virus particles (\(V\))

We will only keep track of the concentration of virus particles in the bloodstream, not the ‘intracellular’ concentrations. Our schematic should look like this:

Healthy T-cells are produced at some constant rate \(s\) by the body, but also decay at some rate \(d\). These healthy cells become infected through contact with virus, in a mass-action manner, with coefficient \(k\). These infected T-cells then decay at some rate \(\delta\). New virus particles are produced when infected T-cells die, with \(N\) giving the average number of virus particles produced by each cell. The virus then also decays at some rate \(c\). The system can therefore be expressed with the following set of ODEs:

(42)\[\begin{align} \frac{dT}{dt} &= s-dT-kVT\label{hiv1}\\ \frac{dT^*}{dt} &= kVT-\delta T^*\label{hiv2}\\ \frac{dV}{dt} &= N\delta T^*-cV\label{hiv3} \end{align}\]

Looking at the equations, we can see that there is an HIV-free equilibrium for \((T,T^*,V)=(s/d,0,0)\). We have another example here of a 3-dimensional system. As we saw in the free-living parasite model, assessing stability here can be a bit more complicated, often reying on using the Routh-Hurwitz criteria. As it turns out, for the HIV-free equilibrium things simplify fairly nicely to the point that we can actually directly calculate one of the eigenvalues, and then just use the trace-determinant conditions (themselves the Routh-Hurwitz criteria for a 2-dimensional system) to determine the sign of the other two.

The general Jacobian is,

\[\begin{align*} J=&\left( \begin{array}{ccc} -d-kV & 0 & -kT\\ kV & -\delta & kT\\ 0 & N\delta & -c \end{array} \right). \end{align*}\]

Subsituting in the equilibrium values this becomes,

\[\begin{align*} J=&\left( \begin{array}{ccc} -d & 0 & -ks/d\\ 0 & -\delta & ks/d\\ 0 & N\delta & -c \end{array} \right). \end{align*}\]

To find the eigenvalues we want to write out the characteristic equation - the determinant of the matrix of \( (J-\lambda I_3)\), where \(I_3\) is the identity matrix. The characteristic equation is,

\[\begin{equation*} (-d-\lambda)\left[(-\delta-\lambda)(-c-\lambda)-\frac{ksN\delta}{d}\right]=0\\ \end{equation*}\]

This is already partly factorised, making our life much easier. The first bracket tells us that one eigenvalue is \( \lambda=-d\lt0\). To find the remaining two eigenvalues we look at what is inside the square bracket, wich can be re-written as,

\[\begin{equation*} \lambda^2+(c+\delta)\lambda+\delta c-\frac{ksN\delta}{d}=0\\ \end{equation*}\]

We could use the quadratic formula to solve for the eigenvalues explicitly, but we can also note that this is just the characteristic equation for a two dimensional system, with \( -(c+\delta)\) being the trace and \( \delta c-\frac{ksN\delta}{d}\) being the determinant. We can then see that we definitely have a negative trace, and that the determinant is positive - meaning the equilibrium is stable - if \( c\gt ksN/d\). In other words, the decay rate of the virus must be high, and the infection rate and production rate must be low.

Of course, when a patient is to undergo treatment, they will not be HIV-free! We are therefore more interested in the pre-treatment infected steady-state. Setting the equations for \(dT^*/dt\) and \(dV/dt\) to zero, yields,

(43)\[\begin{equation} N\delta T^*_0=cV_0 \text{ and } kV_0T_0=\delta T^*_0 \implies T_0=\frac{c}{Nk} \end{equation}\]

Solving dT/dt=0, then gives,

(44)\[\begin{equation} V_0=\frac{sN}{c}-\frac{d}{k} \end{equation}\]

and, by substituting back in,

(45)\[\begin{equation} T^*_0=\frac{kV_0T_0}{\delta} \end{equation}\]

Treatment

Assume a patient arrives with their virus and T-cell concentrations at the pre-treatment steady state. We then wish to put them on a course of treatment. One of the most effective early treatments for HIV was using ‘protease inhibitors’. These drugs did not prevent infection, but meant that any new virus produced inside infected T-cells were non-infectious. The system therefore becomes:

(46)\[\begin{align} &\frac{dT}{dt} = s-dT-kV_IT\label{hiv11}\\ &\frac{dT^*}{dt} = kV_IT-\delta T^*\label{hiv12}\\ &\frac{dV_I}{dt} = -cV_I\label{hiv13}\\ &\frac{dV_{NI}}{dt} = N\delta T^*-cV_{NI}\label{hiv14} \end{align}\]

In the early stages of treatment, it is reasonable for us to assume that the number of healthy T-cells stays roughly constant at it’s pre-treatment level \(T(t)=T_0\), so dT/dt can be ignored. We also know that before treatment, all virus particles were of the infective type, such that \(V_I(0)=V_0\) and \(V_{NI}(0)=0\). Given this, it is easy to show that,

(47)\[\begin{equation} V_I(t)=V_0e^{-ct}\label{vt} \end{equation}\]

Given that \(latex T(t)=T_0\), we now find that \( dT^*/dt\) is also now linear,

(48)\[\begin{equation} \frac{dT^*}{dt} +\delta T^*= kV_IT_0 \end{equation}\]

It may not initally be obvious that this is linear, since the right-hand side had \( V_I\times T_0\). However, \( T_0\) is a constant and we already have an expression for \( V_I=V_0e^{-ct}\). This can therefore be solved using an integrating factor. Let’s go through the working for this.

Taking an integrating factor of \( e^{\delta t}\) we can write this as,

(49)\[\begin{equation} \frac{d}{dt}\left[T^*e^{\delta t}\right] = kV_0T_0e^{-ct}e^{\delta t} \end{equation}\]

Integrating both sides with respect to \( t\) we then get,

(50)\[\begin{equation} T^*e^{\delta t} = \frac{kV_0T_0}{\delta-c}e^{-ct}e^{\delta t}+C \end{equation}\]

where \( C\) is a constant of integration. After rearranging we reach our initial result,

(51)\[\begin{equation} T^* = \frac{kV_0T_0}{\delta-c}e^{-ct}+Ce^{-\delta t}. \end{equation}\]

Earlier we found the pre-treatment equilibrium for \( T^*_0=kV_0T_0/\delta\), so we can say that when \( t=0\) we have,

(52)\[\begin{equation} \frac{kV_0T_0}{\delta} = \frac{kV_0T_0}{\delta-c}+C. \end{equation}\]

We then rearrange to find the constant,

(53)\[\begin{equation} C=\frac{-ckV_0T_0}{\delta(\delta-c)} \end{equation}\]

and then substitute back in to get our final solution,

(54)\[\begin{equation} T^*(t)=\frac{kT_0V_0}{\delta(\delta-c)}\left(\delta e^{-ct}-ce^{-\delta t}\right). \end{equation}\]

Phew! Now it’s your turn….

Have a go

Using similar methods show that,

(55)\[\begin{equation} V_{NI}(t)=\frac{cV_0}{\delta-c}\left(\frac{c}{\delta-c}\left(e^{-\delta t}-e^{-ct}\right)+\delta te^{-ct}\right) \end{equation}\]

Click for solution

Firstly we can write this equation out in the correct form for using an integrating factor,

(56)\[\begin{equation} \frac{dV_{NI}}{dt}+cV_{NI} = N\delta T^*=\frac{NkT_0V_0}{\delta-c}\left(\delta e^{-ct}-ce^{-\delta t}\right). \end{equation}\]

The integrating factor is therefore \( e^{ct}\), which leads us to,

(57)\[\begin{equation} \int\frac{d}{dt}\left(V_{NI}e^{ct}\right)dt = \int\frac{NkT_0V_0}{\delta-c}\left(\delta-ce^{-\delta t+ct}\right)dt. \end{equation}\]

Computing these two integrals we find,

(58)\[\begin{align} V_{NI}e^{ct} &= \frac{NkT_0V_0}{\delta-c}\left(\delta t-\frac{c}{c-\delta}e^{-\delta t+ct}\right)+A\\ \implies V_{NI}&=\frac{NkT_0V_0}{\delta-c}\left(\delta te^{-ct}-\frac{c}{c-\delta}e^{-\delta t}\right)+Ae^{-ct} \end{align}\]

where \( A\) is the constant of integration.

At the start of treatment we’d have \( V_{NI}(0)=0\). Substituting this in tells us,

(59)\[\begin{equation} A=\frac{NkT_0V_0}{\delta-c}\left(\frac{c}{c-\delta}\right), \end{equation}\]

and so we have (after a little rearranging),

(60)\[\begin{equation} V_{NI}=\frac{NkT_0V_0}{\delta-c}\left(\delta te^{-ct}+\frac{c}{\delta-c}(e^{-\delta t}-e^{-ct})\right). \end{equation}\]

Finally we use our assumption that \( T_0=c/(Nk)\) to reach the final solution.

This can be added to the soluion for \(V_I(t)\) to find the overall density of virus particles in the bloodstream at any time point after treatment has begun. After estimating the values of each of the model parameters, Perelson and his colleagues plotted the model’s predictions against data of individual patient’s virus concentrations. The figure below shows that the model provides a remarkably good fit.

Figure: Virus concentrations in the bloodstream of three patients compared to the model predictions. Taken from Perelson & Nelson (1999).

3 key points

1. A patient will only be disease-free if the virus decay rate is high, and its infection rate and burst size is low.

2. Under a simplifying assumption the treatment model becomes linear, so we can solve the model explicitly.

3. The treatment model predicts exponential decay of the virus, and matches well to real data.