Understanding the antiviral effects of RNAi-based therapy in HBeAg-positive chronic hepatitis B infection

The RNA interference (RNAi) drug ARC-520 was shown to be effective in reducing serum hepatitis B virus (HBV) DNA, hepatitis B e antigen (HBeAg) and hepatitis B surface antigen (HBsAg) in HBeAg-positive patients treated with a single dose of ARC-520 and daily nucleosidic analogue (entecavir). To provide insights into HBV dynamics under ARC-520 treatment and its efficacy in blocking HBV DNA, HBsAg, and HBeAg production we developed a multi-compartmental pharmacokinetic–pharamacodynamic model and calibrated it with frequent measured HBV kinetic data. We showed that the time-dependent single dose ARC-520 efficacies in blocking HBsAg and HBeAg are more than 96% effective around day 1, and slowly wane to 50% in 1–4 months. The combined single dose ARC-520 and entecavir effect on HBV DNA was constant over time, with efficacy of more than 99.8%. The observed continuous HBV DNA decline is entecavir mediated, the strong but transient HBsAg and HBeAg decays are ARC-520 mediated. The modeling framework may help assess ongoing RNAi drug development for hepatitis B virus infection.

www.nature.com/scientificreports/ To better understand the effect of RNAi therapies, additional information regarding the host-virus-drug dynamics and therapy outcomes are needed. In this study, we developed mathematical models that best reproduce observed HBV DNA, HBsAg and HBeAg kinetics following a single dose of ARC-520 in five HBeAg-positive patients from the Heparc-2001 study. Mathematical models of hepatitis B infection have been used to study the dynamics of acute, chronic, and occult HBV infections [25][26][27][28][29] , anti-HBV therapy 14,[30][31][32][33][34][35] , cell-to-cell transmission 36 , intracellular interactions [36][37][38] , cellular immune responses 26,30,[39][40][41] , antibody-mediated immune responses 11,38,42 , HBeAg 38,43,44 , and HBeAb 38 dynamics. We build on previous modeling work, consider the interaction between HBV DNA, HBsAg and HBeAg titers in the presence of a single dose RNAi-based therapy, and use the model to run in silico experiments to predict individual contributions of different drug effects on the dynamics for HBsAg titers.
We develop a mathematical model that considers the interactions between infected hepatocytes, I (in cells per ml); total intracellular HBV DNA, D (in copies per ml); serum HBV DNA, V (in IU per ml); serum HBsAg, S (in IU per ml); and serum HBeAg, E (in IU per ml). We assume that infected cells decay at per capita rate δ , and we exclude cell proliferation (we will relax this assumption later on). We assume intracellular HBV DNA is synthesized at rate α and is lost at constant per capita rate c D . The replication rate α summarizes various steps that are not modeled explicitly, such as the transcription of pregenomic RNA (pgRNA) from cccDNA, and the generation of single stranded DNA by reverse transcription. Intracellular HBV DNA is assembled and released into blood as free virions at rate p which are cleared at rate c. To account for the different units of intracellular and serum virus, we use the conversion factor ξ = 1/5.3 IU/copies 45 . Lastly, we assume HBsAg and HBeAg are transcribed from cccDNA inside infected hepatocytes and then released into blood at rates p S and p E , respectively, and are cleared at per capita rates d S and d E , respectively. The model is given by the following model: Patients were administered daily nucleoside analogous treatment with entecavir starting at day t 0 = 0 . ETV is known to block reverse transcription of HBV DNA, and therefore inhibit HBV DNA synthesis. We model this (see model (5)) as a constant reduction of the HBV DNA synthesis rate α to (1 − ǫ)α , where 0 ≤ ǫ ≤ 1 is the ETV efficacy. Experimental studies in humanized mice have shown that serum HBV DNA declines in biphasic manner while HBV-infected cell are not lost in the first months following NA treatment initiation 46,47 . To account for the biphasic HBV DNA decay in the absence of infected cell killing, we assume that ETV has additional time-dependent inhibitory effects on intracellular HBV DNA synthesis and model it by decreasing α further to α ETV treat = αe −gt (1 − ǫ) , where g ≥ 0 is a constant and t is the time in days post ETV initiation. Moreover, a single ARC-520 dose was administrated at time t 0 = 0 . Unlike ETV, which was given daily, we model the build-up and clearance of ARC-520 pharmacokinetics over time by considering a two-compartment pharmacokinetic model consisting of drug quantity in the plasma and liver, C p and C e , respectively 48 . The inoculum C p (0) = C 0 decays exponentially at rate d = d + k eo , where d is the plasma drug degradation rate and k eo is the absorption into the liver rate. The drug in the liver decays at rate k eo , identical with the absorption rate 49 . Following these assumptions, the pharmacokinetic model has the form: with initial conditions C p (0) = C 0 and C e (0) = 0 . This is a linear model which can be solved to give solutions: www.nature.com/scientificreports/ Lastly, we assume the relationship between the drug quantity in the liver C e (t) and drug efficacy η i (t) to be given by: where η max = 1 is the maximum drug efficacy, EC 50,i are drug quantities that yield half-maximal effects, and i = {1, 2, 3} are the infectious events that are affected by ARC-520 therapy, i.e., the transcription of HBV DNA, the transcription of HBsAg, and the transcription of HBeAg, respectively. The effects of ARC-520 on intracellular HBV DNA, HBsAg and HBeAg are modeled as the reduction of intracellular HBV DNA synthesis α to α ARC treat = (1 − η 1 )α , HBsAg production from p S to p S,treat = (1 − η 2 )p S , and of HBeAg production from p E to p E,treat = (1 − η 3 )p E , respectively. Considered together, models (1) and (4) give the following pharmacokinetics-pharamcodynamics (PK/PD) model: Data fitting. We used published kinetic HBV DNA, HBsAg, HBeAg data in serum measured from five HBeAg-positive, treatment-naive chronic hepatitis B patients as described in the 'Patient data' section.
Parameter values. We assume that, prior to therapy initiation, model (5) describes a persistent chronic infection and is at the quasi-equilibrium, given by the initial values I(0) = I 0 , D(0) = D 0 , V (0) = V 0 , S(0) = S 0 and E(0) = E 0 . Initial values for HBV DNA, V (0) = V 0 ; HBsAg, S(0) = S 0 ; and HBeAg, E(0) = E 0 , are set to the patient data prior to the start of therapy, t −1 = −8 , (day eight prior to the ARC-520 injection). The percentage of HBV-infected hepatocytes is reported to vary between 18 ± 12% in chronic HBsAg carriers 50,51 and 99% in acute infections 26,52 . Without loss of generality, we arbitrary assume that 50% of hepatocytes are infected at the beginning of treatment. Liver contains approximately 2 × 10 11 hepatocytes, which, when distributed throughout 15 liters of extracellular fluid, gives a total hepatocyte concentration T max = 1.4 × 10 7 cells/ml 53 . We set the initial infected hepatocyte population to I 0 = 0.5T max . Lastly, the pre-treatment level of intracellular HBV DNA in HBeAg positive patients is set to D 0 = 225/(I 0 /T max ) = 450 copies/ infected cell, as in 54 .
Since we assume that model (5) is in chronic equilibrium (for the additional assumption δ = 0 ) before the therapy initiation, parameters α , p, p S , p E are fixed according to the following formulas: We start by ignoring the dynamics of infected cells, such as infection of susceptible cells and/or infected cell proliferation (we will relax this assumption in later sections), and assume that infected cells decay due to natural death and immune mediated killing at per capita rate δ = 4 × 10 −3 per day, corresponding to a life-span of 250 days (we will later investigate the effect of increasing the killing rate, to include increased immune mediated killing or RNAi induced toxicity and death). The estimated half-life of intracellular HBV DNA is 24 h 55 , which corresponds to the intracellular HBV DNA decay rate c D = 0.69 per day. ARC-520's half-life has been reported to range between 3 and 5 h 56 , corresponding to decay rates 3.3 < d < 5.5 per day; we fix d = 4 per day. Lastly, we set the initial ARC-520 quantity to the trial dose of C 0 = 4 mg/kg.
The unknown parameters are parm = {g, c, d S , d E , ǫ T , EC 2 , EC 3 , k eo } . Here, (1 − ǫ T ) = (1 − ǫ)(1 − η 1 (t)) accounts for the total drug effect on HBV DNA production. Since preliminary simulations (not shown) indicate that η 1 (t) is time independent, we cannot separate the ETV effects 1 − ǫ from the ARC-520 effects 1 − η 1 (t) . We lump them together, and assume a total drug effect, which ranges between 0.9 < ǫ T < 1 . The other parameter ranges are found as follows. The time-dependent inhibitory effects of treatment on intracellular HBV DNA production, g, was estimated from HBV infected humanized mice treated with NA to range between 0.059 and 0.42 per day. We expand this range by searching over the parameter space 0 < g < 1 . There is a wide range of estimates for the free virus clearance rate in serum: as low as 0.69 per day 25,33,57 ; and as high as 21.7 per day 58 ; we search the entire 0 < c < 100 parameter space. The decay rate of HBsAg is bounded between 0 < d S < 200 per day, containing previous estimates ranging between 0.057 to 0.58 per day 59,60 . In previous modeling work 44,61 HBeAg decay rate d E was set to 0.3 per day. We allow for a larger range 0 < d E < 200 per day, corresponding to half-lives greater than 5 minutes. We assume that the drug absorption rate k eo ranges between 0 < k eo < 1 per www.nature.com/scientificreports/ day. Since ARC-520 was reported to have long lasting effects 56 , we assume a large range for the half-maximal quantity EC i ; between 10 −7 < EC i < 1 mg/kg. These ranges are summarized in Table 1.
Optimization algorithm. We estimate the unknown parameters parm given in Table 1 by minimizing the least squares functional: for each patient. Functional SSQ describes the distance between HBV DNA, HBsAg, and HBeAg titers V data (t i ) , S data (t i ) , E data (t i ) at times t i ( i = {1, . . . , 8} ) and populations V (t i ) , S(t i ) and E(t i ) as given by model (5) at times t i ( i = {1, . . . , 8} ). As described previously (see Eq. (6)), the before treatment titers at t −1 = −8 days are used to determine parameters α , p, p S , p E such that the model's equilibrium matches the titers exactly. Since we assume that the model stays in equilibrium until treatment initiation, we ignore the titers at time t 0 = 0 days. Lastly, it should be noted that we assign the same weight to errors in HBV DNA, HBsAg, and HBeAg. Within the parameter space defined in Table 1, we determine optimal parameter fits for each patient by following four steps (code available upon publication): 1. We create 100 parameter sets using the Latin hypercube samples (LHS) Matlab routine lhsdesign, with random number generator seed two and uniform probability density distribution on each parameter interval.
Since the parameter space spans several orders of magnitude in EC 2 and EC 3 directions, we replace them with EC 2 = 10 EC 2 and EC 3 = 10 EC 3 . Thus, instead of sampling EC 2 and EC 3 in [10 −7 , 1] , we sample EC 2 and Our preliminary work showed that ǫ T ≈ 1 often yields the best results.Therefore, we replace (1 − ǫ T ) = 10 ǫ T and sample ǫ T in the parameter space [−8, −1]. 2. HBV DNA dynamics do not influence HBsAg and HBeAg dynamics. Therefore, we minimize separately over their corresponding parameter sets parm V = {g, c, ǫ T } and We split the LHS into LHS V and LHS S,E containing the respec- Table 1. Variables and parameters in model (5). Parameters indicated by a * are fitted within the given range. www.nature.com/scientificreports/ tive initial parameter guesses and, using Matlab's fmincon routine to minimize SSQ V and SSQ S,E within the parameter space in Table 1, obtain 100 optimal parm V and parm S,E parameter sets. 3. Of the 2 × 100 optimal parameter sets found in part two, we choose the ones yielding minimal SSQ = SSQ V + SSQ S,E , as the overall optimal parameter set for the given patient. 4. To obtain confidence intervals for the optimal parameter estimates p opt for each patient, we employ a bootstrapping technique. We assume that the best fit parameters yield the true dynamics, and that any discrepancy from the data is due to measurement errors. First, we calculate the residuals between simulated functions and measured data at times t i ( i = {1, . . . , 8} ). Next, we create 1000 data sets for the HBV DNA, HBsAg, and HBeAg data at times t −1 , . . . , t 8 , where data at times t −1 and t 0 are as before and data at the remaining times are obtained by adding a randomly drawn residual (with repetition) to the true value at each time, i.e.
where P ∈ {V , S, E} , i = 1, . . . , 8 , and j P,i is drawn at random from {1, . . . , 8} . Lastly, for each data set, we find a new set of optimal parameters by using Matlab's fmincon with initial parameter guess p opt to minimize SSQ V and SSQ SE , as described in (2.). This yields 1000 sets of parameters (one for each data sets), and the confidence intervals on the optimal parameters p opt are obtained as the ranges from the 2.5th percentiles to the 97.5th percentiles of the 1000 parameter values.

Results
Parameter estimates. The best parameter estimates, the respective errors (SSQ) and the the 95% confidence intervals obtained by bootstrapping, are given in Table 2. Numerical solutions for each population versus data are shown in Fig. 1 (see also Figs. 2, 3, and 4 for zoomed in results). Table 3 gives the parameters obtained from equilibrium conditions (6).
Previously reported virus clearance rates range from 0.69 per day 25,33,57 to 21.7 per day 58 . We estimate average virus clearance rates among the five patients c = 3.37 ± 3.38 per day, corresponding to average life-spans of 7.1 h. The fastest free virus clearance rate, c = 9.27 per day (life-span of 2.6 h), occurs in patient 704, who has the lowest pre-treatment virus titer. Assuming 50% of hepatocytes are HBV-infected, we estimate an average intracellular HBV DNA release rate p = 3.21 ± 3.54 per day. Patient 711, who has the highest pre-treatment virus titer, has p = 9.37 per day, 2.9 times higher than the average. Under these estimates, the pre-treatment serum virus production rates, pD 0 , range between 301.5 and 1260 copies/(infected cell×day) for patients 703-710, similar to the 200-1000 copies/(infected cell×day) reported for acute HBV infection 62 . Patient 711, however, has a pre-treatment  (5)  www.nature.com/scientificreports/ serum virus production rate, pD 0 = 4216.5 copies/(infected cell×day), four times larger than in 62 . Intracellular HBV DNA synthesis rates are α = 1755.66 ± 1590.21 copies/(cell× day). As with the serum release rate, patient 711 has 2.6-times higher intracellular HBV DNA synthesis than the average, α = 4526.23 copies/(ml× day). The reported half-life of circulating HBsAg in chronically infected patients is 6.7 days (with a standard deviation of 5.5 days) 59 , which corresponds to HBsAg decay rates 0.057 < d 0,S < 0.58 per day. We estimate average HBsAg decay rates d S = 0.18 ± 0.06 per day, corresponding to HBsAg life-span of 5.6 days for patients 703 and 708-711, and d S = 0.6 per day, corresponding to HBsAg life-span of 1.7 days, for patient 704. The average clearance rates of circulating HBeAg d E = 1.05 ± 0.52 per day, correspond to HBeAg life-spans ranging between   We estimate high efficacy rates, ǫ T > 99.88% , for the combined entecavir and ARC-520 effects in blocking HBV DNA synthesis. The additional time-dependent inhibitory effect on intracellular HBV DNA synthesis is on average g = 0.029 ± 0.018 per day.
The estimated k eo = 0.07 ± 0.021 per day, predicts slow transport of ARC-520 from plasma to liver. The halfmaximal quantities are small, with average log 10 (EC 2 ) = −3.38 ± 0.22 and log 10 (EC 3 ) = −2.98 ± 0.0.33 for the ARC-520 effects on HBsAg and HBeAg, respectively. This implies that the effects of ARC-520 are long-lived, as suggested by Schluep et al. 56 who found that RNA inhibitors persist and induce antiviral effects for longer than the drug's life-span.   www.nature.com/scientificreports/ Pharmacokinetic-pharmacodynamic model dynamics. The predicted HBV DNA populations as given by model (5) for the estimated parameters follow a biphasic decay with short and sharp first phase corresponding to the removal of HBV DNA followed by long and slow second phase decay due to time dependent treatment induced inhibition of intracellular HBV DNA synthesis and infected cell loss. HBsAg and HBeAg decay at steep rates during the first 24.67 ± 10.2 and 7.64 ± 3.95 days, respectively. After reaching minimum values, on average 1.57 ± 0.19 and 1.6 ± 0.33 orders of magnitude smaller than their initial levels, HBsAg and HBeAg rebound (see Figs. 3 and 4). Once the effects of ARC-520 have completely waned, HBsAg and HBeAg decay at rate δ.
In-silico knockout experiments. We are interested in understanding the individual and combined effects of ETV and one-dose of ARC-520 on the dynamics of HBV DNA, HBsAg and HBeAg as given by model (5). We consider the following about the combined ETV and ARC-520 effects on reducing intracellular synthesis, ǫ T : we either attribute it to ETV alone, ǫ T = ǫ ETV T ; or split it between the two effects, ǫ T = ǫ both T . Using the parameters obtained from fitting the combination therapy model (5)    The equations for HBeAg is given by: and for HBeAg is given by: Note that both S(t) and E(t) are independent of ǫ T . HBV DNA follows a biphasic decay with short and sharp first phase corresponding to the removal of free virus followed by a slow second phase decay due to time dependent treatment induced inhibition of intracellular HBV DNA synthesis and removal of infected cells (see Fig. 6, dashed curves). Serum antigen levels remain elevated for all three populations (see Figs. 7 and 8 , dashed curves). When we consider that the treatment that blocks intracellular HBV DNA synthesis, ǫ T , comes from both ETV and ARC-520, we recover the solutions of model (5) for combination therapy given by η 2 = η 3 � = 0 , g = 0 , and ǫ T = ǫ both T � = 0 . Both HBsAg and HBeAg decay at a steep rate during the first 22.7 ± 8.5 and 7.6 ± 4.1 days, respectively. After reaching minimum values, on average 1.5 ± 0.2 and 1.6 ± 0.4 orders of magnitude smaller than their initial levels, HBsAg and HBeAg rebound to their respective ETV monotherapy levels (see Figs. 7 and 8, solid curves). Patient 711 Figure 6. Short-term HBV DNA dynamics under ETV monotherapy (dashed curves), and combined ETV and ARC-520 therapy (solid curves), as given by model (5). Parameters are given in Tables 1 and 2. Additionally, g = 0, ǫ T = ǫ ETV T = 0 and η 2 (t) = η 3 (t) = 0 for ETV monotherapy. Note that both axes are plotted on log scale and that the two graphs overlap. www.nature.com/scientificreports/

Sensitivity of model predictions with respect to changes in the infected cell population's initial condition.
Previous estimates for the percentage of HBV-infected hepatocytes vary between 18 ± 12% in chronic HBsAg carriers 50,51 and 99% in acute infections 26,52 . We have derived our results by assuming that during chronic HBeAg-positive cases half of the liver is infected. Here, we investigate how changes in the size of the initial infected cell population alter our predictions. Analytical investigations show that the dynamics of the viral proteins HBsAg and HBeAg are not influenced by the initial size of the infected cell population, I 0 . After treatment initiation I(t) = I 0 e −δt , and p S = d S S 0 /I 0 and p E = d E E 0 /I 0 (based on the equilibrium assumption (6)). Therefore, the equations for S and E: and are independent of I 0 . Moreover, for p = cV 0 /(ξ D 0 I 0 ) and D 0 = 225/(I 0 /T max ) we find that intracellular HBV DNA D depends on I 0 (see Fig. 9) but HBV DNA in serum does not.

Long-term predictions and the need for uninfected hepatocyte dynamics
We assumed above that infected hepatocytes have a fixed life-span of 250 days. In this section, we are relaxing this assumption and investigate long-term HBV DNA and HBsAg dynamics when increased hepatocyte loss (due to either drug toxicity, or immune-mediated killing) is being considered. When we model it by increasing the infected cell death rate δ in (5) we obtain the following: long-term dynamics of S and E under ETV monotherapy predict that HBsAg decreases below 1 IU/ml 5.32 ± 0.54 months for δ = 7 × 10 −2 per day, 4.21 ± 0.35 years for δ = 7 × 10 −3 per year, and 7.35 ± 0.61 years for δ = 4 × 10 −3 per day, following the initiation of therapy. Since ETV and other nucleoside analogues do not trigger cccDNA removal (and consequently HBsAg and HBeAg removal), the fast loss of HBsAg predicted by model (5) for higher killing rates δ is not realistic.
In this section, we include the dynamics of uninfected and infected cell populations and investigate changes in predictions for increased killing rate δ We incorporate uninfected hepatocytes T which get infected by free virus at rate β , as modeled previously 26,39,63 . Note that we ignore the age of the infection and assume that once a cell becomes infected, it is producing virus (for a PDE model extension in a hepatitis C virus infection, see 64,65 ). Both  www.nature.com/scientificreports/ uninfected and infected hepatocytes proliferate according to a logistic term with maximal growth rate r T and r I and carrying capacity T max . In chronic HBV infections, cccDNA persist under long-term nucleoside analogues treatment 66 . Since the average cccDNA number of untreated HBeAg positive patients is 2.58 copies per infected cell 54 , infected hepatocytes may have two infected off springs. On the other hand, it has been suggested that cccDNA is destabilized by cell division or even lost during mitosis 66 . We account for this by assuming that a fraction of proliferating infected hepatocytes have one infected and one uninfected offspring, and the remaining infected hepatocytes have two infected offsprings. The new model is given by: Liver regenerates rapidly after injury. To account for fast proliferation during chronic disease, we assume that hepatocytes' maximum proliferation rate is r T ≤ 1 per day, and r I = 1 per day, corresponding to doubling time of (up to) 16 h 26,67 . The infectivity rate is at the lower end of previously fitted values 11 , β = 10 −9 IU/(ml× day); we include a death rate for the uninfected hepatocyte population, d T = 4 × 10 −3 per day 68 , identical to that in model (5); and set the fraction of infected hepatocytes that have one uninfected and one infected offspring to = 0.05 . Initial conditions of uninfected and infected hepatocytes are set such that the model is in equilibrium prior to treatment with D 0 = 450 , and V 0 , S 0 , and E 0 as in Table 1. This leads to almost all hepatocytes being infected. Without loss of generality, we investigate the dynamics for patient 703 under combination therapy for a continuum of δ values. Our hypothesis is that NA monotherapy cannot lead to HBsAg loss. In order to obtain infected cell persistence (under NA monotherapy), we need to decrease r T (for a fixed r I = 1 ) as δ increases (a r T − δ threshold required for infected cells persistence is given in Fig. 10). Therefore, HBsAg persistence under increased infected cell killing (as seen in NA treatment) may be explained by high ratio of infected to uninfected cell proliferation. Other events, such as HBV DNA integration, adaptive immune responses, such as cytolytic and non-cytolytic effects, and/or antibody neutralization 11,26 may also explain HBsAg persistence under infected cell (and potentially cccDNA) loss. This is especially true for HBeAg negative patients and NA experienced, HBeAg-positive patients.  Tables 1 and 2

Discussion
Reaching functional cure with current anti-HBV therapies in patients with chronic hepatitis B infection is hindered difficult by the lack of approved direct anti-HBsAg treatment and the presence of large numbers of HBsAg in the blood of infected patients 69,70 . Therapies silencing viral translation through RNA interference technology 17,20,21,71 , inhibiting HBsAg release via nucleic acid polymers [72][73][74] , and inducing neutralization of HBsAg via specific antibodies 75,76 have shown different levels of success 69,70 . Understanding the relative effects in reducing HBV DNA, HBsAg and HBeAg titers of these new approaches alone, and in combination with traditional nucles(t)ide analogues, is particularly important in informing the development of new generation anti-HBsAg therapies.
To help in this endeavor, we developed mathematical models describing the HBV DNA, HBsAg and HBeAg in the presence of a silencing RNAi drug called ARC-520. We used the models and clinical trial data from treatment naive, HBeAg-positive patients that receive a one time ARC-520 injection and daily nucleoside analogue treatment with entecavir 20 , to determine the efficacy of ARC-520 and nucleoside therapies on the short and long-term dynamics of HBV DNA, HBsAg, and HBeAg. To the best of our knowledge, we report for the first time that the time-dependent ARC-520 effects on HBsAg and HBeAg are more than 96% effective around day 1, and slowly wane to 50% in 1.8-3.4 months and 1.5-3.5 months, respectively. The combined ARC-520 and entecavir effect on HBV DNA is constant over time, with efficacy of more than 99.8% , which is similar to other nucleoside analogues trials.
A simplified version of the model, which ignored the dynamics of hepatocyte proliferation and infection, was sufficient to explain the short-term (about 100 days) dynamics observed in five patients in the current study. In the long-term, however, infected cells may die at faster rates, due to either drug toxic effects or increased immune killing. Lowering infected hepatocyte's life-span to 100 (10) days, however, resulted in fast HBsAg removal, with decay below 1 IU/ml in 4.2 years ( 5.3 months). This loss, however, was in contradiction with clinical reports of low percentages of patients clearing HBsAg during long-term nucleoside analogues treatment 6 , suggesting that more complex models are needed for long-term (several years) predictions.
To determine under what conditions increased infected cells death does not spill over into unrealistic HBsAg and HBeAg loss under long-term nucleoside analogue therapy, we extended model (5) to include infected and uninfected cell dynamics. We assumed lower infected cells life-span (100 and 10 days), included division of both infected and uninfected populations, and determined that long-term HBsAg and HBeAg persistence under long-term HBV DNA clearance can be explained by high ratios of infected to uninfected division rates. Therefore, high ratio of infected to uninfected division rates, which correspond to the infection of the entire liver and may be indicative of scenarios where HBsAg seroclearance will not happen. Interestingly, we and others have associated high ratios of infected to uninfected division rates to triphasic HBV DNA decay under treatments with nucleoside analogues, a sign of suboptimal drug response 33,35 . Whether infected hepatocytes indeed proliferate faster than uninfected hepatocytes remains under investigation.
While modeling results suggest that one-dose of ARC-520, in combination of daily entecavir, has limited long-term effects, we did not consider whether a transient reduction of HBsAg and HBeAg leads to the appearance of anti-HBs or anti-HBe antibodies, removal of immune-exhaustion, and eventual functional cure. Recent studies found that large levels of HBsAg might cause dysfunctional programming of HBsAg-specific B cells through persistent stimulation 77 . It has been suggested that therapeutic vaccines containing one (PreS2) or two (PreS1 or PreS2) envelope proteins together with serum HBsAg reducing drug therapies are needed in order to induce high levels of anti-HB antibodies, which may correlate with functional cure [78][79][80] . We ignored the level of immune modulation following RNAi based therapy, such as cytolytic and non-cytolytic T cell functions and antibody responses, which is a model limitation, and therefore, we cannot say whether such effects were induced at higher rates during the transient HBsAg loss.
Our study has limitations. We only used the data on HBeAg-positive patients (cohort 7 in 20 ) since they best responded to ARC-520 therapy. Moreover, we did not model HBV DNA integration, which has been reported as a source of HBsAg production, especially in HBeAg-negative and NA-experienced HBeAg-positive patients with low cccDNA 20 . As kinetic HBV data from next generation RNAi therapy capable of inducing stronger HBsAg reduction in both HBeAg-negative and HBeAg-positive patients becomes available 21,81,82 , we aim to adapt our modeling framework to include HBV DNA integration.
In conclusion, we developed a mathematical model and used it together with patient data, to estimate the time-dependent ARC-520 efficacies in blocking HBsAg and HBeAg productions. Additional data and theoretical efforts are needed to determine whether RNAi therapies have a feedback effect on the reversal of immune exhaustion, immunomodulatory immune responses, and potential functional cure. www.nature.com/scientificreports/