Advanced psoriasis treatment: Development and optimization of resorcinol-loaded nanoemulgel with lavender oil

1.   Introduction

Hepatitis B virus (HBV) infection is a major public health burden with high endemic areas in South East Asia, China, and sub-Saharan Africa ^[1]; and approximately 240 million chronically infected people worldwide ^[2]. HBV infects a subset of liver cells (i.e. hepatocytes) ^[3] and can lead to either acute or chronic disease. About 90% of perinatally and 20-30% of childhood acquired HBV infections become chronic ^[4], while healthy adults clear the infection in 95% of the cases ^[3]. Severe complications, such as liver cirrhosis and hepatocellular carcinoma (HCC), follow chronic infections ^[3].

Chronic HBV spans five distinct disease stages which are built around the dynamics of a serological marker called hepatitis B e antigen (HBeAg), a secretory protein that is not required for viral replication or hepatocyte infection ^[5] and has been described as a downregulator of the cellular immune response, therefore acting as an immune tolerogen ^[6]. The stages are not clearly separated or sequential ^[2]. The first four are distinguished by the presence/absence of HBeAg and by disease (hepatitis), as described below ^[2]. The first phase, called HBeAg-positive infection (formerly known as the immune tolerant phase), lasts between 10–30 years ^[4,2]. Viral DNA levels are high and HBeAg is detectable, while alanine aminotransferase (ALT) (a marker of liver disease) is normal, indicating lack of liver cell damage ^[4,7]. The second stage, called HBeAg-positive hepatitis (formerly known as immune clearance or immune active phase), is characterized by high and fluctuating viral DNA levels, the presence of HBeAg, elevated ALT and moderate to severe liver damage ^[2,4]. The third phase, called HBeAg-negative infection phase (formerly the inactive carrier phase), is marked by HBeAg seroconversion, i.e. hepatitis B e-antibody (HBeAb) production and subsequent HBeAg loss, low viral DNA levels, normal ALT and no liver disease ^[2,8]. The fourth phase, called HBeAg-negative hepatitis, is characterized by undetectable HBeAg, detectable HBeAb, moderate to high viral DNA levels, elevated ALT and liver disease ^[2]. The fifth phase, which is no longer determined by HBeAg, is called hepatitis B surface antigen (HBsAg) negative phase. It is marked by normal ALT levels and usually undetectable viral DNA levels ^[2]. Reaching this phase before the onset of cirrhosis significantly reduces the risk of liver damage (either by cirrhosis or HCC) ^[2].

For chronically infected HBV patients two main groups of HBV drugs are available: interferon-$\alpha$ (IFN$\alpha$) or its pegylated form PegIFN$\alpha$, and five nucles(t)ide analogues (NAs): lamivudine, telbivudine, adefovir, entecavir, and tenofovir ^[9]. IFNs have both immunoregulatory and antiviral effects ^[9]. NAs have only antiviral effects, such as inhibition of HBV replication ^[9,10]. The goal of HBV treatment is to reduce the risk of disease progression and HCC development ^[2]. The optimal treatment endpoint is seroclearance of HBsAg, associated with a very low risk of viral relapse and progression to HCC ^[2,11].

The severe side effects associated with IFN$\alpha$ and PegIFN$\alpha$ limit the amount of time when they can be administrated: usually 48 weeks, but in certain cases up to 96 weeks ^[2,12]. Therefore, the more widely used treatment options are NAs. Older NAs caused viral resistance, but this risk was significantly reduced by the newest generation of NAs: entecavir and tenofovir, which allow for long-term (indefinite) treatment ^[2,13,14,15]. In spite of these advances in NAs therapy, HBsAg seroclearance is reached in only 1% of treated patients ^[13]. Hence HBeAg seroconversion, in addition to viral remission to an undetectable viral DNA level, is often considered a more realistic treatment endpoint ^[11,16]. In particular, for HBeAg-positive patients that undergo HBeAg seroconversion and viral remission during treatment, stopping therapy after some consolidation phase is recommended ^[2]. A systematic review of NA treatment studies has reported that after initial seroconversion about 95 (92, 88)% of patients stay HBeAg-negative for 6 (12, 24) months post treatment cessation, and that about 73 (62, 53 51)% of initially HBeAg-positive patients remain in viral remission 6 (12, 24, 36) months after the end of therapy.^[13].

The role of HBeAg seroconversion has been investigated in several clinical studies. It has been reported that patients that undergo HBeAg seroconversion have a better prognosis than those that are consistently seropositive, such as slower disease progression and regression of the fibrosis ^[17,18]. Other studies have shown that persistently HBeAg-positive patients have a higher risk of developing HCC ^[19], and liver cirrhosis ^[20,21,22]. Additionally prolonged HBeAg-positive hepatitis or higher age at HBeAg seroconversion was associated with a higher risk for liver cirrhosis ^{[21,23,24,25]}.

Drug therapy was correlated with faster HBeAg seroconversion. A meta-analysis reported increased HBeAg seroconversion rates after one year of NA treatment regardless of NAs efficacy ^[26]. Another meta-analysis reported that the rates of HBeAg seroconversion after one year of treatment are greatest for tenofovir ^[27] and telbivudine ^[28]. However, for longer treatment periods (of 3, 4, 5 years) the rate of HBeAg seroconversion is reduced compared to spontaneous seroconversion obtained without treatment ^[28,29], indicating that a prolonged duration of treatment with, in particular the highly efficient NAs entecavir and tenofovir, is deleterious in achieving HBeAg seroconversion. Given the importance of HBeAg in HBV pathogenesis and the contradictory reports regarding spontaneous and drug-induced HBeAg seroconversion, we propose a mathematical modeling approach for studying the dynamics of HBeAg loss under various hypotheses. The developed models will investigate the role of antibody formation, their role in the disease transition from HBeAg-positive to HBeAg-negative infections, and the trade-off between virus loss during therapy and HBeAg seroclearance.

Over the past decades, mathematical models have been developed to study the dynamics of acute, chronic, and occult HBV infections ^[30,31,32], drug therapy ^{[33,34,35,36,37,38,39,40,41]}, cell-to-cell transmission ^[42], intracellular interactions ^[42,43,44], cellular immune responses ^{[31,34,45,46,47]}, antibody-mediated immune responses ^[44,48,49], HBeAg ^[44,50], and HBeAb ^[44]. We build on the previous modeling work and consider the interaction between HBeAg, cellular immune responses, HBeAb levels and drug efficacy. We hypothesize that B cells mature into HBeAb-producing plasma cells during the HBeAg-positive stages of HBV infection, investigate various modulation mechanisms for HBeAb dynamics, and use the models to predict the differences in seroconversion times under treatment and in the absence of treatment.

This paper is structured as follows. In Section 2, we develop an in-host model of hepatitis B infection in the absence of therapy which focuses on the function of HBeAg in disease progression. In Section 3, we investigate the model analytically and numerically and predict the interplay between cellular and antibody responses on HBeAg seroclearance. In Section 4, we investigate the role of NAs treatment on HBeAg seroclearance. We conclude with a discussion.

2.   Model development

We model the interactions between uninfected hepatocytes, $T$; infected hepatocytes, $I$; hepatitis B virus, $V$; effector cytotoxic T lymphocytes (CTLs), $E$; HBeAg, $e$; HBeAb, $A$; and HBeAg-HBeAb immune complexes, $X$. To incorporate the ability of the liver to regenerate after cell loss ^[51], we assume that uninfected hepatocytes follow a logistic growth with maximum proliferation rate $r$ and carrying capacity $T_m$. They become infected with HBV at rate $\beta$. Infected hepatocytes are killed by effector cells at rate $\mu$, and produce virus and HBeAg at rates $p$ and $\pi$, respectively. Furthermore, we assume that infected cells proliferate with maximum proliferation $r$ ^[32,37,45]. Since there is evidence, at least in acute HBV infections ^[32,45,52], that covalently closed circular DNA can be lost during cell proliferation, we assume that an infected cell produces one uninfected and one infected offspring. HBV is cleared at rate $c$. Effector cells are recruited at constant rate $s_E$ and, after interaction with infected cells, expand at maximum rate $\alpha$ and carrying capacity $E_m$ ^[46]. It has been reported that HBeAg suppresses the cellular immune responses ^[6,53]. We model this by decreasing the effector cell recruitment at rate $\sigma$. Effector cells die at rate $d_E$. We only model the cytolytic effect of the effector cells and ignore their non-cytolytic function ^[31]. HBeAgs decay at per capita rate $\delta_e$ and bind HBeAb to form complexes at rate $k_p$. Complexes dissociate at rate $k_m$ and decay at per capita rate $c_X = \delta_e + d_X$, where $\delta_e$ is the decay rate of HBeAg and $d_X$ is increased removal due to phagocytosis. Previous papers ^[54] have presented detailed models of B lymphocyte proliferation and differentiation into plasmablast, antibody producing plasma cells and memory cells after they encounter antigen. For simplicity, we ignore the details of B-lymphocyte dynamics and differentiation into antibody producing cells and assume that free HBeAb, $A$, is produced at rate $s$ proportional to the HBeAg. Moreover HBeAb is maintained after HBeAg clearance through antigen-independent homeostatic proliferation of memory B cells and long-lived plasma cells, which we model through a logistic term with maximum proliferation rate $s_A$ and carrying capacity $A_m$ ^[49]. The corresponding system of equations is given by system (2.1) and a schematic representation is shown in Figure 1.

All parameters in this model are positive. Moreover the initial conditions are

3.   Results

3.1.   Analytical results

System (2.1) has four virus clearance equilibria. Equilibrium $S_2 = \left(T_m, 0, 0, 0, 0, 0, 0 \right)$ where clearance happens in the absence of immune responses, equilibrium $S_4 = \left(T_m, 0, 0, 0, \frac{E_m(s_E-d_E)}{s_E}, 0, 0 \right)$ where clearance is achieved due to CTL responses, and equilibrium $S_8 = \left(T_m, 0, 0, 0, 0, A_m, 0 \right)$ where HBeAb responses are maximal and CTL responses are absent are unstable. Lastly, equilibrium $S_{10} = \left(T_m, 0, 0, 0, \frac{E_m(s_E-d_E)}{s_E}, A_m, 0 \right)$ where CTL responses are present and HBeAb responses are maximal is locally asymptotically stable when $(T_m \beta p s_E)/(E_m \mu c (s_E-d_E)) < 1$ and $s_E > d_E$. This means that when virus and infected cells are cleared faster than produced, while CTLs are produced faster than cleared, infection dies out.

Similarly, positive virus equilibria can be obtained in the absence of any immune responses ($S_5$), in the presence of anergic CTL responses ($S_6$), under positive but inefficient antibody responses ($S_{11}$), and in the presence of both anergic CTLs and inefficient antibody responses ($S_{12}$). The exact description of these equilibria, other equilibria that are not mentioned here, and their stability can be found in the supplementary material (Appendix 1).

3.2.   Numerical results

3.2.1.   Parameter values

The parameter values used in the simulations are given in Table 1. Liver makes up one fiftieth of a persons weight ^[55]. Each gram of liver contains $14 \times 10^7$ hepatocytes ^[56]. Therefore, an average 70 kg person has about $2 \times 10^{11}$ hepatocytes. As in ^[45], distributing this throughout 15L of extracellular fluid results in a liver capacity of $T_m = 13.6 \times 10^6$ hepatocytes/mL.

Liver cells have the potential to regenerate fast. We account for it by setting the hepatocyte proliferation rate to $r = 1$/day, as in ^[31,49]. The estimates for the infected hepatocytes' half-lives range between 10–100 days ^[30]. Given that we assume a maximum effector cells' level of $E_m = 10^3$ cells/mL, and an infected hepatocytes' half-life of 11 days, the infected cells clearance rate becomes $\mu = \ln{2}/(11\times E_m) = 6 \times 10^{-5}$ mL/(cells$\times$day). The estimates for virus' half-life is at most 4.4 hours ^[57]. We assume a half-life of 4 hours, corresponding to a decay rate of $c = 4.2$/day. The viral infectivity rate during acute HBV infection was estimated to range between $10^{-10}-1.8 \times 10^{-9}$ mL/(virus$\times$day) ^[49] and the virus production rate ranges between 200–1000/day ^[58]. Since virus levels are lower during chronic disease, we consider a one-fold reduction in the infectivity rate $\beta = 4 \times 10^{-11}$ mL/(virus$\times$day) and keep the viral production levels as in acute cases, $p = 400$/day. Similar results can be obtained for high (acute level) infectivity and reduced production rates (not shown). As in ^[50], we assume that HBeAg is degraded at rate $\delta_e = 0.3$/day. The half-life of effector cells is short, therefore we assume a decay rate of $d_E = 0.5$/day ^[47,50,59].

The dissociation rate of HBsAg-HBsAb immune complexes is $k_m = 10$/day ^[49]. We assume the same is true for HBeAg-HBeAb dissociation rate. Furthermore, we assume that HBeAb is mostly IgG and that the avidity for HBeAg-HBeAb binding is similar to that for HBsAg-HBsAb binding, which is set at $K = 10^{7} \text{M}^{-1}$ as in ^[60]. For an IgG molecular weight of $150\text{kDa} = 150 \times 10^3 \text{g} / \text{mol}$ ^[49], and a conversion of $1 \text{ IgG in mg/mL} = 13.43 \text{ IgG in IU/mL}$ ^[61], we obtain a binding rate

The remaining parameters are chosen as follows. We assume that infected cells produce HBeAg at rate $\pi = 10^{-4}$ IU/(cells$\times$day), and activate effector cells at rate $\alpha = 2.1$/(cells$\times$day). Further, effector cells are activated in an infected cell independent manner at rate $s_E = 0.1$/day and their production is inhibited by HBeAg at rate $\sigma = 10^4$ mL/IU. HBeAb is produced in an HBeAg-independent manner at rate $s_A = 10^{-10}$/day and has a carrying capacity of $A_m = 10$ IU/mL. HBeAg-dependent HBeAb production rate $s$ will be varied throughout our investigations. Lastly, the HBeAg-HBeAb immune complex removal rate is $c_X = 1.2$/day, four times higher than the clearance of free HBeAg ^[62].

3.2.2.   Initial conditions

While free HBeAb is not detected by assays during most chronic HBeAg-positive HBV infections, it is reasonable to believe that antibody specific for HBeAg are present in immune complexes before the free antibody can be detected. This has been shown in HIV infections where immune complexes have been detected three weeks prior to free antibody detection ^[63]. However, since the number of immune complexes is small, we model this by assuming HBeAb is initially negligible, i.e. $A = 0$, and $X = 0$. Under this assumption HBeAb does not influence the dynamics of the remaining variables and system (2.1) reduces to

Furthermore, during HBeAg-positive infection, the CTL immune responses are suppressed. We model this by assuming an immune tolerant equilibrium, in which the CTL responses are non-existent or reduced. Asymptotic analysis of model (3.2) (see supplementary material, Appendix 1) shows that there are at most five equilibria in which virus population is non-zero: a no CTL state in which the entire liver is infected

and at most four CTL-inefficient infectious states

where

and $\bar{I}$ is a root of the fourth degree polynomial

The coefficients $C_0, C_1, C_2, C_3$, and $C_4$ are defined in Appendix 1.

Equilibrium $S^{noA}_5$ represents completely absent CTL responses, while $S^{noA}_6$ represents inefficient (exhausted) CTL responses. Studies have shown that during chronic HBeAg-positive infections, cellular immune responses are anergic rather than completely absent ^[64]. We therefore assume that the equilibrium in the absence of HBeAb is given by $S^{noA}_6$, where the CTL responses are non-zero but inefficient. Numerically we find that this is the only stable equilibrium of system (3.2) for the parameter values given in Table 1. Our goal is to investigate how the emergence of antibodies affects HBeAg dynamics and how much this event contributes to the reverting of T cell exhaustion. To address this, we investigate the dynamics of system (2.1) under initial conditions given by the T cell exhaustion state (3.4), together with a small initial free HBeAb concentration $A_0 = 10^{-6}$ IU/mL. This means that at time $t = 0$ we perturb system (2.1) from its unstable equilibrium $(S_6^{noAb}, 0, 0)$ by introducing a small number of free HBeAb. These initial conditions are summarized in Table 2.

3.2.3.   The dynamics of HBeAg and CTL populations

For the parameter values in Table 1, and HBeAg-dependent HBeAb production rate $s\; = \; 6$ mL/(IU$\times$d), the system's dynamics are shown in Figure 2. The solid and dashed curves show the dynamics in the presence and absence of HBeAbs. If the effects of the HBeAb are negligible, then the HBeAg levels do not change and the CTL responses are not strong enough to cause viral remission (see Figure 2, dashed curves). Contrarily, spontaneous HBeAb production leads to virus suppression below the threshold level of HBeAg-negative infections, corresponding to $10^{4}$ cp/mL ^[11], which for the remainder of our study will be called low level virus concentration. The HBeAg population drops below $0.1$ IU/mL, corresponding to the lower limit for HBeAg quantification assays, which for the remainder of our study will be called HBeAg seroclearance level (see Figure 2, solid curves).

We investigated which immune factors are responsible for HBeAg seroclearance, as well as for the infected cells' decay. Both HBeAg and infected cells decrease in biphasic manner (see Figure 3). The first slope decay is steeper for HBeAg than for infected cells, while the second slope decays are the same (see Figure 3, grey versus black lines). The additional removal of free HBeAg during the first phase decay is due to antibody binding, and the formation and removal of immune-complexes by phagocytes. Following the initial antibody responses, a decrease in HBeAg levels leads to a decrease in their tolerogenic pressure on CTLs, which get activated and kill infected cells. As a result, infected cells do get removed by potent cellular immune responses, but HBeAg production slows down as well. During the second phase decay, CTL responses outweigh the antibody responses, hence, the slopes for both HBeAg and infected cells' decay are identical.

Together, this analysis leads to the prediction that a combination of antibody and cellular responses is needed to drive the system into a state of low level virus concentration and undetectable HBeAg. Antibody responses act first, by reducing the immune tolerant effects of HBeAg, while cellular responses control the later stages of the infection.

As shown in Figure 2, under HBeAb responses which grow in an HBeAg-dependent manner at rate $s = 6$ mL/(IU$\times$d), HBeAg is cleared in 2.8 years. We quantified the time of HBeAg seroclearance, $\tau(s)$, as a function of the HBeAb production rate $s$ for $5.51 \leq s \leq 50$ mL/(IU$\times$d) (see Figure 4). For $s < 5.51$ mL/(IU$\times$d), HBeAg seroclearance is never reached. The time to HBeAg seroclearance (below $0.1$ IU/mL) decreases from $41.9$ years for small $s = 5.51$ mL/(IU$\times$d) to $1.2$ years for $s > 50$ mL/(IU$\times$d) (see Figure 4). As expected, higher antibody expansion rates lead to shorter times to HBeAg serocleareance. The time to HBeAg seroclearance is dependent on the initial HBeAb level, with faster HBeAg clearance for high initial HBeAb levels (see Figure S1). In our simulations we fixed $A_0 = 10^{-6}$ IU/mL to avoid any numerical problems.

3.3.   Core and precore mutations

Host-virus predictors for HBeAg seroconversion are limited, especially in chronic infections with high virus loads ^[65]. Two positive events that lead to HBeAg loss are virus mutations in the core promoter ^[66] and precore region ^[67], which prevent the coding of HBeAg and lead to loss of HBeAg-positive virus. Clinical studies have found that at least one core/precore mutation can be found in $89\%$ of HBeAg-negative patients and $56\%$ of HBeAg-positive patients ^[68]. Another study ^[69] has reported even higher numbers of HBeAg-positive patients with core/precore mutations. We, therefore, assume that core/precore mutations precede the HBeAg seroclearance. However, the quantitative relationship between mutation rates and time to HBeAg seroconversion is not known. Core/precore mutations in our investigations can be seen as a proxy for any event leading to an infected cell's loss of HBeAg production. Here, we investigate the role of mutations in HBeAg clearance as follows. We modify system (2.1) to include both wildtype HBeAg-positive virus $V_w$, and mutant HBeAg-negative virus $V_m$. Cells are either infected by wiltdtype or mutant viruses, $I_w$ and $I_m$. At time $t_1$, a fraction $\Phi$ of the viruses produced by $I_w$ are HBeAg-negative and the rest $1-\Phi$ are HBeAg-positive. We ignore back mutations based on the observation that HBeAg-negative virus is mostly inactive. By contrast, all viruses produced by $I_m$ are HBeAg-negative. Wildtype and mutant viruses replicate at rates $p_w$ and $p_m$, and are cleared at rates $c_w$ and $c_m$. Furthermore, hepatocytes get infected by wildtype and mutant viruses at rates $\beta_w$ and $\beta_m$ and are killed by CTLs at rates $\mu_w$ and $\mu_m$. We assume that during division an infected hepatocyte of either wild or mutant-type ($I_w$ or $I_m$) produces one uninfected and one infected offspring of the same hepatocyte type. The system describing these interactions is

where $\Phi(t) = \begin{cases} & 0, \quad \text{ if } t < t_1\\ & \Phi, \quad \text{ if } t \geq t_1 \end{cases}$.

All parameters in this model are positive. Moreover the initial conditions are $T(0) = T_0, I_w(0) = I_0, I_m(0) = 0, V_w(0) = V_0, V_m(0) = 0, e(0) = e_0, E(0) = E_0, A(0) = A_0, X(0) = X_0,$ where $T_0, I_0, V_0, e_0, E_0, A_0, X_0$ are defined as in the case of system (2.1).

3.3.1.   Analytical results

We assume that the wildtype and mutant viruses are identical in everything but their ability to produce HBeAg, i.e. $\beta_w = \beta_m = \beta$, $\mu_w = \mu_m = \mu$, $p_w = p_m = p$, $c_w = c_m = c$, and analyze the long-term behavior of the system. System (3.6) has the following non-infectious equilibria: viral clearance in the absence of antibody responses, $S^{mut}_2 = \left(T_m, 0, 0, 0, 0, 0, 0, 0, 0 \right)$, viral clearance due to CTL responses, $S^{mut}_4 = \left(T_m, 0, 0, 0, 0, 0, \frac{E_m(s_E-d_E)}{s_E}, 0, 0 \right)$, virus clearance under maximal antibody responses, $S^{mut}_6 = \left(T_m, 0, 0, 0, 0, 0, 0, A_m, 0 \right)$ which are all unstable; and viral clearance in the presence of CTL and maximal antibody responses, $S^{mut}_8 = \left(T_m, 0, 0, 0, 0, 0, \frac{E_m(s_E-d_E)}{s_E}, A_m, 0 \right)$, which is locally asymptotically stable if and only if $s_E > d_E$ and $\frac{T_m \beta p s_E}{E_m \mu c (s_E-d_E)} < 1$.

The model has two equilibria in which only the mutant virus persists and the wildtype virus goes extinct: mutant persistence in the presence of anergic CTL responses $S^{mut}_{11}$, which is unstable; and mutant persistence due to a combination of CTL and maximal antibody responses, $S^{mut}_{12}$. Furthermore, the system has two hyperplanes of equilibria, where HBeAg can take on any value, and in which mutant and wildtype virus coexist: coexistence in the absence of any immune responses, $S^{mut}_9$, which is biologically relevant if $e < T_m \pi / \delta_e$; and coexistence under antibody responses, $S^{mut}_{10}$. Coexistence of wildtype and mutant virus results results in infection of the entire liver. For $S^{mut}_{9}$, $S^{mut}_{10}$, and $S^{mut}_{12}$ we did not perform stability analysis (see supplementary material, Appendix 2 for details).

3.3.2.   Numerical results

We are interested in the relationship between the time of HBeAg seroclearance and the time of precore/core mutations. First, we want to see how the timing of HBeAg seroclearance under mutation (with and without concomitant free HBeAb expansion) compares with the timing of HBeAg seroclearance in the absence of mutations and presence of antibodies. We plot the dynamics of the total virus in the presence of antibody as given by the system without mutations (2.1) for $s = 6$ mL/(IU$\times$d) and $A_0 = 10^{-6}$ IU/mL (see Figure 5, solid curves) and the dynamics of the wildtype virus given by the system with mutations (3.6) for $\Phi = 0.12$ and $A_0 = 0$ IU/mL (see Figure 5, dashed curves). In both scenarios, HBeAg clearance takes 2.8 years (see Figure 5, HBeAg ($e$) panel, solid versus dashed curves). Furthermore, if $A_0 = 10^{-6}$ IU/mL at the time of core/precore mutations in system (3.6), HBeAg clearance takes $1.4$ years, twice as fast as in either of the single event scenarios (see Figure 5, dotted curves). Additionally, model (3.6) predicts the asymptotic loss of wildtype virus, cells infected with the wildtype virus, HBeAg, and HBeAg-HBeAb immune complexes. The total virus population is reduced from $\bar{V}_{total} = 863.5$ cp/mL in model (2.1) to $\bar{V}_{total} = \bar{V}_n = 161.6$ cp/mL in model (3.6). This is independent of antibody help (see Figure 6).

Clinical studies predicted that initial core/precore mutations are followed by the appearance of antibodies against the HBeAg ^[70]. We test the effect that the two sequential events, mutations followed by HBeAb expansion, have on HBeAg loss as follows. Let $\tau_1(t_0, \Phi)$ and $\tau_2(t_0, \Phi)$ be the times between spontaneous HBeAb expansion and HBeAg seroclearance when mutations with frequency $\Phi$ occur at $t = 0$ and $t = t_0 > 0$. In both cases $A_0 = 10^{-6}$ IU/mL antibodies are introduced at time $t = t_0$ (see Figure 7). In other words, $\tau_1$ measures the time between HBeAb expansion (at time $t_0$) and HBeAg seroclearance when core/precore mutations start $t_0$ days before HBeAb expansion (i.e. at time $0$), while $\tau_2$ measures the time between HBeAb expansion (at time $t_0$) and HBeAg seroclearance when core/precore mutations and HBeAb expansion start concomitantly at time $t_0$. Since the system is at equilibrium in the absence of mutations (i.e. in the context of $\tau_2$ from time 0 to $t_0$), $\tau_2(0, \Phi) = \tau_2(t_0, \Phi) = \tau_1(0, \Phi)$ for all $t_0$. The difference between the two seroclearance times, $\tau_1(t_0, \Phi) - \tau_2(t_0, \Phi)$, are shown in Figure 8.

We only consider ranges of $t_0$ and $\Phi$ where $\tau_1(t_0, \Phi)$ is positive, based on the assumption that HBeAb expansion happens before HBeAg seroclearance. For any fraction of mutations $\Phi$, we have $\tau_1(0, \Phi) = \tau_2(0, \Phi)$, i.e. $\tau_1(0, \Phi) - \tau_2(0, \Phi) = 0$. The model predicts that for $\Phi > 0.1$, as $t_0$ increases, $\tau_1(t_0, \Phi)$ and consequently $\tau_1(t_0, \Phi)-\tau_2(t_0, \Phi)$ decrease for small and large $t_0$, and increase for intermediate $t_0$ values (see Figure 8, right panel, for $\tau_1(0, 0.5) - \tau_2(0, 0.5) = 0$). Very large $t_0$ values result in HBeAg seroclearance that is exclusively driven by mutations (see Figure 8, white region). For fixed $\Phi < 0.1$, $\tau_1(t_0, \Phi)-\tau_2(t_0, \Phi)$ decreases as $t_0$ increases leading to $\tau_1\sim\tau_2$ for small $\Phi$ and all $t_0$ within $0 \leq t_0 < 40$ months.

Hepatocyte turnover. Changes in HBeAg are associated with the activation of adaptive immune responses, fluctuations in HBV DNA and ALT levels, and minimal to high liver damage. Events such as spontaneous HBeAg seroconversion and emergence of core/precore mutations, which are considered positive events in the natural course of infection of 70–80% of chronically infected patients ^[71], can precede or follow the HBV DNA and ALT dynamics. Spontaneous HBeAg seroclearance followed by the recovery of HBV-specific T-cell functions, however, may result in hepatitis and liver disease ^[72], through T lymphocyte cytotoxic function ^[73], cytokine production, such as IL-6, IL-12 and TNF-$\alpha$ ^[74], and natural killer-induced INF-$\gamma$ production ^[72].

Since the order of sequential events leading to HBeAg loss is unknown (and patient dependent), we use model (3.6) to determine: (ⅰ) the amount of time it takes to reach HBeAg seroclearance following HBeAb expansion; (ⅱ) the amount of monthly liver turnover, which is defined as the average liver loss due to immune response mediated killing each month ^[75] (here, the time period from HBeAb expansion to HBeAg seroclearance), for various fractions of core/precore mutations. For both questions core/precore mutations start at time $t = 0$ and HBeAb matures at time $t = t_0 \geq 0$.

As defined in the previous section, $\tau_1(t_0, \Phi)$ is the time between HBeAb expansion and HBeAg seroclearance (see Figure 7, left panel). We only consider ranges of $t_0$ and $\Phi$ where $\tau_1(t_0, \Phi)$ is positive, based on the assumption that HBeAb expansion happens before HBeAg seroclearance. Under these assumptions, we set

to be the average monthly hepatocyte turnover over the ($t_0$, $t_0+\tau_1(t_0, \Phi)$) time interval. Heat maps for $\tau_1(t_0, \Phi)$ and $HL(t_0, \Phi)$ are presented in Figure 9.

For large core/precore mutations, $\Phi > 0.5$, HBeAg seroclearance happens before antibody formation for large $t_0$ (see Figure 9, left panel, white region), and in the first $12$ months post HBeAb expansion ($t_0+12$ months following the start of mutations) for small $t_0$ (see Figure 9, left panel, blue region). HBeAg loss is (almost) exclusively due to mutations. As $\Phi$ decreases, $\tau_1(t_0, \Phi)$ becomes less sensitive to the delay $t_0$ between the start of mutations and antibody expansion, indicating that antibodies gain more influence on the progression to HBeAg seroclearance. Finally, for small core/precore mutations, $\Phi < 0.01$, the time between HBeAb expansion and HBeAg seroclearance, $\tau_1(t_0, \Phi)$, is almost constant for a fixed $\Phi$ regardless of $t_0$ (in the range $0\leq t_0 \leq 40$ months). This means that the initial HBeAg loss is almost exclusively antibody-driven and only the second phase decay is influenced by mutations. Previously, monthly hepatocyte turnover was estimated during acute HBV infections to range between $0.12-1$ $T_m$/month ^[32]. Our model predicts a monthly hepatocyte turnover of $0.4$ $T_m$/month (see Figure 9, right panel, red region) for large $\Phi$ and small $t_0$ and $0.1$ $T_m$/month for $\Phi < 0.01$, regardless of $t_0$ (in the range $0\leq t_0 \leq 40$ months).

We next investigate how these results are altered in the presence of antiviral therapy.

4.   HBeAg dynamics during antiviral therapy

4.1.   Treatment in the absence of mutations

We investigate how the dynamics presented in Figure 2 change when we consider the effects of nucleos(t)ide analogous antiviral treatment. Unless noted otherwise, we assume that treatment and HBeAb expansion occur concomitantly at time $t = 0$. NAs suppress viral replication^[9,10]. We model this effect by reducing the virus production rate $p$ by the treatment efficacy $0.5 < \epsilon < 1$. Hence, $p_{treat} = (1-\epsilon)p$.

Under high treatment efficacy $\epsilon = 0.9$, virus undergoes a triphasic decay to a new equilibrium of 2 cp/mL: an instantaneous $0.94\log_{10}$ drop in the first day, followed by a slower $0.5\log_{10}$ decay over the next $4$ months and finally a $6.4\log_{10}$ decay over the next $9.6$ months. The HBeAg does not decay during the first day, and follows a triphasic decay that reaches the limit of detection within $7.9$ months following the start of therapy (see Figure 10, dashed curve). The HBeAg first phase decay is due to treatment induced virus loss, while the second and third decay phases present qualitative dynamics that are similar to the ones obtained in the absence of treatment: where HBeAb efficiently binds HBeAg, leading to the removal of the tolerogenic effect on the CTLs, which eventually control the HBeAg (infected cells) loss (see Figure 10, dashed vs. solid curves). During the third phase, the decay rates of virus and HBeAg levels are steeper during treatment due to CTL-mediated removal of lower virus (infected cell) populations. Treatment efficacies (in the range considered, $0.5 < \epsilon < 1$) influence the time of HBeAg seroclearance, with fast HBeAg seroclearance for high $\epsilon$ (see Figure 11). In the absence of treatment, HBeAg has a seroclearance time of $2.83$ years. By contrast, HBeAg seroclearance time decreases to $11.8$ and $7.2$ months for $\epsilon = 0.5$ and $\epsilon = 0.99$ (see Figure 11).

4.2.   Treatment in the presence of mutations

We next examine how treatment with various drug efficacies $\epsilon$ influences: (ⅰ) the amount of time it takes to reach HBeAg seroclearance after HBeAb expansion; and (ⅱ) the monthly amount of liver turnover between HBeAb expansion and HBeAg seroclearance $HL$, under various fractions of core/precore mutations (see Figure 12). For both questions we assume that mutations start at time $t = 0$, and treatment and HBeAb expansion occur concomitantly at time $t = t_0\in [0, 40]$ months. As before, we define $\tau_1(t_0, \Phi)$ to be the time between HBeAb expansion and HBeAg seroclearance.

As in the absence of treatment, for large $\Phi > 0.5$ and large $t_0$, HBeAg seroclearance happens before antibody formation for all $\epsilon$ (see Figure 12, for $\epsilon = 0.7$). HBeAg loss is exclusively due to mutations. As $\Phi$ decreases, $\tau_1(t_0, \Phi)$ becomes less sensitive to the delay $t_0$ between the start of mutations and antibody expansion, indicating that antibodies gain more influence on the progression to HBeAg seroclearance. Finally, for small core/precore mutations, $\Phi < 0.01$, the time between HBeAb expansion and HBeAg seroclearance, $\tau_1(t_0, \Phi)$, is almost constant regardless of $t_0$ (in the range $0\leq t_0 \leq 40$ months). This means that the initial HBeAg loss is almost exclusively antibody-driven and only the second phase decay is influenced by mutations. After HBeAb expansion, HBeAg is cleared faster under treatment with $\epsilon = 0.7$ than in the absence of treatment (compare Figures 12 and 9, left panels). Furthermore, our model predicts that for $\epsilon = 0.7$ the monthly hepatocyte turnover is greatest for large $\Phi$ and small $t_0$, up to $0.18$ $T_m$/month (see 12, right panel, red region), lower than the $0.4$ $T_m$/month in the absence of treatment. As $\Phi$ decreases, and $\tau_1(t_0, \Phi)$ increases, the monthly hepatocyte turnover decreases to $0.12$ $T_m$/month for $\Phi < 0.01$ regardless of $t_0$ (in the range $0\leq t_0 \leq 40$ months), higher than the $0.1$ $T_m$/month observed in the absence of treatment. We find that low drug efficacy yields an increase in monthly hepatocyte turnover for $\Phi < 0.01$. For example, for $\epsilon = 0.5$, the monthly hepatocyte turnover is $0.16$ $T_m$/month, while for $\epsilon = 0.99$, the average monthly hepatocyte turnover is $0.1$ $T_m$/month, i.e. equal to the average monthly hepatocyte turnover in the absence of treatment. However, since HBeAg seroclearance is faster under treatment, the total hepatocyte turnover is significantly reduced.

5.   Discussion and conclusion

We developed a mathematical model describing the host-pathogen interactions during HBeAg-positive chronic hepatitis B virus infections, with a focus on the effects of HBeAb expansion on disease progression and in particular HBeAg seroclearance. HBeAg, a secretory protein described as an immune tolerogen which downregulates the cellular immune responses ^[6], is removed when HBeAb is produced. The moment at which HBeAb production starts, its mechanistic interactions with HBeAg and other immune responses, or the timing of HBeAg seroclearance are not well understood. We employed a mathematical modeling approach to provide insight into antibody formation and HBeAg seroclearance. Previous mathematical models have investigated the interplay between hepatocytes, infected hepatocytes, HBV virus ^[30], CTL-mediated immune responses ^{[31,34,45,46,47]}, humoral-mediated immune responses ^[48,49], a class of subviral particles ^[49], and HBeAg ^[50]. To our knowledge, however, only one paper, in the context of occult HBV infection, has incorporated the dynamics of HBeAb ^[44].

We used the model to determine the interplay between HBeAb and CTL immune responses on HBeAg seroclearance in the absence and presence of treatment. Our model predicted that, in the absence of treatment, a combination of CTL and antibody responses are needed to achieve HBeAg seroclearance and viral remission. We found that HBeAg seroclearance follows several key stages: (1) newly matured HBeAb expand, bind HBeAg to form HBeAg-HBeAb immune complexes, which are removed via phagocytosis; (2) HBeAg population decays to levels that no longer affect CTLs activity, leading to CTLs activation and expansion; (3) infected cells are killed by CTLs; and, consequently, (4) less HBeAgs are produced and eventually decay below $0.1$ IU/mL (called seroclearance level). We investigated the relationship between HBeAg seroclearance and the rate of HBeAb production and found that bigger HBeAb production rates lead to shorter times to HBeAg seroclearance.

Host-virus predictors for HBeAg seroclearance are limited. Mutations of HBV in its core or precore region result in loss of expression of HBeAg by hepatocytes infected with the mutated virus ^[66,67]. We assumed that these mutations are followed by antibody formation and eventual HBeAg seroclearance, and investigated how core/precore mutations influence the time between HBeAb expansion and HBeAg seroclearance and the hepatocyte turnover. We found that large fractions of mutations result in fast HBeAg seroclearance without the help of HBeAbs and high monthly hepatocyte turnover. Intermediate fractions of mutations result in longer seroclearance times and lower monthly hepatocyte turnover. For small fractions of mutations HBeAg seroclearance is driven almost exclusively by HBeAbs, the seroclearance times are longest and monthly hepatocyte turnover is lowest. Furthermore, we found that mutations can clear HBeAg even in the absence of antibodies.

Increased HBeAg seroconversion rates after one year of NA treatment, regardless of NAs efficacy, have been observed in clinical studies ^[26]. We predicted that, regardless of the treatment efficacy, the time from HBeAb expansion to HBeAg seroclearance is significantly shortened by treatment, ranging between $7$ month for high efficacy treatment and $11$ months for low efficacy treatment, i.e. within one year, compared to $2.8$ years in the absence of treatment. Our model does not, however, inform why a prolonged duration (several years) of high efficacy treatment is deleterious to the achievement of HBeAg seroconversion ^[28].

Under treatment, as in the absence of treatment, mutations influence the length of time between HBeAb expansion and HBeAg seroclearance as well as the amount of monthly hepatocyte turnover. The qualitative results are the same as in the untreated case, with fastest HBeAg seroclearance and highest monthly hepatocyte turnover for large fractions and slowest HBeAg seroclearance and lowest monthly hepatocyte turnover for small fractions of mutations. While the time to HBeAg seroclearance under treatment is always shorter than in the absence of treatment, the monthly hepatocyte turnover is decreased only for large fractions. Small fractions of mutations yield an increased monthly hepatocyte turnover under treatment.

When we modeled core/precore mutations, we assumed similar characteristics for the wildtype and mutant virus strains and their corresponding infected cell classes. Studies have reported different half-lives for the wildtype and mutant virus strains, although the exact values vary. Dandri et al. reported average half-lives of $46$ and $2.5$ minutes ^[76]. Ribeiro et al. reported half-lives of $25.2$ and $13.1$ hours and shorter half-lives of cells infected by mutant compared to wildtype virus, $12.1$ days versus $16$ days ^[77]. Lastly, it has been suggested that cells that express mutant virus can lose the ability to express wildtype virus ^[78]. This can be modeled by assuming that a fraction $z$ of cells infected with the wildtype virus transitions into the population of cells infected with the mutant virus. We have investigated how differences in the two strains affect our results by considering a $10$-fold increase in the mutant clearance rate, combined with a $1.3$-fold increase in the clearance rate of cells infected with the mutant virus, and a $z = 0.01$ transitioning rate from wildtype to mutant infection. We found changes in the timing of HBeAg seroclearance, with the longest time occurring when mutant clearance rates are increased and the transition is zero, and the shortest time occurring when the mutant clearance rates are kept at the wild-type levels and the transition is non-negative (see Appendix 3, heterogeneity in infected and virus populations section, Figures S15 and S16).

Our model assumes that more than $75\%$ of the liver is infected during chronic HBV infections. While such large values have been reported in acute infections ^[31,79], it is possible that fewer hepatocytes get infected during chronic disease ^[80]. We performed a sensitivity analysis on the infectivity parameter $\beta$ which showed that our results are sensitive to changes in $\beta$, and, consequently, to having a lower amount of the liver infection in the first 2 years (Appendix 3, sensitivity analysis section, Figure S4). Further work is needed to determine how the results will change if we assume lower level of liver infection. There are many other parameters for which we have limited quantitative information. Most of these parameters were kept fixed in our model. The model's sensitivities to a number of these parameters are shown in Figures S2–S11 (Appendix 3, sensitivity analysis section). Another limitation of our model is that the production rate of virus is unaffected by the non-cytolytic effects of the effector cell population. Effector cells produce cytokines such as IFN-$\gamma$ that inhibit the HBV replication and thus, the HBV virions production rate. Further work is needed to determine how the results will change if we assume that effector cells inhibit the virus production.

In conclusion, we have built a mathematical model of HBeAg seroclearance and investigated how HBeAb and CTL work together to secure the transition into an HBeAg-negative HBV infection and how the efficacy of drug therapy affects the timing to HBeAg seroclearance. Such results are important for understanding this milestone event in HBV natural history and can be used to inform human interventions.

Acknowledgments

SMC and SK acknowledge funding from National Science Foundation grant No. 1813011. We thank the reviewers for the valuable comments and suggestions.

Conflict of interest

The authors declare there is no conflict of interest.

Supplementary

Appendix 1. Stability analysis of the system without HBeAb. System (3.2) has the following non-negative equilibria: no liver, no CTL responses

infection free state

no liver under CTL responses

clearance due to CTLs

infection in the absence of CTL responses

and infection under immune anergy. There are up to four solutions of this form

where

and $\bar{I}$ is a root of the fourth degree polynomial

with

Equilibria $S^{no A}_1$ and $S^{no A}_3$ are solutions in which the total hepatocyte population is zero, corresponding to the death of the patient.

Theorem 5.1. Equilibria $S^{no A}_1$ and $S^{no A}_3$ are unstable.

Proof. The Jacobian of system (3.2) evaluated at $S^{no A}_1$ or $S^{no A}_3$ has eigenvalue $\lambda_1 = r > 0$. Hence $S^{no A}_1$ and $S^{no A}_3$ are unstable equilibria.

Equilibria $S^{no A}_2$ and $S^{no A}_4$ are solutions in which the infection is cleared. In $S^{no A}_4$ a population of memory cells persists over time, while in $S^{no A}_2$ CTLs go extinct.

Theorem 5.2. Equilibrium $S^{no A}_2$ is unstable.

Proof. Using Maple we find that the Jacobian of system (3.2) evaluated at $S^{no A}_2$ has the eigenvalue

which is always positive. Hence $S^{no A}_2$ is an unstable equilibrium.

Theorem 5.3. Equilibrium $S^{no A}_4$ exists and is locally asymptotically stable iff $s_E > d_E$ and $\frac{T_m \beta p s_E }{E_m c \mu (s_E-d_E)} < 1$.

Proof. Using Maple we find that the Jacobian of system (3.2) evaluated at $S^{no A}_4$ has two eigenvalues $\lambda_1 = -r$ and $\lambda_2 = -\delta_e$ which are always negative. A third eigenvalue $\lambda_3 = d_E - s_E$ is negative iff $s_E > d_E$. The remaining two eigenvalues are

If $s_E > d_E$, then $(c s_E +E_m \mu (s_e-d_E)) > 0$ and hence $Re(\lambda_5) < 0$. For $\lambda_4$ we find

In equilibria $S^{no A}_5$ and $S^{no A}_6$ virus is not cleared. $S^{no A}_5$ represents a state in which CTLs have vanished and the entire liver is infected, while $S^{no A}_6$ is a state in which CTLs are ineffective.

Theorem 5.4. Equilibrium $S^{no A}_5$ is locally asymptotically stable iff $\frac{\delta_e(T_m \alpha + s_E)}{d_E(T_m \pi \sigma + \delta_e)} < 1$.

Proof. Using Maple we find that the Jacobian of system (3.2) evaluated at $S^{no A}_5$ has four eigenvalues $\lambda_1 = -r$, $\lambda_2 = -\delta_e$, $\lambda_3 = - c$, and $\lambda_4 = -\frac{T_m \beta p}{c}$ which are always negative. The fifth eigenvalue

is negative iff $\frac{\delta_e(T_m \alpha + s_E)}{d_E(T_m \pi \sigma +\delta_e)} < 1$.

We did not attempt to analytically study the asymptotic stability of $S^{no A}_6$.

Stability analysis of the system with HBeAb. System (2.1) has the following non-negative equilibria: no liver, no CTLs

infection free state

no liver under CTL responses

clearance due to CTL responses

infection in the absence of any immune responses

infection during anergic CTL responses. There are up to four solutions of this form

$S_1$ through $S_6$ are the equilibria of (3.2) with additional variables $\bar{A} = 0$ and $\bar{X} = 0$.

Theorem 5.5. Equilibria $S_1$ through $S_4$ are unstable.

Proof. The Jacobian of system (2.1) evaluated at $S_1$ through $S_4$ has eigenvalue $\lambda_1 = s_A > 0$. Hence $S_1$ through $S_4$ are unstable equilibria.

Theorem 5.6. Equilibrium $S_5$ is locally asymptotically stable iff $\frac{\delta_e(T_m \alpha + s_E)}{d_E(T_m \pi \sigma + \delta_e)} < 1$ and $\frac{(c_X+k_m)(T_m s \pi + \delta_e s_A)}{T_m k_p c_X \pi} < 1$.

Proof. Using Maple we find that the Jacobian of system (3.2) evaluated at $S^{no A}_5$ has four eigenvalues $\lambda_1 = -r$, $\lambda_2 = -\delta_e$, $\lambda_3 = - c$, and $\lambda_4 = -\frac{T_m \beta p}{c}$ which are always negative. The remaining three eigenvalues are

which is negative iff $\frac{\delta_e(T_m \alpha + s_E)}{d_E(T_m \pi \sigma +\delta_e)} < 1$, and

We find $Re(\lambda_{6, 7}) < 0$ iff

and

We find that (5.14) implies (5.13), and hence all eigenvalues have negative real part and $S_5$ is locally asymptotically stable iff $\frac{\delta_e(T_m \alpha + s_E)}{d_E(T_m \pi \sigma + \delta_e)} < 1$ and $\frac{(c_X+k_m)(T_m s \pi + \delta_e s_A)}{T_m k_p c_X \pi} < 1$.

Theorem 5.7. If $S_6^{no A}$ is unstable in system (3.2), then $S_6$ is unstable in system (2.1).

Proof. For $S_6$ the equilibrium levels of $A$ and $X$ are zero. When plugging $A = 0$ and $X = 0$ into the Jacobian, $J$, of system (2.1) yields

where $J^{no A}$ is the Jacobian of system (3.2). Hence, if $J^{no A}$ has an eigenvalue with positive real part, so does $J$.

The remaining equilibria have an antibody response component. They are: no liver and maximal antibody responses

infection free due to maximal antibody responses

no liver under anergic CTLs and maximal antibody responses

infection free due to combined CTLs and antibody responses

Theorem 5.8. Equilibria $S_7$ and $S_9$ are unstable.

Proof. The Jacobian of system (2.1) evaluated at $S_7$ or $S_9$ has eigenvalue $\lambda_1 = r > 0$. Hence $S_7$ and $S_9$ are unstable equilibria.

Theorem 5.9. Equilibrium $S_8$ is unstable.

Proof. Using Maple we find that the Jacobian of system (2.1) evaluated at $S_8$ has eigenvalue

which is always positive. Hence $S_8$ is an unstable equilibrium.

Theorem 5.10. Equilibrium $S_{10}$ is locally asymptotically stable iff $s_E > d_E$ and $\frac{T_m \beta p s_E}{E_m \mu c (s_E-d_E)}$.

Proof. Using Maple we find that the Jacobian of system (2.1) evaluated at $S_{10}$ has two eigenvalues $\lambda_1 = -r$ and $\lambda_2 = -s_A$ which are always negative. A third eigenvalue $\lambda_3 = d_E - s_E$ is negative iff $s_E > d_E$. The remaining four eigenvalues are

and

If $s_E > d_E$, then $Re(\lambda_5), \ Re(\lambda_6), \ Re(\lambda_7) < 0$, and

Infection with no CTLs and positive antibody responses includes up to two states, given by

where $\bar{e}$ satisfies the quadratic equation $C_2 \bar{e}^2 + C_1 \bar{e} +C_0 = 0$ with

We did not attempt to analytically study the asymptotic stability of $S_{11}$

Lastly, infection during anergic CTLs and inefficient antibody responses is given by:

where

and $\bar{y}$ is the root of a degree seven polynomial whose expression is messy and will not be provided here (Maple file available upon request). We did not attempt to analytically study the asymptotic stability of $S_{12}$.

Appendix 2. Stability analysis of the system with mutations. We determine system (3.6)'s equilibria and their stability under the assumption that $\beta_{m} = \beta_{w} = \beta$, $\mu_m = \mu_w = \mu$, $p_m = p_w = p$, and $c_m = c_w = c$. System (3.6) has the following non-negative equilibria.

where $e$ can take on any value less than or equal to $\frac{T_m \pi}{\delta_e}$ (for biological relevance, i.e. non-negativity).

where $e$ can take on any value, and

where $\bar{I}_n$ is a root of the fourth degree polynomial

with

where $\bar{I}_n$ is defined as in $S^{mut}_{11}$.

Theorem 5.11. Equilibria $S^{mut}_1-S^{mut}_7$, and $S^{mut}_{11}$ are unstable.

Proof. The Jacobian of system (3.6) evaluated at $S^{mut}_1, \ S^{mut}_3, \ S^{mut}_5$, and $S^{mut}_7$ has eigenvalue $\lambda_1 = r > 0$. The Jacobian of system (3.6) evaluated at $S^{mut}_2$, $S^{mut}_4$, and $S^{mut}_{11}$ has eigenvalue $\lambda_1 = s_A > 0$. The Jacobian of system (3.6) evaluated at $S^{mut}_6$ has eigenvalue $\lambda_1 = -\frac{c}{2} + \frac{\sqrt{4 T_m \beta p + c^2}}{2} > 0$. Thus, $S^{mut}_1-S^{mut}_7$, and $S^{mut}_{11}$ are unstable equilibria.

Theorem 5.12. Equilibrium $S^{mut}_8$ is locally asymptotically stable iff $s_E > d_E$ and $\frac{T_m \beta p s_E}{E_m \mu c (s_E-d_E)} < 1$.

Proof. Using Maple we find that the Jacobian of system (3.6) evaluated at $S^{mut}_8$ has two eigenvalues $\lambda_1 = -r$ and $\lambda_2 = -s_A$ which are always negative. A third eigenvalue $\lambda_3 = d_E-s_E$ is negative iff $s_E > d_E$. The remaining six eigenvalues are

and

If $\lambda_3 < 0$, i.e. $s_E > d_E$, then $-(c s_E +E_m \mu (s_e-d_E)) < 0$ and hence $Re(\lambda_5), Re(\lambda_7) < 0$. Further $Re(\lambda_8), Re(\lambda_9) < 0$, and

Since $\Phi \geq 0$, if $Re(\lambda_4) < 0$ then $Re(\lambda_6) < 0$. Hence all eigenvalues have negative real part and $S^{mut}_8$ is locally asymptotically stable iff $s_E > d_E$ and $\frac{T_m \beta p s_E}{E_m \mu c (s_E-d_E)} < 1$.

Proving stability for $S^{mut}_{9}, S^{mut}_{10}$, and $S^{mut}_{12}$ is messy and will not be presented here.

Appendix 3. Sensitivity Analysis. We analyzed the time-dependent semi-relative sensitivity of model (2.1)'s dynamics to small changes in various parameters. For times $0$ to $2$ years and parameters $P = \{\alpha, A_m, \beta, k_p, \mu, \pi, s, s_A, s_E, \sigma\}$ we look at the semi-relative sensitivity curves $P \frac{\partial Y}{\partial P}$, where $Y = \{T, I, V, e, E, A, X\}$. The results are presented below. As expected, we find that $\{\alpha, A_m, \mu, s, s_A, s_E \}$ have negative effects on $\{I, V, e\}$ and positive effects on $E$ (see Figures S2, S3, S6, S8, S9, S10), while $\{\beta, k_p, \pi, \sigma\}$ have the opposite effects (see Figures S4, S5, S7, S11). Furthermore, we observe that virus, HBeAg, and effector cell populations are most sensitive to changes in $k_p$ and $s$ (see Figures S12–S14, left panel). While HBeAg and effector cells are most sensitive around day $125$ regardless of the parameter, the maximal sensitivity of the virus population to $\beta, A_m$ and $\mu$ occurs later (see Figures S12–S14, right panel).

Heterogeneity in infected and virus populations. Studies have reported different half-lives for the wildtype and mutant virus strains, although the exact values vary. Dandri et al. reported average half-lives of $46$ and $2.5$ minutes, corresponding to $20$-fold higher clearance rate of mutant compared to wildt-ype strains ^[76]. Ribeiro et al. reported half-lives of $25.2$ and $13.1$ hours, corresponding to $2$-fold increase in clearance of mutant compared to wildtype virus ^[77]. Moreover, Ribeiro et al., reported shorter half-lives of cells infected by mutant compared to wildtype virus, $12.1$ days versus $16$ days (a $1.3$-fold increase in the clearance rate of $I_m$ compared with $I_w$) ^[77]. In this study we assumed that both viruses have half-lives of $4$ hours, corresponding to the clearance rates $c_w = c_m = 4.2$/day and both infected cell types are removed (by the effector cells) at rate $\mu$. To address how differences in the two strains half-lives affect our results, we change their clearance rates as follows: we keep $c_w = 4.2$/day, and assume a $10$-fold increase in the mutant clearance, $c_m = 42$/day, corresponding to $24$ minutes half-life. Moreover, we keep the killing rate of the cells infected with the wildtype virus at $\mu_w = \mu = 6 \times 10^{-5}$ and increase $\mu_m = 1.3\times \mu = 8 \times 10^{-5}$, corresponding to half-lives of $11.5$ and $8.7$ days under maximal CTL response $E_m$. Lastly, it has been suggested that cells that express mutant virus can lose the ability to express wildtype virus ^[78]. This can be modeled by assuming that a fraction $z$ of cells infected with the wildtype virus transitions into the population of cells infected with the mutant virus. The system of equations expressing these changes is

We solved the model for four parameter combinations: (ⅰ) $\Phi = 0.1$, $z = 0$, $c_w = c_m$, $\mu_w = \mu_m$; (ⅱ) $\Phi = 0.1$, $z = 0$, $c_w < c_m$, $\mu_w < \mu_m$; (ⅲ) $\Phi = 0.1$, $z\neq 0$, $c_w = c_m$, $\mu_w = \mu_m$; (ⅳ) $\Phi = 0.1$, $z\neq0$, $c_w < c_m$, $\mu_w < \mu_m$. We predict the longest time to HBeAg seroclearance in case (ⅱ), when cells infected by the wildtype keep their ability to produce wildtype virus ($z = 0$), and the mutant families $V_m$ and $I_m$ have faster clearance rates ($10$ and $1.3$-fold increase, respectively) (see Figure S15, solid lines). When a fraction $z = 0.01$ of cells infected with the wildtype transition into cells producing only mutant virus (cases ⅲ and ⅳ), we obtain faster HBeAg seroclearance times (see Figure S15, dashed and dotted lines). As before, however, increasing the mutant virus and cells infected with mutant virus clearance rates (case ⅳ), decreases the overall HBeAg seroclearance time (see Figure S15, dotted lines). Interestingly, the dynamics of the model now allow for long-term coexistence between the mutant and wildtype virus strains.

We next investigate how these results change if we vary $\Phi$ and $z$. We let $\tau_{same}(z, \Phi)$ be the time to HBeAg seroclearance assuming equal clearance rates for wildtype and mutant populations, $\tau_{diff}(z, \Phi)$ be the time to HBeAg seroclearance assuming a $10$ and $1.3$ fold increased clearance rate of mutant virus and cells infected with the mutant virus, and

be the average monthly hepatocyte turnover. Here $\tau$ is either $\tau_{same}$ or $\tau_{diff}$. Note, that we assume that mutations, transitions, and HBeAb formation all start concomitantly at time $t = 0$. We find that for all $10^{-4} < z < 10^{-1}$ and $10^{-2} < \Phi < 1$, the $10$ and $1.3$-fold increases in clearance rates of mutant virus and cells infected with mutant virus result in slower HBeAg clearance ($\tau_{same}(z, \Phi) - \tau_{diff}(z, \Phi) < 0$ in Figure S16, left panel), but decreased monthly hepatocyte turnover ($HL_{same}(z, \Phi) - HL_{diff}(z, \Phi) > 0$, in Figure S16, right panel). Furthermore, the negative impact that increased mutant clearance rates have on the time to HBeAg seroclearance is most prevalent for intermediate fractions of mutations $\Phi$ and intermediate transition rates $z$ (see Figure S15, left panel, blue region).