A Mathematical Model for early HBV and -HDV Kinetics during Anti-HDV Treatment

Hepatitis delta virus (HDV) is an infectious subviral agent that can only propagate in people infected with hepatitis B virus (HBV). HDV/HBV infection is considered to be the most severe form of chronic viral hepatitis. In this contribution, a mathematical model for the interplay between HDV and HBV under anti-HDV treatment is presented. Previous models were not designed to account for the observation that HBV rises when HDV declines with HDV-specific therapy. In the simple model presented here, HDV and HBV kinetics are coupled, giving rise to an improved viral kinetic model that simulates the early interplay of HDV and HBV during anti-HDV therapy.


Introduction
Globally, ~15-72 million people are currently infected with the hepatitis delta virus (HDV) [1,2]. HDV is an infectious subviral agent that can only propagate in individuals infected with hepatitis B virus (HBV), which supplies the necessary envelope proteins needed to assemble infectious HDV progeny virions. HDV was recognized as a distinct agent in 1977 [3,4], and is the smallest known viral genome that infects humans. HDV causes an accelerated course of liver disease compared to HBV alone [1] and is the most severe form of chronic viral hepatitis in humans [5][6][7].
There is currently no FDA-approved therapy for HDV. Pegylated interferon-α (IFN), which is the recommended treatment by the American Association for the Study of Liver Diseases (AASLD) is associated with very low response rates and high relapse rates, even when administered for 5 years [8,9]. There is little information on the kinetics of HDV infection or the critical interplay between HDV, HBV, and the host, which are needed to understand the mechanism of action and optimize emerging antiviral therapies for HDV. Importantly, the recent development of new HDV experimental systems and anti-HDV drugs that are now in clinical trials provide a unique opportunity to study, for the first time, the dynamic interaction between HDV and HBV. Herein, we describe a simple model for the HDV and HBV interplay under anti-HDV treatment.

Background
A mathematical model for HDV kinetics during antiviral treatment was described in [10,11] without considering HBV dynamics. Subsequently, mathematical models were developed that account for both HDV and HBV [12,13]. These mathematical models for HBV-HDV dynamics aimed to theoretically simulate acute infection or antiviral treatment with nucleos(t)ides analogs or IFN after chronic infection established in silico [12,13]. The existing models were not designed to explain anti-HDV drugs that solely affect HDV.
The model of de Sousa and Cunha [12] was the first HDV-HBV model, which included viral load dynamics in response to therapy by lamivudine and/or IFN, where both HBV and HDV viral loads decrease post-therapy as depicted in . The model contains six equations: Zakh et al. Page 2 dv 1 dt = (1 − ε)[k 1 y 1 (t) + k 3 y 3 (t)] − u 1 v 1 (t) (5) dv 2 dt = k 2 y 3 (t) − u 2 v 2 (t) (6) where x(t) is the number of uninfected cells at time t, y 1 (t) is the number of HBV-infected cells at time t, y 2 (t) is the number of HDV-infected cells at time t, y 3 (t) is the number of infected cells with both HBV and HDV at time t, v 1 (t) is the HBV viral load at time t, and v 2 (t) is the HDV viral load at time t. Both the variables and all of the parameters of the model are shown in Table 1.

x(t) is the number of uninfected cells at time t, y(t) is the number of HBV-onlyinfected cells at time t, z(t) is the number of HDV-only-infected cells at time t, w(t) is the number of HBV-HDV-coinfected cells at time t, and T(t) is the sum of all these cells at time t.
Both the variables and all of the parameters of the model are shown in Table 2.
At the time both the above models were published, therapy to suppress HDV alone was not available. Bulevirtide (also known as Myrcludex B or Hepcludex), an HDV entry inhibitor, can profoundly suppress HDV RNA and was recently approved by the European regulatory bodies for the treatment of chronic hepatitis delta [14]. In addition, the prenylation inhibitor Lonafarnib (LNF), an investigational oral HDV drug in clinical trials, was shown to suppress HDV in a proof-of-concept phase 2A study. Our new mathematical model, presented here, incorporates the observation that HBV viral loads can increase under LNF monotherapy [15].

Materials and Methods
Based on newly published empirical results [15], we present a new mathematical model for the interaction between HBV and HDV that can explain the coinfection pattern of HBV increasing as HDV decreases. A schematic diagram for the model is depicted in Figure 1.
The model comprises only two equations: dB dt = p 2 I 0 n −4 D 0 D(t) n − cB(t) (13) where D is the HDV viral load and B is the HBV viral load. Both the variables and all of the parameters of the model are shown in Table 3. We keep the number of infected cells, I 0 , constant, i.e., ignoring the dynamics of infected cells, such as infection of susceptible cells and/or infected cell proliferation and/or death (assuming a long-infected cell life-span as previously performed [16,17]). Treatment is assumed to solely block HDV production rate p 1 by (1 -ε), where 0 ≤ ε ≤ 1 represents its efficacy. To account for the biphasic HDV RNA decay in the absence of infected cell dynamics, we assume that treatment has additional time-dependent inhibitory effects on HDV production and we model it by decreasing p 1 further by e −gt where g ≥ 0 is a constant and t is the time in days post-treatment initiation, as previously performed [16,18] To accommodate the steady-state solution, the following two equations hold: p 2 = cB 0 I 0 (15) where it should also be noted that, prior to treatment onset, ε = 0, g = 0, n = 1.
It is possible to obtain an analytical solution for D(t). The exact solution can be derived as follows: dD dt = (1 − ε)p 1 I 0 e −gt − cD(t) (16) Multiplying both sides of the equation by e ct , one obtains: dD dt e ct + cD(t)e ct = (1 − ε)p 1 I 0 e −gt e ct (17) Therefore, the following equation is obtained: Integrating both sides yields: De ct = (1 − ε)p 1 I 0∫ e −gt e ct dt (19) The solution for the integral is: It then follows that the general solution for D is: In order to find C 1 , at time t = 0, D = D 0 , one plugs in the initial condition and C 1 is obtained: Therefore, the exact solution for D can now be written in full by the following expression: It could be possible to extract the exact solution for B(t) as well, by a formidable derivation using hypergeometric functions, as was preliminary checked with Mathematica by Wolfram Research, but the result is likely less intuitive, and this is left for future work. In practice, the numerical solution can be used for simulations with the model as performed next.

Results
Given the data shown in Figure 5 by Yurdaydin et al. [15], we digitized the data points with [19] and performed simulations with our new model. Unlike applications such as in [20,21] where more sophisticated numerical methods are needed to solve the differential equations, here, the model is simple enough to utilize a standard Runge-Kutta scheme of the fourth order where numerical simulations are performed using Berkeley Madonna version 8.3.18. More specifically, Patient 1 ( Figure 2) was treated with LNF 200 mg twice a day (BID) and Patient 2 ( Figure 3) was treated with LNF 300 mg BID ( Figures 5A,B in [15], see digitized data in the supplementary material Tables S1 and S2). Model parameter values for Patient 1 and Patient 2 show in Table 4. Figures 2 and 3 depict the HBV/HDV levels after LNF treatment initiation and are in agreement with [15]. While HDV decreases, HBV increases, which is not captured in previous models [12,13].

Discussion and Conclusions
A new model is presented for HBV-HDV interaction, which captures the kinetics seen during anti-HDV treatment with LNF. These observations are in agreement with data showing that the HBV viral load is suppressed by the presence of HDV [22]. The model is a simplified one and contains only two equations for HDV and HBV dynamics. The HDV equation contains an exponent not included in the HBV equation (i.e., e −gt ), which represents the time-dependent blocking of viral production that explains the second-phase decline observed with LNF. The HBV equation is inversely dependent on HDV as per the empirical results, and includes patient-specific variables. The clearance rate of HDV and HBV from circulation is considered similar to that in [11]. We assume that the number of infected cells is constant, a feature that we intend to relax in future models.
Most studies in chronically infected patients show the levels of both viruses to be in a steady state. Typically, three HDV/HBV steady-state profiles are observed, including predominant HDV replication, similar rates of replication for both viruses, and less commonly, predominant HBV replication. Patients exhibiting predominant HDV replication before treatment initiation (as the two patients shown in Figures 2 and 3) comprise the majority of cases [23,24].
The model presented here was made possible by the development of new antiviral agents for HDV. The model provides an important first step toward understanding the dynamic interaction between HDV and HBV, which is needed to achieve a detailed understanding of viral kinetics and treatment response in patients with HDV/HBV infection. In future work, our simplified model can be augmented to address additional features such as immune response and uninfected and infected cell dynamics. Moreover, the previous models [12,13] can be further developed to account for the interplay between HDV and HBV incorporated in the new model presented in the current study.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Funding:
This research was funded by the NIH grants R01AI144112 and R01AI146917.

Data Availability Statement:
The data presented in this study are available in the supplementary material (digitized data is provided in Tables S1 and S2).  The variables and the parameters in the model of de Sousa and Cunha [12].