September 22, 2017, Friday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Two-phase approximate solution to an acute viral infection model

From MathBio

Jump to: navigation, search

Smith, A. M., F. R. Adler and A. S. Perelson. (2010) An accurate two phase approximate solution to an acute viral infection model. J. Math. Biol. 60:711-726.[1] This paper reviewed by Zach Vaughan and JaeAnn Dwulet.

Executive Summary

This paper developes a target cell limited influenza model and finds analytical approximations for the two exponential phases of the viral infection. Phase I represents the exponential growth of the infection while phase II is the phase when infection resolution occurs in the form of exponential decay. Modified systems of ordinary differential equations are evaluated for each phase. There is a gap between Phase I and Phase II that Smith et al. do not analyze in this paper due to its non-linearity. It is the subject of their current work.

We reproduce their results using a simbiology model. We also used simbiology for parameter estimation to determine if the program could find the parameter values that would produce equal or better curves for fitting of the data. \beta and p were fixed because the best fit values of these two parameters were unphysiological. The values found through simbiology's parameter fitting could not produce a better fit of the data than the origional values.


The two exponential phases of Influenza A Virus (IAV) infection is analyzed. This is an acute viral infection which means it lasts about 7-10 days. The models are intended for the comparison and analysis of viral properties. This can be used to create more efficient drug therapies and individual drug therapies based on the specific infection. The data collected for the viral load came from the human nasal wash of six patients. [2]


Acute viral infections occur when the subject is initially infected by a virus, the virus levels rise, reach a peak and then decline as the body fights back. We have seen several mathematical models used to study the infection kinetics of various viruses. More often then not, these models are nonlinear and analysis is limited to parameter estimation and numerical simulations. The goal of this article is to find approximate analytically tractable solutions which could provide an easier interpretation of the virus, and investigate how varying the parameters effects the infection kinetics.

Mathematical Model

Basic Model

The basic model found below, describes a viral infection with three different states. T represents target cells, I represents infected cells, and V represents the free virus cells.

(1) {\begin{cases}{\frac  {\partial T}{\partial t}}=s-dT-\beta TV\\{\frac  {\partial I}{\partial t}}=\beta TV-\delta I\\{\frac  {\partial V}{\partial t}}=pI-cV\end{cases}}

In system (1), the uninfected cells are supplied at a constant rate s, die at rate d, and become infected at rate \beta V. The Infected cells I are either self degraded through apoptosis, or removed by immune cells, at rate \delta . Last, the virions, V are produced at rate p, and cleared at rate c.

This type of model can be very helpful under certain conditions, but a more accurate model may be desired. This leads us to the two phase limited target cell model.

Two Phase Exponential Model

First, the basic model is generalized to include the amount of time it takes for an infected cell to produce new virions. These cells are represented by I_{1}, and are cells that the virus has attacked but have yet to produce new virions because this process takes . For the new model, we assume that the non-producing infected cells I_{1} enter the I_{2} state (producing) at rate k, and are not subject to cell death until they are actively producing. Once active, the cells are removed at rate \delta , but the natural apoptosis cell death rate d is now assumed to be zero, because we are looking on such a short time scale. Last, the target cell production rate s is set to zero under the assumption that very few epithelial cells are actually regenerated by the time the virus is cleared. These assumptions lead us to the new model which we will be focusing on for the rest of the summary.

(2) {\begin{cases}{\frac  {\partial T}{\partial t}}=-\beta TV\\{\frac  {\partial I_{1}}{\partial t}}=\beta TV-kI_{1}\\{\frac  {\partial I_{2}}{\partial t}}=kI_{1}-\delta I_{2}\\{\frac  {\partial V}{\partial t}}=pI_{2}-cV\end{cases}}

A numerical solution to system 2 of each of the different states is shown below in figure 1. Figure 2 shows just the virion cell population on a log scale with the actual viral data of patient 4.

Figure 1 Fig1ajz2.jpg
Figure 2 Fig2bjz.jpg

All of the parameter values used to produce the numerical solution were estimated by Baccam et al. (2006) for each of the individual 6 patients. These values are found in the table below.

Table 1 Table1jz.jpg

Phase I

One of the disadvantages of this model is that it is non-linear and so only a numerical solution can be found. In order to combat this, the authors generalized the model into two separate phases, and could therefore make valid assumptions to transform the system from being non-linear to linear. As one can see in Fig. 1, the target cell population starts at the initial concentration and then rapidly declines to zero. Knowing this, the authors create two phases, with the initial phase having a constant T population T_{0}, and phase two having a T population of 0. Making T constant, reduces the system to only linear terms, and so an analytical solution can be found. The new model for phase 1 now becomes

(3) {\begin{cases}T_{0}=4\times 10^{8}\\{\frac  {\partial I_{1}}{\partial t}}=\beta T_{0}V-kI_{1}\\{\frac  {\partial I_{2}}{\partial t}}=kI_{1}-\delta I_{2}\\{\frac  {\partial V}{\partial t}}=pI_{2}-cV\end{cases}}

Which produces the following solution

(4) {\begin{cases}T(t)=T_{0}-I_{1}(t)-I_{2}(t)\\I_{1}(t)=\kappa _{1}e^{{\lambda _{1}t}}+\kappa _{2}+e^{{\lambda _{2}t}}+\kappa _{3}e^{{\lambda _{3}t}}\\I_{2}(t)=\kappa _{1}\xi _{{12}}e^{{\lambda _{1}t}}+\kappa _{2}\xi _{{22}}e^{{\lambda _{2}t}}+\kappa _{3}\xi _{{32}}e^{{\lambda _{3}t}}\\V(t)=\kappa _{1}\xi _{{13}}e^{{\lambda _{1}t}}+\kappa _{2}\xi _{{23}}e^{{\lambda _{2}t}}+\kappa _{3}\xi _{{33}}e^{{\lambda _{3}t}}\end{cases}}

with the eigenvalues defined as

\lambda _{1}={\sqrt[ {3}]{-{\frac  {q}{2}}+{\sqrt  {{\frac  {q^{2}}{4}}+{\frac  {u^{3}}{27}}}}}}+{\sqrt[ {3}]{-{\frac  {q}{2}}-{\sqrt  {{\frac  {q^{2}}{4}}+{\frac  {u^{3}}{27}}}}}}-{\frac  {B}{3}} \lambda _{{2,3}}=(-{\frac  {1}{2}}\pm i{\frac  {{\sqrt  {3}}}{2}}){\sqrt[ {3}]{-{\frac  {q}{2}}+{\sqrt  {-{\frac  {q^{2}}{4}}+{\frac  {u^{3}}{27}}}}}}+(-{\frac  {1}{2}}\mp i{\frac  {{\sqrt  {3}}}{2}}){\sqrt[ {3}]{-{\frac  {q}{2}}-{\sqrt  {{\frac  {q^{2}}{4}}+{\frac  {u^{3}}{27}}}}}}-{\frac  {B}{3}}

The following table shows the actual calculated eigenvalues for each individual patient

Table 2 Table2jz.jpg

Phase II

Phase II is characterized by a slowed growth rate of the free virus cells which leads to a peak in both the free virus and infected cells. The target cells sharply decline and the virus then begins to exponentially decay.

Due to the behavior, we let T be approximately zero and the following conclusions are made: Z\gg T which implies that I_{1}\approx Z-I_{2}

This reduces system (2) to the following:

(5) {\begin{cases}{\frac  {\partial I_{1}}{\partial t}}=-kI_{1}\\{\frac  {\partial I_{2}}{\partial t}}=kI_{1}-\delta I_{2}\\{\frac  {\partial Z}{\partial t}}=-\delta I_{2}\\{\frac  {\partial V}{\partial t}}=pI_{2}-cV\end{cases}}

Phase II begins at time t_{2}, the time at which T(t) is sufficiently small or when it is 10% of its initial value and can be approximated to be zero. The gap between phase I and phase II cannot be approximated by this model and therefore there is not an explicit solution for t_{2} because the values for the parameters at the end of phase I are not the values at the beginning of phase II. The values of t_{2} for each of the patients can be seen in table 2, which were found numerically by satisfying

T(t_{2})=0.1T_{0}=4\times 10^{7}

The solution to this system is a summation of exponentials involving the parameters k,c and \delta . The solution takes on a different form when k=c,k=\delta and/or c=\delta which can be solved using L'Hopital's rule

Alternate Phase II

The phase II solution is only valid when the number of target cells is small at time t_{2} and continues to approach zero as t\rightarrow \infty . As seen in Patient 5, the phase II solution does not fit the approximate solution. Patient 5's data shows a much larger number for the final target cell count. This type of behavior is due to a small R_{0} value, the reproductive ratio. When R_{0}\lesssim 7, the final cell count for target cells cannot be approximated at zero (Figure 5), instead the number of target cells approaches a steady state making the phase II solution inaccurate. This alternate approximation occurs at a later time, {\tilde  {t}}_{2} and does not include the viral peak. This can occur when the virus or the infected cells are cleared by the immune cells fast enough to prevent further infection. In this case T_{{\infty }}, or the final target cell count can be approximated using the classic epidemic model by Kermack and Mekendrick model by satisfying ln(T_{{\infty }}/T_{0})=R_{0}(T_{{\infty }}/T_{0}-1).

Because the target cell count approaches a non zero constant the same approach can be used to solve for this alternate phase II as with phase I when the total number of cells, Z(t) remained constant. Therefore we can use the system of equations from the phase I solution with an alternative value of Z(t), and T_{0} can be replaced by T_{{\infty }}

Figure 3

Simbiology Model

We used simbiology to find numerical solutions of system (2) using the following set up.


We then used a nonlinear least squares fit to see if we could find parameters that better fit patient 4's data. Initially we varied \beta ,k,\delta ,p, and c. While this produced a better fit, the parameters \beta and p were negative, which is an unphysiological situation. We deduced that this could be because both \beta and p are several magnitudes smaller then the other parameters. We then just varied k,p and \delta , and produced the following numerical solution of the virus population.

Figure 4

As one can see, this fit is not as accurate as the original. We calculated the new parameters as c=1.2779,\delta =1.1891, and k=1.1873. By limiting ourselves to varying only a select number of parameters, we subject the solution to not being as accurate as possible. This also shows that \beta and p while small, still have a large effect on the interacting system.


  • Phase I

The leading eigenvalue, \lambda _{1}, determines the exponential viral growth of phase I. \lambda _{1} is very sensitive to changes in k. The eigenvalue, \lambda _{1} converges to the leading eigenvalue r_{0} as k\rightarrow \infty , which can be seen in Figure 5.

Lambda 1.jpg
Figure 5

Using the basic model, Nowak et al. [3] found that when T is constant, viral growth is described by e^{{r_{0}t}} where r_{0} is the leading eigenvalue that solves the equation r_{0}^{2}+(c+\delta )r_{0}-c\delta (R_{0}-1)=0. R_{0} is the basic reproductive ratio, R_{0}=\beta ps/(s\delta d). The eigenvalue, \lambda _{1} converges to the leading eigenvalue r_{0} as k\rightarrow \infty , which can be seen in (Figure 6). In physical terms, when rate k is extremely fast, I_{1} immediately enter the I_{2} state. The equations then look like the basic model which has only one infected class and therefore it makes sense that the leading eigenvalue \lambda _{1} approaches the leading eigenvalue, r_{0} of the model with only one infected class. For physically acceptable values of k,\lambda _{1} is significantly lower than r_{0}.

The other eigenvalues \lambda _{2} and \lambda _{3} determine the initial dip in virus levels that occurs at the beginning of phase I. After this occurs the positive \lambda _{1} drives viral growth.

  • Phase II

The rate of decay depends on the smallest of the parameters, k,c, and \delta and the peak of the virus is determined by all three of these parameters. Smith et al. believes that \delta , the infected cell death rate, often determines the rate of decay.

  • Alternate Phase II

The value for the target cells left late in the reaction can be determined by T_{{\infty }}\approx T_{0}e^{{-R_{0}}}. This indicated a small value for R_{0}, the reproductive ratio. When R_{0}\lesssim 7 the Phase II model fails and the alternative phase II can be approximated with the equations for phase I using the steady state target cell count.


  • The Phase I equation for the leading eigenvalue can be used to evaluate the differences between hosts, and individual virus strains. The eclipse phase length can lead to very different viral growth rates. Fast viral growth or large \lambda _{1} leads to rapid target cell depletion and an earlier peak in viral cell count meaning that the first phase of the virus is shorter. On the other hand slow viral growth leads to a longer Phase I. These conclusions can be used to determine the most effective time at which a drug such as Tamiflu should be administered. This can lead to individualized drug therapies especially for a virus such as H1N1. In conclusion, Phase I would be the most effective time to treat a patient with a drug.
  • When k\gg \delta and c\gg \delta , the slope of the a linear fit to Phase I can be used to estimate the value of \lambda _{1} and this can be used to compare host and/or viral strains. This can be seen in the following plot, where a log linear regression was run on patient 4's data.
Lin fit.jpg
FIgure 6
  • The value of the smallest parameter k,c, or \delta can be estimated from the slope of a linear fit of phase II. This approach is only useful if one knows or has a good idea of which parameter is the smallest.
  • The gap between Phase I and Phase II is still under investigation due to its non-linear characteristics.

Biological Applications and Questions

  • Biologists can use this model by determining the parameters to compare host and/or viral strains. The ability to find a linear fit for phase I gives a simple way for biologists to compare host and viral strains by estimating \lambda _{1} from known parameters. With this information can biologists create individualized drugs and drug treatments?
  • Another way that biologists can use this without knowing the parameters specifically is to determine the extent of damage to epithelial cells and virus transmissible. These quantities can be evaluated in terms of parameters by the integrals of I_{2}(t) and V(t).
  • Biologists can also compare virus strains without knowing the parameters specifically by using estimated parameters that are reasonable chosen. They can be chosen by comparing data of similar viruses. For example genetically engineered influenza strains may have similar parameter values to the IVA strain analyzed in this paper. As long as reasonable values for the other parameters are chosen the phase behaviors of each strain can be tested.


These are the MATLAB codes we used to regenerate their plots

  • Plot of the leading eigenvalue approaching r_{0} as k\rightarrow \infty
clear all 
close all

beta = 5.3*10^-6;
p = 0.13;
K = 0:1:400;
c = 3.8;
delta = 3.8;

T0 = 4*10^8;
V0 = 4.9;
R0 = beta*p*T0/(c*delta);

%rind ro 


B = K;
C = K;
D = K;
lambda1 = K;

for ii = 1 : length(K);
B(ii) = K(ii)+delta+c;
C(ii) = K(ii)*delta+K(ii)*c+c*delta;
D(ii) = -K(ii)*c*delta*(R0-1);

u = C(ii) - B(ii)^2/3;
q = D(ii)+ (2*B(ii)^3-9*B(ii)*C(ii))/27;

lambda1(ii) = (-q/2+sqrt(q^2/4+u^3/27))^(1/3)+(-q/2-sqrt(q^2/4+u^3/27))^(1/3)-B(ii)/3;
%lambda2 = (-.5+1i*(sqrt(3)/2))*(-q/2+sqrt(q^2/4+u^3/27))^(1/3) + (-.5-1i*(sqrt(3)/2))*(-q/2-sqrt(q^2/4+u^3/27))^(1/3)-B/3;
%lambda3 = (-.5-1i*(sqrt(3)/2))*(-q/2+sqrt(q^2/4+u^3/27))^(1/3) + (-.5+1i*(sqrt(3)/2))*(-q/2-sqrt(q^2/4+u^3/27))^(1/3)-B/3;


plot(K, lambda1)
hold on
plot(K, r01)
  • This generates the plot with linear fit of each phase
 clear all 
close all

t_data = [1,2,3,4,5,6,7];
V_data = [10^3.4, 10^5.5, 10^6.5, 10^5.6, 10^3.5, 10^4, 10^0];

hold on;

t1 = t_data(1:3);
V1 = log(V_data(1:3));

t2 = t_data(3:7);
V2 = log(V_data(3:7));

p1 = polyfit(t1, V1,1);
p2 = polyfit(t2, V2,1);

y1 = t1*p1(1)+p1(2);
y2 = t2*p2(1)+p2(2);

plot(t1, y1);
hold on
plot(t2, y2);


  1. Smith, A. M., F. R. Adler and A. S. Perelson. (2010) An accurate two phase approximate solution to an acute viral infection model. J. Math. Biol. 60:711-726.
  2. Baccam P.,Beauchemin C., Macken C.A., Hayden F.G., Perelson A.S. (2006) Kinetics of influenza A virus infection in humans. J Virol 80(15):7590-7599.
  3. Nowak M.A., Lloyd A.L. Vasquez G.M., Wiltrout T.A., Wahl L.M. Bischofberger N. Williams J., Kinter A., Fauci A.S. Hirsch V.M. Lifson J.D. (1997) Viral dynamics of primary viremia and antiretroviral therapy in simian immunodeficiency virus infection. J Virol 71(10):7518-7525.