
MBW:Twophase approximate solution to an acute viral infection modelFrom MathBioSmith, 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:711726.^{[1]} This paper reviewed by Zach Vaughan and JaeAnn Dwulet. ContentsExecutive SummaryThis 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 nonlinearity. 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. and 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. ContextThe two exponential phases of Influenza A Virus (IAV) infection is analyzed. This is an acute viral infection which means it lasts about 710 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]} BackgroundAcute 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 ModelBasic ModelThe basic model found below, describes a viral infection with three different states. represents target cells, represents infected cells, and represents the free virus cells. In system (1), the uninfected cells are supplied at a constant rate , die at rate , and become infected at rate . The Infected cells are either self degraded through apoptosis, or removed by immune cells, at rate . Last, the virions, are produced at rate , and cleared at rate . 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 ModelFirst, 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 , 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 nonproducing infected cells enter the 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 , 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.
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. 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. Phase IOne of the disadvantages of this model is that it is nonlinear 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 nonlinear 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 , 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 Which produces the following solution with the eigenvalues defined as
Phase IIPhase 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: which implies that This reduces system (2) to the following: Phase II begins at time , the time at which 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 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 for each of the patients can be seen in table 2, which were found numerically by satisfying The solution to this system is a summation of exponentials involving the parameters and . The solution takes on a different form when and/or which can be solved using L'Hopital's rule Alternate Phase IIThe phase II solution is only valid when the number of target cells is small at time and continues to approach zero as . 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 value, the reproductive ratio. When , 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, 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 , or the final target cell count can be approximated using the classic epidemic model by Kermack and Mekendrick model by satisfying . 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, remained constant. Therefore we can use the system of equations from the phase I solution with an alternative value of , and can be replaced by
Simbiology ModelWe 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 and . While this produced a better fit, the parameters and were negative, which is an unphysiological situation. We deduced that this could be because both and are several magnitudes smaller then the other parameters. We then just varied and , and produced the following numerical solution of the virus population. As one can see, this fit is not as accurate as the original. We calculated the new parameters as and . 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 and while small, still have a large effect on the interacting system. Results
The leading eigenvalue, , determines the exponential viral growth of phase I. is very sensitive to changes in . The eigenvalue, converges to the leading eigenvalue as , which can be seen in Figure 5. Using the basic model, Nowak et al. ^{[3]} found that when T is constant, viral growth is described by where is the leading eigenvalue that solves the equation . is the basic reproductive ratio, . The eigenvalue, converges to the leading eigenvalue as , which can be seen in (Figure 6). In physical terms, when rate is extremely fast, immediately enter the state. The equations then look like the basic model which has only one infected class and therefore it makes sense that the leading eigenvalue approaches the leading eigenvalue, of the model with only one infected class. For physically acceptable values of is significantly lower than . The other eigenvalues and determine the initial dip in virus levels that occurs at the beginning of phase I. After this occurs the positive drives viral growth.
The rate of decay depends on the smallest of the parameters, and and the peak of the virus is determined by all three of these parameters. Smith et al. believes that , the infected cell death rate, often determines the rate of decay.
The value for the target cells left late in the reaction can be determined by . This indicated a small value for , the reproductive ratio. When 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. Discussion
Biological Applications and Questions
AppendixThese are the MATLAB codes we used to regenerate their plots
clear all close all clc 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 aa=1; bb=c+delta; cc=c*delta*(R01); r01=(bb+sqrt(bb^24*aa*cc))/(2*aa) r02=(bbsqrt(bb^24*aa*cc))/(2*aa) 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*(R01); u = C(ii)  B(ii)^2/3; q = D(ii)+ (2*B(ii)^39*B(ii)*C(ii))/27; lambda1(ii) = (q/2+sqrt(q^2/4+u^3/27))^(1/3)+(q/2sqrt(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) + (.51i*(sqrt(3)/2))*(q/2sqrt(q^2/4+u^3/27))^(1/3)B/3; %lambda3 = (.51i*(sqrt(3)/2))*(q/2+sqrt(q^2/4+u^3/27))^(1/3) + (.5+1i*(sqrt(3)/2))*(q/2sqrt(q^2/4+u^3/27))^(1/3)B/3; end plot(K, lambda1) hold on plot(K, r01)
clear all close all clc 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]; plot(t_data,log(V_data),'o'); 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); References
