Modeling in-host Hepatitis C infection

( previous, home, next )

Facts about HCV

Simple model of viral infection in-host

HCV life cycle (Nature Reviews Gastroenterology & Hepatology 6, 403-411 (July 2009))

HCV life cycle
HCV life cycle

Patient data: HCV viral load in plasma following interferon therapy

Viral load for patient 1e Viral load for patient 3d

Data from Neumann et al. (1990).

Patient 1e (left)

[ Data : hide , shown as table , shown as CSV shown as QR ]

# time (days), viral load (log10 per mL)
#
# HCV viral load data following interferon therapy
# Patient 1E from Neumann et al. (1998)
# http://www.jstor.org/stable/2897780
#
0.0000,100
0.0833,92
0.1666,140
0.2916,110
0.4166,98
0.5833,41
0.7916,12
1.0000,7
1.2083,7.5
1.4167,7
2,8
3,2.4
4,3.9
5,2.6
7,1.2
9,0.4
11,0.45
14,0.05

Patient 3d (right)

[ Data : hide , shown as table , shown as CSV shown as QR ]

# time (days), viral load (log10 per mL)
#
# HCV viral load data following interferon therapy
# Patient 3D from Neumann et al. (1998)
# http://www.jstor.org/stable/2897780
#
0,38
0.08333,44
0.29167,45
0.41667,6
0.58333,3.8
0.79167,1.7
1,1
1.2083,0.65
1.4167,0.5
2,0.15
3,0.2
5,0.05
7,0.05
9,0.01

Simple model, chronic infection

Linear compartmental model of viral dynamics
Linear compartmental model of viral dynamics

Let \(I\) represent the number of infected hepatocyte cells at time \(t\) and \(V\) the number of virus, in a single mL. Which can be written in ODE form as \[\begin{align} \frac{dI}{dt} &= \beta(T) V - \delta I \\ \frac{dV}{dt} &= p I - c V \end{align}\]

Having a model is great, but just the start of the story. In order to actually use this model, we have to get estimates for the various parameters we have introduced, \(c\), \(p\), \(\delta\), and \(\beta(T)\). That’s allot harder to do than you might first expect when dealing with a micro-organism inside a human body. In fact, before the discovery of some effective treatments, nobody knew how to do it.

Separation of time-scales and the steady-state approximation

At this point, we will apply a very powerful-but-sometimes-incorrect idea to simplify this system of equations – the separation of time-scales. The idea, in rough terms, is that different parts of large systems often change at different rates. Some parts change quickly, while other parts change slowly. If these speed differences are large enough, then we can sometimes use them to simplify the dynamics of a system. When a component of a system oscillates very quickly, we can try to replace that component with it’s average behavior. And when a component of a system changes very slowly, we can try to replace that component with its steady-state behavior. These are core concepts in the field of applied mathematics known as “asymptotic approximation theory”. While relatively simple to grasp, the mathematics of applying these ideas can be quite complicated, and sometimes leads to serious mistakes if carelessly applied.

Here in our HCV model, though, the idea turns out to be straight-forward and very useful. In long-term, chronic infection, prior to the start of treatment, we expect the number of target cells, infected cells, and virus particles to be roughly constant, day-to-day. Once treatment starts, the number of infected cells and virus will change rapidly, but the number of target cells will not change very quickly because the processes in the human body that control target-cell creation are relatively slow processes. Thus, we can hypothesize that the number of target cells will be at approximately the same steady-state value for the duration of our experiments. At steady-state prior to treatment, the time-derivatives of both the number of infected cells and the amount of virus should vanish, so \[\begin{align} 0 &= \beta(T) V - \delta I, \\ 0 &= p I - c V. \end{align}\] Solving these equations, we find \(\beta(T) = c\delta/p\). Substituting back in, our potentially non-linear system for viral dynamics becomes a linear 2x2 system. \[\begin{align*} \frac{dI}{dt}&= \frac{c\delta}{p} V - \delta I, \\ \frac{dV}{dt}&= pI-cV. \end{align*}\]

Treatment

Estimate parameters using data following initiation of therapy. But… what does interferon therapy do?

As it turns out, interferon therapy blocks the ability of cells to produce virus particles – the more effective it becomes, the fewer virus particles get produced. We can parameterize interferon treatment efficacy, then, as the fraction of viral production blocked, \(\epsilon\) with the natural constraint that \(0\leq \epsilon\leq 1\). Introducing this into our linear system, then, \[\begin{align} \frac{dI}{dt} &= \frac{c\delta}{p} V - \delta I \\ \frac{dV}{dt} &= (1-\epsilon)pI - c V \end{align}\] or in matrix form, \[\begin{gather} \frac{d}{dt}\left(\begin{array}{c}I\\V\end{array}\right) = \left(\begin{array}{cc}-\delta & c\delta/p \\(1-\epsilon)p&-c\end{array}\right) \left(\begin{array}{c}I\\V\end{array}\right) \end{gather}\]

We can solve this system to get \(V(t) = C _ 1 e^{-\lambda _ 1 t} + C _ 2 e^{-\lambda _ 2 t}\), where \(\lambda_1\), \(\lambda_2\) are the eigenvalues of the matrix \(\left(\begin{array}{cc}-\delta & c\delta/p \\(1-\epsilon)p&-c\end{array}\right)\),

\[\lambda_{+,-}=-\frac{c+\delta\pm\sqrt{(c-\delta)^2+4c\delta(1-\epsilon)}}{2}\]

An effective therapy should have \(\epsilon\) near one, making \(1-\epsilon\) small, in which case \[\lambda_{+,-}=-c,\;\;-\delta.\]

Note that our data appears biphasic: two phases of exponential decay (y-axis is on a log scale). Our modeling predicts that the slopes of the decay phases will give rough estimates for the infected cell death rate \(\delta\) and the viral clearance rate \(c\):

Pat1e_fits

[Show code]
Pat1e=loadtxt('NeumannEtAl1998_Pat1e.csv',delimiter=',')

# Start by showing data for viral decay phase
figure(1)
plot(Pat1e[:,0], log10(10**4*Pat1e[:,1]), 'b.', markersize=15)
axis([0, 1.1, 4.5, 6.5])
xlabel('Time (Days)',fontsize=24)
ylabel('log$_{10}$[Viral load (per mL)]',fontsize=24)
legend(['Data'],loc=1,fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
plt.tight_layout()
savefig('FirstPhaseData_1e.png')

# do linear fit for viral clearance
coeff = polyfit(Pat1e[4:8,0], log(10**4*Pat1e[4:8,1]), 1)
t=linspace(0,14,1000)
Phase1=exp(coeff[1])*exp(coeff[0]*t)
plot(Pat1e[4:8,0], log10(10**4*Pat1e[4:8,1]), 'g.', markersize=15)
plot(t,log10(Phase1),'g--',linewidth=2)
legend(['Data','First phase data', 'First phase fit'],loc=1,fontsize=20)
savefig('FirstPhaseData_EstimateC_1e.png')

# Now show second phase data
axis([0, 14.2, 2, 6.5])
savefig('SecondPhaseData_1e.png')

# do linear fit for infected cell death rate
coeff2 = polyfit(Pat1e[8:,0], log(10**4*Pat1e[8:,1]), 1)
t=linspace(0,14,1000)
Phase1=exp(coeff2[1])*exp(coeff2[0]*t)
plot(Pat1e[8:,0], log10(10**4*Pat1e[8:,1]), 'r.', markersize=15)
plot(t,log10(Phase1),'r--',linewidth=2)
legend(['Data','First phase data', 'First phase fit','Second phase data', 'Second phase fit'],loc=1,fontsize=20)
savefig('SecondPhaseData_EstimateDelta_1e.png')


### Plot lambda1, lambda2 as a function of epsilon assuming c=0.39 and delta=2.57 as a function of epsilon
figure(2)
c=-coeff2[0]
delta=-coeff[0]
epsilon=linspace(0,1,100)
lambda1=-(c+delta)/2+sqrt((c+delta)**2-4*c*delta*epsilon)/2
lambda2=-(c+delta)/2-sqrt((c+delta)**2-4*c*delta*epsilon)/2
plot(epsilon,lambda1,'r-',linewidth=2)
#plot(epsilon,-0.39*ones_like(epsilon),'r--',linewidth=2)
plot(epsilon,lambda2,'b-',linewidth=2)
#plot(epsilon,-2.57*ones_like(epsilon),'b--',linewidth=2)
xlabel('IFN efficacy $\epsilon$',fontsize=24)
ylabel('Eigenvalues',fontsize=24)
legend(['$\lambda_1$; $\lambda_1(1) = -\delta$','$\lambda_2$; $\lambda_2(1) = -c$'],loc=2,fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
plt.tight_layout()
savefig('Eigenvalues_1e.png')

For patient 1e the estimates are \(\delta=0.36\) per day and \(c=4.57\) per day. Based on our model we estimate that in patient 1e, infected cells die after an average of \(1/\delta\approx 3\) days.

Model prediction: how many virus per day? We can use our estimates to predict how many virus are produced in an infected individual per day. In our model that quantity is \(pI\). In the presence of therapy, \(\dot{V}=0=pI-cV\), so the quantity we’re looking for is \(=cV\). We prefer \(cV\) because we have an observations of \(V\): the average viral load before therapy. For patient 1e, the average is \(V_0=5.5\times 10^5\) virus per mL. Thus we predict that \(2.6\times 10^9\) virus produced per day per mL. Assuming 15L extracellular fluid (average human) that’s \(\mathbf{3.9\times 10^{10}}\) HCV virions produced per day in an infected individual!

( previous, home, next )