An age-structured model of hepatitis B viral infection highlights the potential of different therapeutic strategies.

Hepatitis B virus (HBV) is a global health threat, and its elimination by 2030 has been prioritised by the World Health Organisation. Here we present an age-structured model for the immune response to an HBV infection, which takes into account contributions from both cell-mediated and humoral immunity. The model has been validated using published patient data recorded during acute infection. It has been adapted to the scenarios of chronic infection, clearance of infection, and flare-ups via variation of the immune response parameters. The impacts of immune response exhaustion and non-infectious subviral particles on the immune response dynamics are analysed. A comparison of different treatment options in the context of this model reveals that drugs targeting aspects of the viral life cycle are more effective than exhaustion therapy, a form of therapy mitigating immune response exhaustion. Our results suggest that antiviral treatment is best started when viral load is declining rather than in a flare-up. The model suggests that a fast antibody production rate always leads to viral clearance, highlighting the promise of antibody therapies currently in clinical trials.

Introduction

Hepatitis B virus (HBV) is a major global health problem with a high risk of causing chronic infection and death due to cirrhosis and liver cancer. Despite pioneering use of recombinant DNA technology to create a safe and cheap prophylactic vaccine against HBV[1,2] it is not universally deployed allowing the virus to infect 1.5 million people each year[3]. These add to the very significant burden of the 296 million chronically infected individuals, that are the consequence of the failure of their immune systems to clear a primary infection[3]. Treatment of these patients relies on generic replicase inhibitors that rapidly elicit resistance and interferon to boost antiviral responses. These treatments are largely ineffective leading to 820,000 additional HBV-related deaths annually, the largest cause of liver cancer worldwide[3]. HBV infected cells produce a significant excess of non-infectious particles, such as subviral Dane particles (SVPs) that occur at 1000–100,000-fold higher concentration than the infectious virion, allowing them to act as immune system decoys by sequestering antiviral antibodies[4,5].

HBV infected individuals that mount an adequate immune response tend to clear the infection. During the early stages of the infection the first contribution of the immune response comes from innate immunity which reduces viral spread and facilitates an adaptive immune response. Adaptive immunity is made of two components, cell-mediated immunity (CD4 and CD8 T cells) and humoral immunity (antibodies). CD4 T cells, also known as helper T cells, assist the activity of other immune cells by releasing cytokines. CD8 T cells are not only responsible for killing of infected cells but also induce the noncytolytic “cure” of such cells while antibodies neutralise virus particles and prevent infection of cells[6,7].

Mathematical modelling provides new insights into the various aspects of viral infections and the impact of the immune response on their clearance. Ciupe et al. presented a model to study the role of cytolytic and noncytolitic immune responses and the time lag associated with effector cell activation and expansion during an acute HBV infection. It is hypothesised that cured cells and their progeny are less likely to get infected, preventing reinfection[8,9]. Later, Ciupe et al. looked into the dynamics of antibodies and showed that retaining a strong cell-mediated immune response is crucial for the control of early infection in unvaccinated individuals[5]. Fatehi et al. developed a mathematical model to take into account contributions from innate and adaptive immune responses, as well as cytokines. The model investigates the roles of different components of the immune response on viral dynamics[10]. These models are based on systems of ordinary differential equations (ODE) or delay differential equations (DDE). Nelson et al. presented an age-structured model of human immunodeficiency virus (HIV) infection to study the impact of variations in the virion production rate and the death rate of infected cells over the course of the infection[11]. Although age-structured models are more complicated, they can provide more realistic dynamics[12].

Experimental studies have shown that persistent stimulation of effector cells may result in immune impairment, e.g., immune exhaustion[13-15]. In HBV infections, persistent antigen presentation by infected cells and exposure to high antigen loads plays an important role in CD8 T cell exhaustion[16]. In order to analyse the impact of T cell exhaustion on viral dynamics, Johnson et al. introduced a variable that captures the antigenic stimulus, called the level of exhaustion[17]. They assumed that T cells are inactivated dependent on the level of exhaustion and modelled it as a Hill function with a half-maximal constant called the exhaustion threshold. They showed that the exhaustion threshold has a significant impact on the ability of the immune response to control an infection[17]. Later, Conway and Perelson included T cell exhaustion into an HIV infection model and showed that the strength of cytotoxic T lymphocytes in killing productively infected cells, and the level of latently infected cells, determine the post-treatment outcome of the infection[18].

We recently introduced a model of intracellular HBV infection dynamics and used it for comparative analysis of different therapeutic strategies[19]. The model reveals a two-phase behaviour in the release of non-infectious SVPs. Shortly after infection, a cell starts secreting SVPs. When the first intact virions are released, after 90 hours post-infection, SVP secretion stalls until the onset of a second secretion phase. In order to analyse the impacts of this behavior on the immune response dynamics, and in particular their interaction with anti-HBsAg antibodies, we have developed an age-structured model using infection kinetic parameters derived from our intracellular model. This model has been parameterised such that its outputs fit data recorded from patients undergoing an acute infection. We have adapted to the scenarios of a chronic infection, immune clearance of the infection, and infection flare-ups via variation of the immune response parameters. The role(s) of cell-mediated and humoral immunity in an acute HBV infection, including the effects of immune response exhaustion, have been analysed. The impacts of direct-acting antivirals (DAA) treatments targeting various steps of the virus life cycle are compared to blocking T cell exhaustion, which is known as exhaustion therapy. Optimal treatment regimes, such as the most effective treatment start times for the reduction of viral load, have been identified. Our analysis highlights the potential benefit of antibody therapies currently in clinical trials, and allows modelling of the impacts of DAAs under development.

Results

Model derivation

In order to analyse HBV infection dynamics with emphasis on  the roles of the immune response, we developed a detailed intercellular age-structured model of HBV infection. This model includes uninfected target cells, T(t), which are assumed to be created at rate and to die at a rate d[18,20]. Each individual infected cell is characterised by its age, a, since infection. Therefore infected cells are denoted as I(a, t)[11]. An infected cell produces complete virions, which are infectious, and incomplete particles which occur in three distinct forms: RNA containing particles; empty virions; and subviral particles (SVPs) which are either filamentous or spherical[4]. The release dynamics of these particles from an infected cell has been studied precisely[19]. The variable indicates the number of complete virions which infect target cells at a rate . All types of incomplete particles are non-infectious and present HBV surface antigens (HBsAg), enabling them to act as decoys for the immune response by depleting HBsAg-specific antibodies, A(t). We therefore only include the variable into the model, which represents the total number of incomplete particles[4]. The production rate of these particles depends on the age of an infected cell. The functions and show the production profiles of complete and incomplete particles, respectively, from infected cells of age a[11]. We use the functions that are presented in Fig. 2c of[19] as the base functions and (see Supplementary Fig. S1 online). Since it has been shown that changing the intracellular model parameters will change the total number of released particles, we scale functions and by and , respectively. We assume that and are changeable parameters and fit them to patient data. Complete and incomplete particles are cleared, based on their half-life, at rates and . HBsAg-specific antibodies, which are a component of the adaptive immune response called humoral immunity, are produced at rate proportional to antigen load, and are degraded at rate . We also add a logistic growth to the antibody equation, which is maintained through homeostatic proliferation of memory B cells after infection. Antibodies can bind to complete and incomplete particles at rate to create complete virion-antibody complexes, , and incomplete particle-antibody complexes, , respectively. These complexes disassociate at rate , or degrade at rate [5].

The other component of an adaptive immune response is cell-mediated immunity. This occurs via the action of CD8 T cells, also referred to as effector cells, E(t). These cells kill or cure infected cells. Since it has been hypothesised that cured cells lose their resistance to productive infection at a slow rate of the order of per day[8], we model these two impacts in one reaction, where effector cells remove infected cells at rate . Effector cells target infected cells in response to the display of viral peptide-major histocompatibility complex (MHC) on the cell surface, which is dependent on the expression of viral proteins in the cells. As the intracellular model has shown, these viral proteins are produced and secreted shortly after infection. We therefore assume that the removal rate of the infected cells is age-independent[19]. In addition, we assume infected cells die at rate due to infection or the action of the innate immune response, which is not included in the model directly[5,10]. Effector cells are assumed to be produced at rate and removed at rate in the absence of infection[8]. During infection, proliferation of effector cells happens in an infected cell density dependent manner with a time delay, i.e. it depends on the density of antigen, with maximum proliferation rate , and on , the level of infected cells at which proliferation is half-maximal[8,17,21]. It has been argued that persistent exposure to high antigen loads is associated with CD8 T cell exhaustion[16]. Therefore, functional effector cells are lost at maximal rate , depending on the level of exhaustion, Q, which is implemented as a Hill function with coefficient n and half-maximal constant [17]. The level of exhaustion, Q, is measured by integrating over the antigenic stimulus () times the parameter which is called “blockade parameter”, and reduces exponentially with coefficient , which is the rate of immune response recovery from exhaustion. The reduction of the blockade parameter is assumed to simulate the blockade of interaction between the inhibitory receptor programmed death 1 (PD-1) and its ligand (PD-L1) which restores CD8 T cell function[22].

With the above assumptions, the model, as illustrated in Fig. 1, takes on the form:where and .

Figure 1

The components and inferred interactions of the age-structured model of the immune response to an HBV infection (1), enabling comparison of different treatment options. Blue and cyan circles show uninfected (T) and infected (I) cells, respectively, green circles indicate the immune response (antibodies (A) and effector cells (E)), and gray circles show complete () and incomplete () viral particles. Pink circles represent complete virion-antibody () and incomplete particle-antibody () complexes. Yellow shows the level of exhaustion (Q). Double arrow-headed lines indicate natural clearance, bar-headed and single arrow-headed lines indicate destruction and production/proliferation, and forward/backward arrows represent binding/unbinding events.

All parameter descriptions and values adapted from the literature are presented in Table 1 together with a pointer to the references they have been adapted from. The remaining parameters are estimated by fitting the model (Fig. 1).

Table 1

Overview of model parameters. Parameters that are fixed in the model are either adapted from the references indicated, or are fitted.

Parameter	Description	Value	Reference
α	Effector cell growth rate	mlday^-1 (varies)	–
A_m	Antibody carrying capacity	3.4×10^12perml	[5,26]
β	Infection rate	mlday^-1 (varies)	–
δ	Infected cells death rate	day^-1 (varies)	–
d	Death rate of target cells	0.01day^-1	[20,25]
d_A	Antibody degradation rate	0.033 day^-1	[5]
d_c	Virus clearance rate	0.67 day^-1	[5]
d_E	Removal rate of effector cells	0.5 day^-1	[9]
d_i	Incomplete particle clearance rate	0.67 day^-1	[5]
d_q	Reduction rate of exhaustion level	0.1 day^-1	[17]
d_x	Complex degradation rate	2.7 day^-1	[5]
κ	Growth rate of exhaustion level	1 day^-1	[17]
k_b	Antibody dissociation rate	10 day^-1	[5]
k_f	Antibody binding rate	10^-10mlday^-1	[5]
λ	Target cells production rate	4×10^5perml	[8,23,24]
λ_E	Effector cell production rate	10 perml	[9]
μ	Removal rate of infected cells by effector cells	mlday^-1 (varies)	–
n	Hill function coefficient for effector cell exhaustion	3	[17]
ϕ	Hill function half-maximal constant for effector cell growth	10^6perml	Fitting
p_A	Antibody production rate	moleculeday^-1 (varies)	–
q_c	Hill function half-maximal constant for effector cell exhaustion	10	[17]
ρ_1	Rate modifying for virus production rate	varies	–
ρ_2	Rate modifying for incomplete particle production rate	varies	–
r_A	Antibody proliferation rate	day^-1 (varies)	–
τ	Effector cell growth time delay	day (varies)	–
ξ	Effector cell maximal exhaustion rate	1day^-1	Fitting

Within-host viral dynamics

The model (Fig. 1) was fitted to viral load data from 6 patients who were recorded during the acute stage of infection (Fig. 2)[5,8]. In this work, we also let the initial viral load () be a variable. In patient number 2 (Fig. 2 P2) the peak in viral load occurs earlier than in other patients, which shows the initial viral load could be higher in this patient. We observe that for patient numbers 1, 3, 4, 5 and 6, virion per ml as suggested by Ciupe et al.[5], but for patient number 2 virions per ml. The parameter values derived from fitting of are presented in Table 2. Our model captures the essential features of viral loads in all patients, including the first rapid decline in the level of viral load after the peak, followed by a slower decline, i.e. the biphasic viral decay after the peak[8]. In all patients viral load eventually decreases to zero, matching the clinical outcomes in these patients with acute HBV infections.

Figure 2

Paramaterising the model against patient data. The fitting of as given by model (Fig. 1) (red line) to patient data () indicates a biphasic decline in viral load and sometimes a slower decline in the level of infected cells compared with virions (green line).

Table 2

Parameter best estimates

Patient	β×10^-9	μ×10^-4	δ	ρ_1	ρ_2	p_A×10^-4	r_A	α	τ	V_c(0)	RSS
1	0.35	1.00	0.0130	8.50	2.50	3.00	0.366	2.20	6.50	0.33	0.29
2	3.57	2.00	0.0210	3 .00	0.70	1.00	0.550	1.50	4.50	10.00	0.14
3	1.20	0.90	0.0170	3.20	0.70	1.00	0.391	1.80	5.50	0.33	0.29
4	1.70	2.00	0.0400	3.00	0.70	1.00	0.370	1.80	5.00	0.33	0.24
5	1.30	1.80	0.0180	3.50	0.70	1.00	0.430	1.80	6.00	0.33	0.22
6	0.35	1.00	0.0130	8.50	2.50	3.00	0.366	2.20	8.70	0.33	0.13
Median	1.25	1.40	0.0175	3.35	0.70	1.00	0.381	1.80	5.75	0.33	–
Mean	1.41	1.45	0.0203	4.95	1.30	1.67	0.412	1.88	6.03	1.94	–
Std	1.19	0.54	0.0101	2.76	0.93	1.03	0.072	0.27	1.49	3.95	–

The average fraction of infected cells at the peak of infection is around 74.5%, with a range of 60%-82%. This is in good agreement with previously estimated values[8,9], and the experimental observation that more than 75% of hepatocytes were hepatitis B core protein antigen positive () in infected chimpanzees[6]. Previous models suggested that the level of infected cells has a biphasic decline and decreases faster than the viral load[8,9]. However, it has been observed that the level of cccDNA (infected cells) has a slower decline compared with the level of free virus and remains detectable for more than two years after infection[27]. Our model predicts that infected cells decline slower compared with viral load (green lines Fig. 2), in good agreement with Fig. 2 in Wieland et al.[27]. We fitted parameter values from Table 2 to a Gaussian distribution, inferred 100 parameter sets and plotted the outcome for the 100 inferred parameter sets (see Supplementary Fig. S2 online). Consistent with the fact that these patients are acute cases, we observe acute infection in 60% of the simulations, a periodic solution in 25% of simulations, while the remainder exhibit a stable chronic infection steady state with a low level of free viruses (around virions per ml).

Adaptive immune response dynamics

The level of serum alanine aminotransferase (ALT) is used as proxy for the dynamics of effector cells in HBV infection[8,9]. It has been observed clinically that the peak in serum ALT occurs after the peak in viral load and the time lag between these two peaks is 3–4 weeks[8,28]. Our model also predicts that the average time lag in these 6 patients is 26.18 days (Fig. 3) which is in good agreement with previously reported values[8,28]. We observe this delay in all patients (Fig. 3) and in our model (Fig. 4c). This is due to the time lag in the production and activation of the effector cells and their propagation from lymph node to the tissue. This delay creates a window of therapeutic  opportunity to start treatment with maximal effect while viral load is decreasing and effector cells are approaching their maximal value. To provide a more realistic dynamic of effector cells we included their exhaustion into our model. Supplementary Fig. S3 online indicates that in each patient the onset of exhaustion coincides with a step increase in the level of effector cells. However, before the peak in exhaustion level, the immune response controls the infection, so that the level of exhaustion starts declining before the peak of effector cells is reached.

Figure 3

Time course of the immune response to the HBV infection in terms of effector cell and anti-HBsAg levels. The maximum in effector cells (black curve) occurs 3–4 weeks after the peak in viral load. The red lines show the best fit to patient data, and the black dashed lines show the level of anti-HBsAg in mIU/ml unit as predicted by the model (Fig. 1).

Figure 4

Stability analysis of the disease-free and chronic infection states, exemplified for different types of infection dynamics: acute infections with immune clearance, and chronic infections with and without infection flare-ups. The parameters used are the median values from Table 2. The gray and yellow areas in (a) and (b) indicate the regions where the disease-free (DF) and chronic infection (CIs) steady states are stable, respectively. Pink is the region where the system shows a stable periodic solution around a chronic infection steady state (CIp). Solid and dashed lines indicate the boundaries of the steady-state and Hopf bifurcation, respectively, and “fH” pinpoints the location of the fold-Hopf bifurcation. Stars indicate the points that are used for model simulations. (c) and (d) correspond to an acute infection. , the maximum proliferation rate of effector cells, in (e) and (f) is reduced by 28% (), indicating periodic oscillations around the chronic infection steady state (hepatitis flare). In (g) and (h), is reduced by 56% (), illustrating the case of a chronic infection.

Figure 3 indicates the level of antibodies against HBsAg (anti-HBsAg) in these patients. Anti-HBsAg levels () are considered as negative[5,29]. In unvaccinated patients, the level of antibodies is under 10 mIU/ml while they are still classed as infected, i.e. their HBV-DNA test is positive[26]. Figure 3 shows that in all patients the level of antibodies is well below 10 mIU/ml during the peak in viral load and the first phase of viral decline. Later, when the level of virus has decreased significantly and effector cells are reduced to their basal level, it starts increasing but is still under 10 mIU/ml[26].

Now we study the impact of varying immune response parameters on viral dynamics. We focus on parameters (growth rate of effector cells) and (the rate at which functional effector cells are lost due to the exhaustion) in the cell-mediated immune response which play crucial roles in the control of the infection[10,28]. The other important parameter is the antibody proliferation rate ()[5]. Figure 4 shows the stability and instability regions of the disease-free and chronic infection states as a function of these parameters. This figure shows that an increase in always stabilises the disease-free state. However, when is low (), increasing or decreasing cannot stabilise the disease-free state, However, it destabilises the chronic steady state due to a Hopf bifurcation, so one observes stable oscillations. Biologically, this is known as “flare-up”, where through the interactions between viral infection and immune response we observe a periodic behaviour[10,30,31]. Examples of different types of infection dynamics are shown in Fig. 4 where Fig. 4c and d indicate acute infection, Fig. 4e and f show a chronic infection with flare-ups, and Fig. 4g and h show a chronic infection without flare-ups (stable chronic infection steady state). Supplementary Fig. S4 online illustrates the impact of the time delay in effector cell activity () and of the removal rate of infected cells by effector cells () on viral dynamics. Supplementary Figs. S4a and b online indicate that a decrease in increases the stability range of the disease-free state with respect to and , respectively. That is, following a decrease in , the disease-free steady state will become stable for lower values of and higher values of . Supplementary Fig. S4c online shows a similar effect with respect to . This figure shows that an increase in results in an increase in the stability range of the disease-free state.

Effects of current clinical therapies

There are currently two approved anti-HBV therapies for clinical use. Interferon (IFN)-based therapy and nucleot(s)ide analogues (NAs)[32], but more therapeutic strategies will be required to meet the Global WHO Challenge to make HBV infection a treatable condition by 2030[16]. We recently compared current DAAs with several possible future treatment options using our intracellular model[19]. The latter include capsid assembly modulators (CAMs), especially based on heteroaryldihydropyrimidines HAPs[33]. Recently we showed that both HCV and HBV regulate the assembly of their nucleocapsids around ssRNA forms of their genomes (e.g. the pgRNA of HBV) at least in part via sequence-specific core protein-genomic RNA interactions at multiple sites, termed Packaging Signals (PSs)[34-36]. CAMs are at an advanced stage of development[37-39], whilst the vital roles of PS-core interactions in HBV nucleocapsid assembly, and their targeting by small ligands are at a much earlier stage. The comparative analysis reveals strong potential synergistic effects between them. These treatment options reduce the viral production rate and their effect is modelled by a factor , where is the drug efficacy and its value is constant from the start of the treatment until its removal[20,40]. This method has also been used in the context of other viruses[41-45]. IFN-based therapy and CAMs also block de novo infections[10,46]. Mathematically, one can represent this effect by a modified transmission rate , where is the drug efficacy[10,19]. Here we refer to these types of treatments as antiviral therapy. The other suggested form of therapy is the recovery of the CD8 T cells from exhaustion (exhaustion therapy)[16]. The impact of this treatment option can be modelled by a modified blockade parameter , where is the efficacy of the treatment[22]. The new equations take the following form:The total drug effectiveness , defined as , allows us to determine a critical drug efficacy, , where the infection gets cleared for [10,47].

IFN-based therapy, which results in higher rates of hepatitis B e antigen (HBeAg) and HBsAg loss compared with NA therapy, is usually administered for around 48 weeks[32]. Thus, 48 weeks after the start of treatment we remove its effect on the model to determine the efficacy required to avoid a viral rebound after the removal of treatment. This will help us to find an optimal efficacy that a new therapy should have to be effective, thus addressing an urgent need identified by Lok et al.[32]. Figure 5 shows the impact of starting antiviral therapy at different times post infection. In Fig. 5a treatment start is 130 dpi, while in Fig. 5b it is 150 dpi. These figures indicate that in a region where the system shows a periodic behaviour, an earlier treatment start during the declining phase of the viral dynamics is more effective than a treatment start during the increasing phase or around the peak of a flare-up (Figs. 5c and d, verified explicitly for treatment start one week before, and after, and the top of the peak). This indicates the co-existence of the disease-free steady state and stable periodic oscillations around the chronic infection steady state, so that depending on the initial conditions the system can either converge to the disease-free case or exhibit a periodic behaviour. Thus, starting treatment at different times can lead to different outcomes. Supplementary Video S1 online shows the minimal total efficacy () that is required to clear the infection following 48-weeks of therapy starting at various times between 50 dpi (representing an early detection of the infection) to 160 dpi (representing a late detection of the infection). It indicates that treatment start at 50-80 dpi (around a month before the first viral peak) is most effective. It also shows that treatment start at 90-120 dpi (close to the first viral peak) is not that effective, in particular when the chronic steady state is stable. After the first viral peak, if the chronic steady state is stable, the treatment start time is not important as the infection is in a stable state. However, in the case of a periodic solution, starting treatment when the viral load is declining requires a treatment of lower efficacy than would be necessary to achieve the same outcome for a treatment start in the increasing phase. These results indicate that a decline in the level of free virus in a patient, can be a sign of an acute infection, or correspond to the declining phase before a flare-up. Therefore, it is crucial to start the treatment, rather than wait to start it after the flare-up. Supplementary Video S2 online shows the minimal efficacy required to clear the infection following 48 weeks of exhaustion therapy starting at in between 50 dpi and 160 dpi. Surprisingly, it shows that this treatment is not as effective as antiviral therapy targeting the viral life cycle and the treatment start point has no significant impact on the outcome. Moreover, our model does not show any significant synergistic effects between these two treatment options. This can be seen from the fact that this treatment option is only effective in a small area of the - parameter plane, corresponding to an area in which antiviral therapy targeting the viral life cycle is also effective even for low drug efficacy.

Figure 5

Start of antiviral therapy when viral load is in a declining phase is more effective. (a) and (b) show the minimal total efficacy () that is required to clear the infection following 48-weeks of therapy starting at 130 and 150 days post infection (dpi), respectively. The white area indicates a stable disease-free state. The red solid curves indicate the onset of the steady-state bifurcation, whereas the red dotted lines of the Hopf bifurcation of the chronic state. The black hatched area indicates the region where an efficacy of is ineffective. The black star (indicated with an arrow) shows the point at which the simulations ((c) and (d) ) have been performed. Magenta and brown vertical lines indicate the start and end of treatment, respectively. The dashed and dotted lines, represent the effects of treatments on viral load, for treatment starting 130 or 150 dpi, respectively.

In Fig. 5 we consider a constant treatment efficacy for the full duration of the treatment, as is typically done in viral infection models[20,40-45]. However, in Supplementary Figs. S5 and S6 online we study the impact of a periodic administration of antiviral treatments of different intensities (efficacies), including also the option of zero efficacy to model the absence of treatment in one of the treatment periods, to find possible synergistic effects for alternating treatments. We consider treatment starts at 130 and 150 dpi, respectively, and different combinations of treatment periods, i.e. a treatment with total efficacy in period 1 is followed by a treatment with total efficacy . Then this scheme is repeated for 48 weeks, which is the recommended time for administration of an HBV therapy[32]. Supplementary Fig. S5 online indicates that alternating antiviral therapy with a total efficacy of 90% for 4 weeks, with 4-week intervals of no treatment over 48 weeks leads to clearance of the infection. To understand the dynamical reason for this, we checked viral load 28 days post treatment stop. We observed that it is below the level of virus in the drug-free (control) case at 186 dpi (130+56) and is in a declining phase. Although for treatment start at 150 dpi the results are more chaotic (see Supplementary Fig. S6 online), cases leading to clearance of the infection show similar behaviour as for the 130 dpi case. Supplementary Fig. S7 online shows the impact of a treatment start at 150 dpi, with periodic administration of antiviral treatments of different efficacy in 4 weeks intervals. This figure shows that in successful cases antiviral treatment is administered during the declining phase of viral load. This is consistent with the main conclusion from our modelling of constant treatment efficacies above, that also shows that the outcome of a treatment depends highly on a start time in the declining phase of viral load. Although Supplementary Figs. S6 and S7 online indicate that drug administration in various efficacies and periods can sometimes lead to clearance of the infection, establishing these time periods and efficacies in a patient is not practicable in the clinic. To address this issue, we studied the impact of a treatment start at each day in one period of the viral load dynamics (i.e. each day from 130 to 212 dpi, covering 83 possible start times). In each case, we either consider a constant efficacy, or a periodic efficacy where the treatment is administered for 28 days and then withheld for 28 days. We chose an 28-day periodicity as this is optimal according to Supplementary Figs. S5 and S6 online. This shows that for a total treatment efficacy of 99% a constant treatment will always leads to viral clearance. However, if the total efficacy is only 90% then a constant treatment is successful only at 52 treatment start times i.e., the chance of having a successful outcome is 62.65%. By contrast, a periodic treatment would increase the chance of success to 96.38%. It is worth noting that all currently available treatments for HCV have an efficacy of more than 99%[20]. Also CAMs and the PS-targeting drugs for HBV, which are under development, have shown promising outcomes with efficacies of more than 95%[19]. Thus, our analysis suggests that a non-periodic treatment option, in line with those previously administered for HBV and other viral infections, appears to be a reasonable strategy. However, different strategies should be tested in clinical trials, as our modelling indicates that under certain circumstances periodic treatment options can be superior.

Monoclonal anti-HBsAg antibody drugs are currently being investigated in clinical trials as a treatment against HBV infections[48-50]. To study the impact of this treatment option we extended our model, modelling these antibodies as an influx into the system (see Supplementary Section S1 online). Supplementary Fig. S8 online indicates that the impact of monoclonal antibody therapy is independent of the treatment start time, leading to clearance of the infection in both the declining and the increasing phase.

Discussion

In this paper we derive an age-structured model for the immune response to HBV infections. It is based on the particle release profiles from a detailed intracellular model[19]. Our age-structured model focuses on two components of the adaptive immune response known as humoral (antibodies) and cell-mediated (CD8 T cells) immunity. To provide more realistic viral dynamics and robust predictions, we fitted our model to patient data. Patient data indicate a biphasic decrease in viral load which is well captured by the model. Previous models have predicted that the number of infected cells declines faster than viral load[8,9], but it has been shown clinically that the level of cccDNA (infected cells) has a slower decay compared with viral load[27]. Interestingly, our model reflects this trend. Regarding the dynamics of the immune response, our model shows that the peak in the level of effector cells occurs around 26 days after the peak in viral load, which is in good agreement with previously reported values[8,28]. The antibody level stays well under 10 mIU/ml during the peak and the first phase of viral decline. Later, it starts increasing, but it remains under 10 mIU/ml[26].

We also performed a stability analysis of the model to study the impact of varying parameter values on the outcome of the infection. The model suggests that increasing (the proliferation rate of antibodies) always stabilises the disease-free state. When is low, decreasing (the rate at which functional effector cells are lost due to exhaustion) or increasing (expansion rate of effector cells) will not stabilise the disease-free state, but the chronic steady state loses its stability and one observes stable oscillations[10,30,31]. Therefore, our model illustrates the importance of the humoral immunity in viral clearance.

Current therapeutic strategies either block the formation of new virions and viral entry (antiviral therapy) or recover the exhausted CD8 cells (exhaustion therapy). We used our model to compare the impacts of these different treatment options. Our results suggest that exhaustion therapy is not as effective as antiviral therapy, and that the start of treatment does not have a significant impact on the outcome of exhaustion therapy. However, the timing of treatment start plays a crucial role in the outcome of antiviral therapy. Our model shows that around one month before the peak in viral load would be the best time to start antiviral therapy, illustrating the importance of identifying cases early into an infection. Since this is not always practicable, we also analysed the impact of a treatment start later during infection. Whilst treatment start time is not significant in the context of stable chronic steady states, it matters for periodic solutions, where a treatment start during the declining phase of viral load is more effective. As declining levels of free virus in a patient can indicate the declining phase before a flare-up, it is crucial to start treatment immediately rather than wait until after a flare-up. Our model moreover illustrates that antibodies play an important role in suppression of viral load, supporting strategies to use monoclonal antibodies as therapy against chronic Hepatitis B viral infections[48-50], and our model suggests that they should be highly efficient ways of combating chronic hepatitis B.

One critical hurdle to such developments has been the restriction of HBV infection to people and other higher primates, restricting the rapid development of new drugs in animal models. The development of humanised transgenic mice, available with and without a human immune system, will hopefully permit much faster experimental development[51,52]. This model and others will be useful for the evaluation of any novel therapeutic strategies being developed.

Methods

Patient data

Our study uses patient data comprising 6 patients whose their viral loads were recorded during the acute stage of infection and are reported in Ciupe et al.[5,8], Webster et al.[28] and Whalley et al.[53]. Patients 1, 2, 4, and 6 are female with an age range of 37 to 71. The data are published in the Supporting Information of Ciupe et al.[5].

Parameter estimation

It has been estimated that the total number of liver cells in an adult is equal to and only 60% of them are hepatocytes[8,23,54,55]. We assume that there are 3 liters of serum in an adult, containing complete and incomplete particles[24]. Therefore, . The death rate of target cells is estimated to be [20,25], so . As in unvaccinated patients with positive HBV-DNA test the antibody level is under [26], we assume that . For parameters (the level of infected cells at which growth is half-maximal) and (the rate at which functional effector cells are lost due to exhaustion) various fixed values have been determined previously[17,18]. We thus set and . These values provide a good fit to the viral load data. All parameter values adapted from the literature are presented in Table 1 together with a pointer to the reference they have been adapted from. The remaining parameters are estimated by fitting the model (Fig. 1) to the viral load data using a modified version of DKLAG6 for the numerical simulations of the model (Fig. 1) in Fortran (see www.radford.edu for more detail)[56] and an implementation of the Nelder-Mead algorithm in Fortran (see people.sc.fsu.edu for more detail)[57]. More details about the method are provided in the supplementary materials online (see Supplementary Section S2 online).

Supplementary Figures.

Supplementary Video 1.

Supplementary Video 2.