THE DYNAMICS OF A DELAY MODEL OF HEPATITIS B VIRUS INFECTION WITH LOGISTIC HEPATOCYTE GROWTH

Chronic HBV affects 350 million people and can lead to death through cirrhosis-induced liver failure or hepatocellular carcinoma. We analyze the dynamics of a model considering logistic hepatocyte growth and a standard incidence function governing viral infection. This model also considers an explicit time delay in virus production. With this model formulation all model parameters can be estimated from biological data; we also simulate a course of lamivudine therapy and find that the model gives good agreement with clinical data. Previous models considering constant hepatocyte growth have permitted only two dynamical possibilities: convergence to a virus free or a chronic steady state. Our model admits a third possibility of sustained oscillations. We show that when the basic reproductive number is greater than 1 there exists a biologically meaningful chronic steady state, and the stability of this steady state is dependent upon both the rate of hepatocyte regeneration and the virulence of the disease. When the chronic steady state is unstable, simulations show the existence of an attracting periodic orbit. Minimum hepatocyte populations are very small in the periodic orbit, and such a state likely represents acute liver failure. Therefore, the often sudden onset of liver failure in chronic HBV patients can be explained as a switch in stability caused by the gradual evolution of parameters representing the disease state.


284
S. EIKENBERRY, S. HEWS, J. D. NAGY AND Y. KUANG who harbor an active HBV infection suffer chronic hepatitis [34].Up to 0.06% of HBV-infected persons are likely to die from complications associated with the disease within the year [21].But mortality is not the only way HBV impacts the human population.All who suffer HBV infection experience significant morbidity, ranging from weeks to months of nausea, fatigue, jaundice, and joint pain associated with acute disease to liver cirrhosis or hepatocellular carcinoma characteristic of late-stage chronic infection.
Primarily because of HBV's significance as a global public health threat, the virus and its associated disease have attracted considerable attention from mathematical and theoretical biologists [10,27,29,30,7,6].The models typically used to study HBV dynamics within the host tend to focus on healthy cells, free virus, and infected cells.Most of these models formulate viral infection as a mass action process between free virions and healthy cells.Some further assume a constant influx of healthy cells from an outside source.For example, the following phenomenological model has been widely used in the study of HBV and other viral diseases, including HIV: ( The number or mass of healthy cells is represented by x, that of infected cells by y, and that of the free virion load by v. Healthy hepatocytes enter the liver from an outside source at rate r, and hepatocytes die at rate d.
The mass action term used to model infection of healthy hepatocytes by free virions is biologically problematic in several ways.Gourley et al. [10] showed that this assumption causes the basic reproductive number, R 0 , to depend upon the homeostatic liver size, r/d.Thus, this model makes the dubious prediction that individuals with smaller livers should be less susceptible to HBV infection.Furthermore, the mass action constant, β, has no clear biological meaning.
Gourley et al. [10] extended model (1.1) by assuming HBV infection of healthy cells is governed not by mass action but by a standard incidence function.In particular, the mass action term in (1.1) is replaced by where T is the total number of liver cells, infected or not.Under their formulation R 0 loses its dependence on liver size.Also, β in (1.2) has a clear biological meaning.
It can be interpreted either as the maximum rate at which virions infect healthy cells, or as the probability that a single virion infects any cell in a healthy liver in a unit interval of time.
A shortcoming of both model (1.1) and the Gourley et al. formulation is the assumption that healthy hepatocytes are replenished at a constant rate.Mathematically, this assumption limits the allowed asymptotic behaviors to only two generic cases [10,27,29,30]: if R 0 < 1 then the homeostatic set-point, (r/d, 0, ..., 0), is globally asymptotically stable, but if R 0 > 1 then a persistent chronic steady state is approached in all cases.
On biological grounds, this assumption of constant hepatocyte influx is clearly inappropriate.It is well established that liver regeneration following injury is achieved by widespread hepatocyte proliferation [26,19,31], not "reseeding" from an outside source.Also, as we argue below, complete liver regeneration following partial hepatectomy in animal models implies that hepatocyte proliferation depends on liver mass.Therefore, we propose an alternative model of HBV dynamics in which hepatocyte proliferation rate is described by a logistic function.We show that this assumption can produce rich, biologically significant dynamics when R 0 > 1.In particular, the asymptotic behavior can switch from a stable chronic equilibrium to an attracting periodic orbit, with an unstable chronic equilibrium, and that this switch is governed by both the virulence of the infection and the rate at which hepatocytes regenerate.We argue that the gradual evolution of parameters describing liver regeneration and HBV virulence over the course of chronic infection can lead to a sudden change in disease dynamics as the stability switch is crossed.Since this switch can occur for realistic parameter values, the model predicts that the onset of acute liver failure can be predicted by the sudden onset of oscillations in viral load.
A further advantage of this model is that all parameters can be estimated directly from empirical biological data.To validate our approach, we model a course of lamivudine therapy and find that empirically determined parameter values yield dynamics in close agreement with clinical data.
2. Full model.The liver has a remarkable ability to regenerate following injury.This regeneration is driven by normal hepatocytes, which unlike the terminally differentiated cells dominating most other tissues, retain an inherent "steminess."That is, mature hepatocytes can still divide and replenish others lost in a fully developed liver.Furthermore, this regeneration occurs at an extremely fast rate.For example, rat livers completely regenerate within a week following experimental hepatectomy of two-thirds of the organ (2/3 PHx) [26].Regeneration is accomplished by several cycles of mitosis involving a large fraction of mature hepatocytes, although progressively fewer hepatocytes participate in each subsequent cycle [26,31].Precisely how these regenerative cycles are initiated and controlled is still under debate, but hemodynamics appears to play a definitive role.A 2/3 PHx exposes hepatocytes to a threefold increase in portal vein blood flow, which appears to lower hepatocyte apoptosis rates [25] and perhaps increases their exposure to growth factors [26].
The results of Lambotte et al. [19] indicate that early changes in hemodynamics following both 2/3 PHx or a "temporary" partial hepatectomy "prime" hepatocytes for proliferation.However, regenerative proliferation occurs only when the actual liver mass is reduced.Following a 1/3 PHx, the regenerative response can increase following further resection, at least 10 hours later.Also, regeneration can decrease in response to a liver graft (increase in liver mass) up to 18 hours after the initial resection.These results strongly indicate that hepatocyte proliferative activity is controlled by the actual mass of the liver and that hemodynamics play a role in preparing hepatocytes for proliferation.Therefore, we conclude that the growth signal for hepatocytes is directly proportional to the total liver cell population and that signal strength is approximately uniform throughout the entire liver.We choose, as a reasonable approximation of the cellular response to such a signal, a logistic term to represent hepatocyte proliferation rate.
Like Gourley et al. [10] we employ a standard incidence function to describe virion infection of healthy cells and explicitly consider a delay between viral infection and production.These considerations produce the following model: where To simplify analysis we define and rewrite equation (2.1) as The number of healthy hepatocytes is x(t), p(t) is the number of healthy cells that have been infected but are not yet producing virions (latent infection), y(t) is the number of infected cells actively producing virions (productively infected cells), and v(t) is the number of free virions.
2.1.Healthy hepatocytes.Equation (2.1) describes the behavior of healthy hepatocytes, which proliferate according to a signal that takes on a logistic form, with maximum per-capita proliferation rate r, and homeostatic liver size K.We assume that hepatocytes die at constant per-capita rate d.They are infected by free virions at maximum rate β, according to a standard incidence function.In equation (2.5), r equates to an ecologist's maximum intrinsic rate of increase, and K is the observed equilibrium mass of a healthy liver.
2.2.Latently infected hepatocytes.Equation (2.2) describes the behavior of hepatocytes that have been infected, but are not yet actively producing virions.We assume that all infected cells initially enter a period of latent infection that lasts exactly τ days.Like healthy hepatocytes, latently infected cells die at background rate d.This assumption is based on an implicit assumption that latently infected cells are not targeted by the immune response.After τ days, latently infected cells become productively infected.Therefore, all cells infected t − τ days ago, where t is the current time, will either transition to the productive class (proportion 1 − e −dτ ) or die in the meantime (proportion e −dτ ).Thus, the delay term can be constructed directly from the underlying biology, although it is also possible to derive this term more rigorously using an age-structured approach, as in [10].
2.3.Productively infected hepatocytes.Equation (2.3) describes the behaviors of hepatocytes that are actively producing virions.The transition from the latent to productive infection has already been described.Productively infected hepatocytes die at rate a, where typically a ≥ d since, in addition to background mortality, productive cells are visible to the immune system.Indeed, both inflammatory and cytotoxic immunity appear to play definitive roles in HBV pathogenesis [12,16].Since we are considering only the case in which the immune response is sub-optimal, leading to chronic infection, we posit no immune-induced mortality beyond a constant rate of attrition.2.4.Free virions.Equation (2.4) expresses the dynamics of free virions.They are produced at rate k per productively infected hepatocyte per unit time, and disintegrate or are cleared by immune attack at rate µ.
3. Model properties.In this section we establish the positivity of model (2.1) -(2.4), and we describe the model equilibria and their local stability properties.We end the section with a computational exploration of the model's behavior.
3.1.Initial data.As in [10], initial data for the system has the form and ds.
The form for p 0 follows from the implicit solution for p(t) The case of T (0) = 0 represents the case of total loss of liver function which is mathematically trivial and of no practical interest.

3.2.
Positivity.We show below that the solution of system (2.1) -(2.4), subject to (3.1), remains bounded (and hence exists for all time) and is nonnegative for all t > 0.
The above contradictions together show that components of the solution of system (2.1) -(2.4), subject to (3.1), are nonnegative for all t ∈ [0, b).This together with the uniform boundedness of T (t) and v(t) on [0, b) imply that b = ∞.This completes the proof of the proposition.

3.4.
Stability.Here we focus on the stability of the biologically relevant equilibria and their bifurcation behavior.In particular, we show that the disease-free state, E f , is asymptotically stable when R 0 < 1, as required by biological intuition.The equilibrium representing chronic infection, E * , is treated computationally.This analysis leads to the most important biological result of the research, viz.that E * undergoes a clinically relevant bifurcation.That is, for clinically relevant parameter values this fixed point experiences a stability switch as certain parameters change.This bifurcation should be clinically evident in certain cases as the sudden onset of oscillations in viral load and liver damage.
3.4.1.Mathematical results.The local asymptotic stability of a steady state can usually be determined from the roots of the characteristic equation, det(P +Qe −λτ − λI).However, due to the form of the standard incidence function, we cannot determine the local stability of E 0 with the linearized equations.We do determine that E f is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1.For our model , and where the roots are λ 1 = −r, λ 2 = −d, and λ 3,4 are given by the solution of Observe that λ = 0 when aµ − kβe −dτ = 0. Rewriting this as Proposition 2. When R 0 < 1, the virus-free equilibrium, E f , is locally asymptotically stable.
Proof.We have that y and v satisfy the differential inequality Therefore, we can apply Theorem 3.2 of [10] and conclude that (y(t), v(t)) → 0 as t → ∞ when R 0 < 1.It is clear from equations (1.1) and (1.2) that if (y(t), v(t)) → 0 then x(t) → K and p(t) → 0. Thus, E f is globally attractive.
Proposition 4. When R 0 > 1 the virus-free equilibrium, E f , is unstable.
3.4.2.Computational results.Numerical solutions show that the stability of the equilibrium representing chronic infection is determined by r and the parameters a, β, d, k, µ, and τ , which collectively determine R 0 .The only parameter that does not have an effect is K, the homeostatic set point for liver size.In other words, the liver's ability to recover from damage, along with the virulence of the virus, determine the final disease state-a reasonable conclusion on biological grounds.When the chronic disease steady state becomes unstable, we observe the formation of a stable, periodic orbit, provided R 0 > 1. Figure 1 shows bifurcation diagrams using r, k, β, and a as bifurcation parameters.The observed bifurcations occur for reasonable values of all parameters (see below).Shown is the stability of the chronic disease steady state and the minimum and maximum values observed for the periodic orbit.Bifurcations are observed for realistic values of r, k, and β over a fairly broad range of other parameter values.However, varying a leads to a bifurcation over a narrower range of other parameters.That is, in order for the death rate of productively infected cells, a, to alter the stability of the chronic disease state, the system must already be "close" to the bifurcation point.This suggests that as long as hepatocytes retain good regenerative capacity, the liver can survive even when the immune response is very weak.It is also likely that the value of a is more important in more virulent infections.Increasing d greatly, indicating a very high rate of natural (i.e., nonpathogenic) hepatic apoptosis or necrosis, can also induce a switch in stability.The switch to a periodic orbit, for every parameter but τ , occurs as the parameters change to represent a more virulent infection or a weaker liver.
Stability is least sensitive to τ ; that is, reasonable values for τ influence stability only when the chronic steady state is already very close to unstable.However, τ has a unique effect on the bifurcation picture.The chronic steady state switches from stable to unstable, and a periodic orbit arises, as τ increases.But, increasing τ further causes the steady state to regain stability.This unusual bifurcation picture is shown in Figure 2.  1).In this model we assume r ≤ 1.0 day −1 .This value is consistent with the results of Ciupe et al. [6] and the observation that hepatocyte proliferation can be blocked in cases of severe liver damage [26].This choice is also supported by experimental liver resection in rats, in which rat livers recover to normal mass within 5-7 days after 2/3 partial hepatectomy.In humans, recovery from similar trauma occurs in 8-15 days [26].
MacDonald [24] measured normal hepatocyte life-span to be between 200 and 400 days in, whereas more recent estimates place it between 150 and 450 days [4].Therefore, we assume the half-life of healthy hepatocytes is 180 days (6 months), giving d = .0039day −1 .The latency period of infection is reported to be between 1 and 2 days [28], and was assumed by Walley et al. to be 1 day [32].We use the same value.The life-span of infected hepatocytes in chronic HBV infection is reported to range between 10 and 100 days in [30], giving a between .0693 and .00693day −1 .
A half-life of 1 day for free virus has been widely used in other models [29,30].In [32], the mean half-life of free virions was estimated to be (1.2 ± 0.6) days with a median of 1.1 days.However, the authors considered this to be an overestimate, and in [28], it was estimated that the half-life of free virus is only 4 hours.The maximum rate of daily virion production during acute HBV infection was measured to be between 200 and 1000 virions per infected cell [32].
Walley et al. [32] estimate the value of R 0 for several HBV patients in the acute phase of infection under the assumption of a 1 day half-life for free virions and a 1 day latency period.Using these estimates of R 0 and the parameter values in Table 1, we can derive an estimate for the value of β. 4. Computational exploration of lamivudine chemotherapy.Figure 3 shows the evolution of two infectious episodes from initial exposure to chronic HBV infection.In the first, r = 1.0.A disease latency period of about 50 days is followed by an acute infection phase, itself leading to damped oscillations converging to a chronic steady state.In the second, r = 0.7, and the solution converges to a periodic orbit representing chronically repeating episodes of liver damage associated with oscillations in viral load.We present a brief validation of the model by comparing the results of lamivudine treatment in a clinical study to model predictions.Lamivudine is a nucleoside analogue that interferes with viral replication within infected cells [8].Typically an oral dose of 100 mg is administered daily over the course of treatment and generally results in a reduction in viral load of 3 to 4 orders of magnitude.However, lamivudine therapy often fails because the virus develops resistance to the drug [8].
Lau et al. [20] present data from a clinical study of patients suffering chronic HBV infection but undergoing a course of lamivudine treatment, and Min et al. [27] fit a basic virus infection model to this data.In Lau's study, treatment lasted 48 weeks, and the patients were monitored for an additional 24.Viral loads were reduced by 4 to 5 orders of magnitude during the course of treatment, but quickly rebounded following cessation of therapy.As lamivudine interferes only with viral replication within the cell it is not expected to affect viral infection of healthy cells.Therefore, within our model, we assume that lamivudine treatment only reduces virion production by infected cells.The modified model follows: where γ represents the effect of treatment.Since patients had chronic hepatitis B, we assume that they were at the chronic disease steady state, E * , before the initiation of treatment.For at least one set of parameter values within the ranges presented in Table 1, namely γ = .9999,a = .011,d = .0039,β = 4.8 × 10 −5 , k = 200, K = 2 × 10 11 , r = 1, τ = 1, and µ = .693,we find excellent qualitative agreement between the model and actual data (Figure 4).

5.
Discussion.The choice of logistic growth for healthy hepatocyte proliferation is somewhat arbitrary, but its qualitative form fits well with biologically realistic liver growth.To place this formal choice on more solid biological grounds, we  posit that it represents a growth signal uniformly distributed throughout the liver, that is somehow inversely proportional to liver mass.The concentration of this signal and the cells' limited ability to respond to it-proliferation rates have upper limits set by cell biochemistry-justify, at least heuristically, the choice of a logistic term.Note that this term explicitly excludes normal hepatocyte mortality, although from a mathematical point of view that detail is irrelevant (compare equation (2.1) with (2.5)).While models using a constant infusion can fit the hepatocyte growth parameter to data, it has little real biological meaning and therefore becomes much more difficult to interpret and validate clinically.
As argued by Gourley et al. [10], the use of mass action for viral infection requires a parameter, namely β in model (1.1), that has a troublesome biological interpretation and also causes the basic reproductive number to depend on healthy liver size [10].Using a standard incidence function to represent infection dynamics eliminates this artifact and gives the infection β a clearer biological meaning, allowing a direct estimate of its value from existing literature, as we have done here.
Introduction of logistic proliferation increases the richness of the resulting dynamics.In addition to the two well known asymptotic behaviors-convergence to a virus-free equilibrium or an equilibrium representing chronic infection-we now have the possibility of convergence to a periodic orbit.Sustained oscillations are a dynamical possibility that has not been admitted by previous models.As in previous models, the virus free equilibrium in our model, E f , is both locally and globally asymptotically stable when R 0 < 1 and unstable when R 0 > 1. Clinically, then, our model predicts that infection will flair whenever aµ < βke −dτ . ( In this model, if relation (5.1) holds then the disease becomes chronic.Depending upon r and R 0 , that chronic infection may be either a steady state or oscillatory.Either asymptotic possibility can occur for biologically realistic parameter values.In the parameter region in which E * is stable, the model predicts chronic disease with a constant viral load.However, when R 0 > 1 and E * is unstable, oscillations in the hepatocyte population would arise, which would manifest clinically as repeated rounds of acute liver damage followed by periods of recovery.In fact, the model predicts that most rounds of liver damage would be profound enough to be diagnosed as acute liver failure (ALF).Clinically, ALF often occurs with little warning, and previously healthy individuals can approach death in a manner of days.ALF is characterized by wide-spread hepatocyte necrosis followed by massive immune activation and viral clearance [22].Following ALF, the liver can spontaneously recover [22,31].The periodic behavior in our model largely mimics these cycles of massive liver failure, viral clearance, and spontaneous recovery.
Since both steady state and oscillatory asymptotic behavior occur in realistic regions of parameter space, it appears possible that a steady state chronic infection can suddenly change to oscillations if parameters wander across the stability switch boundary.If this occurs, the bifurcation diagrams suggest that oscillations could rapidly evolve into increasingly wild swings between ALF and recovery as parameters continued to wander.Possible mechanisms leading to this sudden onset of oscillations include gradual degradation in the regenerative capabilities of hepatocytes (decrease in r), an increase in necrosis or apoptosis of healthy hepatocytes caused by mechanisms besides the virus (increase in d), a weakening of immune response (decrease in a), and an increase in virus virulence (increase of β or k).The first of these possibilities, reduction in the regenerative capabilities of hepatocytes, has in fact been suggested as a mechanism for the onset of acute liver failure in HBV patients.It is well established that severe liver damage impairs hepatocyte regeneration [31].Also, HBV can inhibit hepatocyte regeneration directly [35].However, high levels of proliferative activity have been observed in hepatocytes in fulminant hepatitis [33,11], suggesting that liver failure is more likely driven by hepatocyte necrosis.Changes in blood flow due to an inflammatory immune response may also cause apoptosis and ischemic damage [23] that is only indirectly associated with the virus.
In all, then, a number of mechanisms exist that could cause a change in asymptotic behavior in chronic HBV.In all cases, changes tending to promote oscillations are associated with late-stage hepatitis B. Therefore, we predict that the sudden onset of oscillations in liver damage and viral load may herald the rapid deterioration and ALF of patients who die from liver damage directly induced by HBV or the cytotoxic immune response against it.Figure 5 shows an example of such deterioration that can occur from stepwise decay of the basic hepatocyte proliferation rate, r, over the course of a 5,000-day simulation.
In addition to this steady physiological decay, an increase in HBV virulence driven by selective forces could also cause a switch in stability, leading to the oscillatory behavior and rapid deterioration in health just described, or vice versa.A future line of investigation may examine the effect of two competing viral strains to see if natural selection can induce a stability switch.It could then be determined if, and under what circumstances, the immune response can spontaneously change the dynamics back to a more stable state through its own alterations of the selective pressures.
Increasing the latency period from the time of hepatocyte infection to active virion production always increases R 0 , implying a less virulent infection.Based on a similar result for an HIV model, increasing this delay has been proposed as a novel therapeutic target [36].However, we have found that modifying the delay, τ , does not have such a straightforward effect and that modest increases in the delay can push the dynamics from a chronic equilibrium into more dangerous oscillations; further increases in the delay cause the dynamics to switch back to a less dangerous equilibrium (see Figure 2).Thus, τ influences the dynamics through more than a linear impact on R 0 , and any novel therapy would have to take this into account.We note, however, that changing τ only has a significant effect on the dynamics when the other parameters are such that the system is already "close" to the bifurcation point.Therefore, we expect that the proposed novel therapies would only have a meaningful effect in later stage diseases.
Two aspects of hepatocyte proliferation and liver recovery from damage from HBV infection could be examined using extensions of the model studied here.Following trauma, liver regeneration tends to overshoot the normal mass.Homeostatic liver size is recovered then by a small wave of hepatocyte apoptosis [26].This observation implies a possible delay in the growth signal, which could be modeled by delayed logistic growth, leading to a two delay model.
It is also well known that per-capita hepatocyte proliferation rate is a function of the amount of trauma the liver suffers.Above 2/3 HPx, more trauma yields lower per-capita proliferation rates.In fact, above 90% PHx, hepatocyte proliferation ceases as metabolic demands on surviving hepatocytes become too great [31,11].On the other hand, per-capita proliferation also declines for resection masses below 2/3 PHx [31].So, maximum per-capita hepatocyte recovery occurs when the surviving liver is 1/3 its normal mass.Our model could be improved by altering the proliferation function to accommodate this more complex growth response.
We have modeled the latency period from infection to virion production with a single discrete delay.In reality, any such delay will vary between individual cells, yielding a distributed delay; the discrete delay we have chosen should be understood as the mean of the actual distribution.Confirming that the behavior and predictions of our model are robust with respect to distributed delays could be done using techniques similar to those in [5].
The standard incidence infection term with logistic hepatocyte growth can potentially generate rich dynamics, rivaling those of ratio-dependent predator-prey models [13,18] or simple epidemic models [3,14,15]; we will explore this potential in the near future.These characteristics, together with the time delay, present many challenging mathematical questions on the global dynamics of model (2.1) -(2.4) with complexity surpassing that of [9].
One major weakness of our model lies in how the immune system is handled.In particular, we ignore the adaptive nature of the immune response.Both the cytotoxic response, represented in part by a, and the humoral, represented by µ, are linear in y and v, respectively.A more realistic form would allow both cytotoxic attack of infected cells and humoral clearance of free virus to increase nonlinearly with the number of productively infected cells and free virions, respectively.In fact, our simplification of the immune response is probably the primary reason why our model fails to replicate the most common clinical course seen in HBV patientsnamely acute, self-limiting disease.Therefore, the next generation model should include a more realistic model of adaptive immunity.

Figure 4 .
Figure 4. Comparison of model results to clinical results for a course of lamivudine therapy.

Figure 5 .
Figure 5.An example of a switch in stability resulting from a stepwise reduction in r over the course of 5,000 days (13.7 years).