A Coupled Mathematical Model of the Intracellular Replication of Dengue Virus and the Host Cell Immune Response to Infection

Dengue virus (DV) is a positive-strand RNA virus of the Flavivirus genus. It is one of the most prevalent mosquito-borne viruses, infecting globally 390 million individuals per year. The clinical spectrum of DV infection ranges from an asymptomatic course to severe complications such as dengue hemorrhagic fever (DHF) and dengue shock syndrome (DSS), the latter because of severe plasma leakage. Given that the outcome of infection is likely determined by the kinetics of viral replication and the antiviral host cell immune response (HIR) it is of importance to understand the interaction between these two parameters. In this study, we use mathematical modeling to characterize and understand the complex interplay between intracellular DV replication and the host cells' defense mechanisms. We first measured viral RNA, viral protein, and virus particle production in Huh7 cells, which exhibit a notoriously weak intrinsic antiviral response. Based on these measurements, we developed a detailed intracellular DV replication model. We then measured replication in IFN competent A549 cells and used this data to couple the replication model with a model describing IFN activation and production of IFN stimulated genes (ISGs), as well as their interplay with DV replication. By comparing the cell line specific DV replication, we found that host factors involved in replication complex formation and virus particle production are crucial for replication efficiency. Regarding possible modes of action of the HIR, our model fits suggest that the HIR mainly affects DV RNA translation initiation, cytosolic DV RNA degradation, and naïve cell infection. We further analyzed the potential of direct acting antiviral drugs targeting different processes of the DV lifecycle in silico and found that targeting RNA synthesis and virus assembly and release are the most promising anti-DV drug targets.

Introduction

Dengue virus (DV) is the most prevalent vector-borne virus in the world, with an estimated global number of 390 million new infections annually, thereof 90 million with clinical manifestations, including severe dengue disease (Bhatt et al., 2013). DV poses a huge burden on human populations and health systems in affected countries, and significantly impacts the economy in tropical countries, including the southern United States (WHO, 2012). Fueled by climate change and globalization and accelerated further by viral evolution, the expansion of DV is expected to increase further (Murray et al., 2013). DV is transmitted mainly by female Aedes mosquitos, and with the spread of its vector, DV is spreading as well (Campbell et al., 2015). In consequence, the global incidence of DV infection has already risen 30-fold during the past 50 years. Infection with DV causes flu-like symptoms but is occasionally associated with severe complications. The fatality rate of dengue infection is between 1 and 5%, and below 1% with proper symptomatic treatment (Ranjit and Kissoon, 2011). There is no antiviral therapy available against DV, and the recently approved vaccine has limited efficacy and depends on baseline serostatus of the vaccine recipient (World Health Organization, 2016).

DV infects dendritic cells (DC), B cells, T cells, monocytes, macrophages, but also the liver. DV is an enveloped, positive-sense (+)RNA virus of the Flaviviridae family within the genus Flavivirus, consisting of four distinct serotypes (DV1/2/3/4). Its approximately 10.7 kb genome encodes the three structural proteins capsid (C), precursor membrane (prM), and envelope (E) protein and seven non-structural proteins (NS1, NS2A, NS2B, NS3, NS4A, NS4B, and NS5). Upon entry into the host cell, the viral RNA genome is translated at the rough endoplasmic reticulum (rER) giving rise to a polyprotein, ~3,400 amino acids in length, which is co- and post-translationally cleaved by viral and host proteases into the structural and non-structural proteins (Neufeldt et al., 2018). DV induces membrane alterations at the rER, forming membrane invaginations. The viral RNA genome is amplified in these replication compartments (RC), starting with minus-strand synthesis to obtain a double-stranded RNA (dsRNA) intermediate which then serves as template for further plus-strand synthesis. The newly synthesized viral (+)RNA genomes leave the RC and are then either packaged into virions, which after maturation are released from the infected cell, or are again used for the next round of viral RNA translation (Bartenschlager and Miller, 2008; Rodenhuis-Zybert et al., 2010; Tuiskunen Bäck and Lundkvist, 2013; Screaton et al., 2015).

The host cell's defense against DV is mediated via pattern recognition receptors (PRRs), in case of DV mainly via the endosomal Toll-like receptors (TLR3/TLR7/TLR8) and the cytosolic RNA helicases retinoic acid inducible gene I (RIG-I) and melanoma differentiation associated gene 5 (MDA-5) (Muñoz-Jordán and Fredericksen, 2010; Morrison et al., 2012). TLR3 recognizes dsRNA, while TLR7 and TLR8 recognize viral single-stranded RNA (Xagorari and Chlichlia, 2008). All three TLRs activate signaling cascades that lead to the production of interferon α/β (IFN α/β) and inflammatory cytokines. RIG-I/MDA-5 signals via mitochondrial antiviral signaling protein (MAVS) and tumor necrosis factor receptor-associated factor 3 (TRAF3), activating tank-binding kinase 1 (TBK1) and ultimately phosphorylating interferon regulatory factor 3 (IRF3) and activating nuclear factor kappa B (NF-κB). The subsequent type I (α/β) and type III (λ) IFN production induces the activation of hundreds of IFN stimulated genes (ISGs), bringing the cells into an antiviral state and resulting in an inhibition of DV (Nasirudeen et al., 2011; Tuiskunen Bäck and Lundkvist, 2013; Dalrymple et al., 2015).

DV, however, is not defenseless, and has evolved a number of mechanisms antagonizing the antiviral response of the cell both at the level of activation of the host cell immune response (HIR) and the induced effector phase. For instance, 2'-O-methylation of the DV RNA genome, mediated by NS5, was shown to slow down the activation kinetics of the IFN response (Schmid et al., 2015). In addition, the DV NS2B-3 protease cleaves the stimulator of interferon genes (STING), thus reducing type I IFN production (Diamond and Pierson, 2015). In fact, several groups have shown that the suppression of the early IFN induction by DV is critical for successful virus infection and replication (Shresta et al., 2004; Perry et al., 2009). Moreover, Schmid et al. (2015) have shown that the ability of IFN to control DV spread might be stochastic and “leaky.” While secreted IFN protects surrounding naïve cells from infection, this protection is incomplete with cells infected with DV prior to activation of the IFN response (Schmid et al., 2015). DV replication occurs inside membrane vesicles corresponding to invaginations into the rER lumen, likely shielding viral dsRNA intermediates from recognition by the HIR (Welsch et al., 2009). At the level of the effector phase, DV NS5, which contains the enzymatic activity for capping and amplification of the viral RNA genome, was shown to bind to and induce the degradation of the signal transducer and activator of transcription factor (STAT) 2 via a proteasome-dependent mechanism (reviewed in Neufeldt et al., 2018), thus blocking ISG induction downstream of the IFN receptor. Therefore, the interplay between DV and the innate immune response (IIR) is complex, and its exact magnitude and dynamics likely impact and possibly determine clinical outcome of the infection.

Mathematical modeling is a valuable tool to study complex dynamical systems and has been employed to analyze infection dynamics for a number of different viruses (Zitzmann and Kaderali, 2018). Most previous work on modeling viral infection has built on the basic model introduced by Nowak and Bangham (1996) and Nowak et al. (1996), focusing on the dynamics of susceptible cells, infected cells, and virus at the cell population level. Especially, the within-host dynamics of human immunodeficiency virus (HIV) and hepatitis C virus (HCV) have been studied in detail with simple models on the cell population scale with regard to the antiviral immune response and treatment opportunities (Ho et al., 1995; Wei et al., 1995; Perelson et al., 1996, 1997; Bonhoeffer et al., 1997; Neumann, 1998; Stafford et al., 2000; Perelson, 2002; Rong and Perelson, 2009; Perelson and Ribeiro, 2013; Perelson and Guedj, 2015). In the case of DV modeling, these so-called target cell-limited models have been linked to the adaptive immune response via modeling of antibodies and T cell responses (Clapham et al., 2014; Ben-Shachar et al., 2016), and the IIR by IFN (Ben-Shachar and Koelle, 2014; Schmid et al., 2015). Several authors have furthermore developed intracellular replication models for related viruses (Dahari et al., 2007; Heldt et al., 2012, 2013; Binder et al., 2013; Guedj et al., 2013; Clausznitzer et al., 2015; Laske et al., 2016; Benzine et al., 2017; Aunins et al., 2018), but to our knowledge there is no mathematical model describing the intracellular steps of DV replication to date. In this manuscript, we focus on the highly dynamic initial phase post cell infection and developed a detailed differential equations model capable of quantitatively describing the intracellular infection dynamics. We measured viral replication in two different cell lines, Huh7 cells (with very little HIR) and A549 cells (with high HIR competence). By integrating the main steps of the HIR into the model, we were able to describe the infection kinetics in both cell types. Our investigation focuses on the cell line-specific impact of host factors, which determine RNA synthesis efficiency, and the involvement of the HIR. Using our mathematical model, we identified possible antiviral modes of action of the HIR on the DV lifecycle and the viral countermeasures suppressing DV recognition and activation of the HIR. We further identified the most sensitive processes in the DV lifecycle, which might constitute promising antiviral drug targets, and we evaluated possible antiviral intervention strategies in silico.

Materials and Methods

Measuring DV Replication in Huh7 and A549 Cells

Cell Lines

A549 and Huh7 cells were cultivated at 37°C, 5% CO² in of DMEMcplt (2 mmol/L L-glutamine, non-essential amino acids, 100 U/ml penicillin, 100 μg/ml streptomycin and 10% fetal calf serum).

Kinetics Experiments (120830 and 120921)

2 × 10⁵ A549 and Huh7 cells were seeded 1 day prior infection. Cells were infected with DV reporter virus expressing Renilla Luciferase (Schmid et al., 2015) at a MOI of 10. After 1 h, the inoculum was removed and cells washed thrice with sterile PBS prior addition of DMEMcplt. Cells were incubated at 37°C for the indicated time points.

Infectivity Titers

Supernatants were harvested and filtered through a 0.45 μm-pore size membrane. Supernatants were supplemented with 15 mM HEPES and stored at −80°C. Infectivity titers of virus supernatants were determined by limiting dilution assay using Huh7.5 cells as described elsewhere (Lindenbach et al., 2005).

Interferon Lambda ELISA

Supernatants of infected cells were supplemented with 1% (v/v) Triton X-100 to inactivate DENV infectious particles and subsequently stored at −80°C until further use. Interferon lambda protein release was determined by VeriKine-DIYTM Human Interferon Lambda/IL-28B/29/28A ELISA (PBL Interferon Source, USA) with an assay range of 62.5 to 4,000 pg/ml. ELISA procedure was conducted according to the manufacturer's protocol using high binding 96 well ELISA microplates (Greiner Bio-One, Frickenhausen, Germany). The optical density of each well was determined immediately using a microplate reader set to 450 nm.

RT-qPCR

Cells were lysed for RNA extraction and subsequent qRT-PCR analysis by adding 350 μl RA1 lysis buffer (Macherey-Nagel, Düren, Germany) supplemented with 1% β-mercapto-ethanol and stored at −80°C. Total RNA was extracted using the NucleoSpin RNA II Kit (Macherey-Nagel) as recommended by the manufacturer. RT-qPCR was described elsewhere (Schmid et al., 2015). The following primers were used: DENV2 (forward 5′-GCC CTT CTG TTC ACA CCA TT-3′; reverse 5′-CCA CAT TTG GGC GTA AGA CT-3′); IFIT1 (forward 5′- GAA GCA GGC AAT CAC AGA AA-3′; reverse 5′-TGA AAC CGA CCA TAG TGG AA-3′). GAPDH mRNA (primer forward 5′ - GAA GGT GAA GGT CGG AGT C−3′; reverse 5′ - GAA GAT GGT GAT GGG ATT TC – 3′) was used for normalization of input RNA.

Luciferase Assay

Cells were lysed in 100 μl Luciferase lysis buffer (1% Triton X-100, 25 mM glycyl-glycine, 15 mM MgSO₄, 4 mM EGTA, and 1 mM DTT, pH 7.8) and stored at −80°C. For detection of Renilla luciferase activity, 20 μl lysate was mixed with 100 μl assay buffer (25 mM glycyl-glycine, 15 mM MgSO₄, 4 mM EGTA, and 15 mM potassium phosphate pH 7.8) containing 1.4 μM coelenterazine (P.J.K). All measurements were done in duplicate by using a tube luminometer (Berthold, Pforzheim, Germany). Replication efficiency was determined by normalization to the 2 h values reflecting infection efficiency.

Mathematical Model

We developed a mechanistic mathematical model using ordinary differential equations (ODEs) and mass action kinetics to analyze the interplay between DV and the HIR. We have previously published a detailed intracellular replication model for HCV (Binder et al., 2013). Both DV and HCV are closely related (+)RNA viruses of the same family, and their intracellular replication proceeds in similar steps. We therefore used the mathematical model of HCV replication as the basis for the model structure of the current DV model. This model was then extended to the full virus lifecycle by adding infection and assembly of new virus particles to the model. We then adapted the model to DV where necessary (see details below), refitted model parameters on the Huh7/DV infection data and complemented the replication model by an HIR sub-model that comprises the key components of the IIR. Figure 1 gives an overview of the DV replication sub-model and Figure 2 shows the HIR sub-model; we describe each of the main model components in the following.

FIGURE 1

Figure 1. Graphical illustration of the DV replication model. (1) The Virus (V) attaches to a permissive cell (k_a) and (2) enters via receptor-mediated endocytosis (k_e). (3) The virus in the endosome (V_E) is degraded with rate (μ_VE). The viral and endosomal membrane fuse (k_f) and release the viral RNA genome (R_V), which is degraded with rate μ_RV. (4) Ribosomes (Ribo_DV) bind (k₁) at the viral RNA genome, forming a translation complex (TC), which in turn is degraded with rate μ_TC. (5) The viral genome is translated (k₂) into a long polyprotein (P_P) and (6) subsequently cleaved (k_c) into structural (P_S) and non-structural proteins (P_N), which degrade with rate μ_PS and μ_PN, respectively. During the translation, luciferase (L) is produced as a marker for translation activity and is degraded with rate μ_L. (7) The TC together with non-structural proteins and host factors (HF) initiating the formation (k_Pin) of a replication complex (R_RC). (8) The antisense synthesis (k_4m) leads to production of double-stranded RNA (dsRNA) that (9) binds to non-structural proteins in the Replication Compartment RC (P_RC) forming (k₅) a minus-strand RNA intermediate complex (R_IDS). R_IDS in turn (10) initiates plus-strand RNA synthesis (k_4p). (11) The newly synthesized plus-strand RNA (R_P) can be exported out of the RC into the cellular cytoplasm (k_Pout), (12) undergoes another round of replication (k₃) or (13) is packed into virus particles and released (k_r) from the host cell (14), where it can infect naïve cells (k_re) for a further round of replication. Species in the RC degrade with constant rate μ_RC, (15) except the dsRNA species which get transported out the RC with rate constant μ_RC resulting in an accumulation of dsRNA in the cytoplasm (R_R). Extracellular IFN and ISG proteins have an antiviral effect (AE⊣ see Figure 2) on the DV lifecycle via different mechanisms (k₁, k_re, and μ_RV) while the P_N dependent countermeasures (VC⊣) targeting HIR pathways (k_rig and k_jak, see Figure 2).

FIGURE 2

Figure 2. Graphical illustration of the HIR sub-model. (16) Upon recognition of cytoplasmic dsRNA (R_R) via the RIG-I pathway (k_rig), IFN (F_IN) is produced that (17) is released from the cell (k_s). (18) Extracellular IFN (F_EX) activates the JAK/STAT pathway (k_jak) that leads to the production of ISG mRNA (I_R). (19) Ribosomes (Ribo_HIR) bind to ISG mRNA in order to form (k_IC) a translation initiation complex (IC) that (20) translates (k_t) the ISG mRNA into ISG proteins (I_P). Extracellular IFN and ISG proteins have an antiviral effect (AE⊣) on the DV lifecycle via different mechanisms (k₁, k_re, and μ_RV, see Figure 1) while the P_N dependent countermeasures (VC⊣, see Figure 1) targeting HIR pathways (k_rig and k_jak).

DV Replication

Our DV replication sub-model (Figure 1) is composed of four main processes in the DV lifecycle: (1) Binding of DV particles to the cell surface and viral entry via endocytosis. (2) Uncoating and release of the viral RNA genome into the cytoplasm, followed by translation into a polyprotein that is subsequently cleaved into the structural and non-structural viral proteins. (3) The non-structural viral proteins initiate the formation of a RC, in which the viral RNA replication takes place via a dsRNA intermediate. (4) The newly synthesized (+)RNA can then be used as a template for further RNA replication, for protein production at the ribosomes, or it is used together with the structural proteins to assemble new virus particles, which are then released from the cell. The cycle thereafter starts over again with further infection of naïve cells.

We model the infection process by the following three ODEs, where V describes extracellular virus, V_A is virus that has attached to the host cell but is still at the cell surface, and V_E is virus that has been endocytosed:

$\begin{array}{ll}\frac{dV}{dt}=k_{re}{V}_{INF}-k_{a}V-{\mu }_{V}V-wV& (1)\end{array}$

$\begin{array}{ll}\frac{d{V}_A}{dt}=k_{a}V-k_e{V}_A& (2)\end{array}$

$\begin{array}{ll}\frac{d{V}_E}{dt}=k_e{V}_A-k_f{V}_E-{\mu }_{V_E}{V}_E& (3)\end{array}$

Equation (1) and (2) describe how extracellular virus (V) binds at the cell surface of a permissive host cell with rate constant k_a and degrades with rate constant μ_V (Equation 1). To keep the model simple, we assume that attachment is non-reversible. According to the experimental set-up, the cells were washed to remove unbound virus from the initial infection. This is considered in the model through the term wV, where w is modeled as

$\begin{array}{ll}w={\omega }_s\frac{1}{\sqrt{2\pi {\omega }_d^2}}{e}^{-\frac{{(t-{\omega }_t)}^2}{2{\omega }_d^2}},& (4)\end{array}$

with washing time point ω_t, washing duration ω_d, and washing strength ω_s (for more details, see Supplementary Material). The term k_reV_INF in Equation (1) describes newly released infectious virus particles from previous rounds of virus replication, which can infect naïve cells. Equations (2) and (3) describe cell entry: Attached virus (V_A) enters the cell with rate constant k_e via receptor-mediated endocytosis. Subsequently, virus in the endosome (V_E) undergoes conformational changes of the nucleocapsid, leading to fusion of the viral and endosomal membranes with rate constant k_f (Equation 3). Endocytosed virus V_E decays with rate constant μ_VE.

The steps associated with RNA translation are described by the following ODEs, which are based on our previously published HCV replication model (Binder et al., 2013). Here, R_V describes viral (+)RNA in the cytoplasm, TC corresponds to active translation complexes, P_P describes viral polyprotein, and P_S and P_N are representatives for the structural and non-structural viral proteins, respectively. In our experimental data, P_P is measured via a bicystronic luciferase reporter system, the species L is therefore required for model fitting purposes and gives luciferase protein concentration, which degrades with rate constant μ_L. The following ODEs describe the temporal evolution of these species:

$\begin{array}{l}\frac{d{R}_V}{dt}=k_f{V}_E-k_1{R}_V({Rib{o}_{DV}}_{tot}-TC)+k_{2}TC\\ +k_{Pout}{R}_P-{\mu }_{R_V}{R}_V& (5)\end{array}$

$\begin{array}{l}\frac{dTC}{dt}=k_1{R}_V({Rib{o}_{DV}}_{tot}-TC)-k_{2}TC\\ -k_{Pin}{P}_{N}TC(H{F}_{R{C}_0}-R_{RC})-{\mu }_{TC}TC& (6)\end{array}$

$\begin{array}{ll}\frac{d{P}_P}{dt}=k_{2}TC-k_c{P}_P& (7)\end{array}$

$\begin{array}{ll}\frac{dL}{dt}=k_{2}TC-{\mu }_{L}L& (8)\end{array}$

$\begin{array}{ll}\frac{d{P}_S}{dt}=k_c{P}_P-{\mu }_{P_S}{P}_s-N_{P_S}{v}_p& (9)\end{array}$

$\begin{array}{ll}\frac{d{P}_N}{dt}=k_c{P}_P-k_{Pin}{P}_{N}TC(H{F}_{R{C}_0}-R_{RC})-{\mu }_{P_N}{P}_N& (10)\end{array}$

After cell entry, the viral RNA genome R_V is released into the cytoplasm and is subsequently translated giving rise to viral protein or is degraded with rate constant μ_RV (Equation 5). Free ribosomes (Ribo_DV) reversibly bind to the viral RNA genome (R_V) to form a translation initiation complex (TC, Equation 6) with rate constant k₁. We assume that the total number of ribosomes, Ribo_DVtot, is constant, hence free ribosomes available for translation are given by Ribo_DV(t) = Rib_oDV0 − TC and it is not necessary to introduce a separate equation for the ribosomes. HF_RC0 represents one or more unspecified host cell factor(s) that are required for formation of the RC and is assumed constant (see Table S1). Upon translation of the viral genome into a polyprotein (P_P, Equation 7) with rate constant k₂, the translation initiation complex (TC) dissociates into free viral RNA (R_V) and ribosomes (Ribo_DV). Furthermore, during TC degradation with rate constant μ_TC, ribosomes dissociate from the complex: Ribo_DV + R_V → TC → Ribo_DV + R_V + P_P + L. We measure polyprotein production using a luciferase reporter system, and hence include luciferase (L, Equation 8) into the mathematical model. L is produced with rate constant k₂ and is degraded with rate constant μ_L. The polyprotein (P_P) is cleaved into structural (P_S) and non-structural proteins (P_N) with rate constant k_c, which degrade with rate constants μ_PS and μ_PN, respectively. We note here further that later in the virus lifecycle, the structural proteins (P_S, Equation 9) are packed together with newly synthesized (+)RNA (R_P) into virions, thus the corresponding term involving v_p in equation (9), which is detailed in equation (17) below. Furthermore, the non-structural proteins (P_N), e.g., the RNA-dependent-RNA polymerase, are required for viral RNA synthesis (compare Equation 11).

Equations (11) to (16) describe the viral RNA replication inside the replication vesicles. The species R_RC describes a replication intermediate complex for minus strand RNA synthesis inside the RC, whereas R_IDS is the corresponding intermediate complex for (+)RNA production. R_DS and R_P describe dsRNA and single (+)RNA in the RC. P_RC corresponds to the viral RNA polymerase in the RC that is required for RNA replication. Replication is thus modeled by the following equations:

$\begin{array}{l}\frac{d{R}_{RC}}{dt}=k_{Pin}{P}_{N}TC(H{F}_{R{C}_0}-R_{RC})-k_{4m}{R}_{RC}\\ +k_3{R}_P{P}_{RC}(H{F}_{R{C}_0}-R_{RC})-{\mu }_{RC}{R}_{RC}& (11)\end{array}$

$\begin{array}{ll}\frac{d{R}_{DS}}{dt}=k_{4m}{R}_{RC}-k_5{R}_{DS}{P}_{RC}+k_{4p}{R}_{IDS}-{\mu }_{RC}{R}_{DS}& (12)\end{array}$

$\begin{array}{ll}\frac{d{R}_{IDS}}{dt}=k_5{R}_{DS}{P}_{RC}-k_{4p}{R}_{IDS}-{\mu }_{RC}{R}_{IDS}& (13)\end{array}$

$\begin{array}{ll}\frac{d{R}_R}{dt}={\mu }_{RC}(R_{DS}+R_{IDS})-{\mu }_{R_V}{R}_R& (\mathrm{14a})\end{array}$

$\begin{array}{l}\frac{d{P}_{RC}}{dt}=k_{4m}{R}_{RC}-k_3{R}_P{P}_{RC}(H{F}_{R{C}_0}-R_{RC})\\ -k_5{R}_{DS}{P}_{RC}+k_{4p}{R}_{IDS}-{\mu }_{RC}{P}_{RC}& (15)\end{array}$

$\begin{array}{l}\frac{d{R}_P}{dt}=k_{4p}{R}_{IDS}-k_3{R}_P{P}_{RC}(H{F}_{R{C}_0}-R_{RC})\\ -k_{Pout}{R}_P-v_p-{\mu }_{RC}{R}_P& (16)\end{array}$

By analogy to HCV, we assume here that the initiation of RNA replication occurs from freshly translated viral RNA, hence Equations (6) and (11) model RNA import into the RC from TC (instead of R_V) with rate constant k_Pin. At the same time, non-structural proteins (P_N), required for RNA synthesis, and an unspecified host factor HF_RC, required for the formation of the RC, are imported into the RC. These assumptions directly match those made in our published HCV replication model (Binder et al., 2013). As the total amount of host factor HF_RC is assumed constant during the time scales considered in the model, a separate equation for HF_RC is not necessary. We furthermore note that, while numerous replication vesicles can be observed during DV replication in every single cell (Welsch et al., 2009), we assume in the model that the sum of all replication vesicles is regarded as a single, large replication compartment, compare (Dahari et al., 2007).

After formation of the replication initiation complex R_RC, minus-strand RNA synthesis is initiated with rate constant k_4m, leading to the production of dsRNA (R_DS, Equation 12) and liberation of viral proteins which remain in the RC (P_RC, Equation 15). For the synthesis of (+)RNA (R_P, Equation 16), dsRNA (R_DS) binds again to P_RC with constant rate k₅ in order to form a minus-strand RNA intermediate complex (R_IDS, Equation 13). The minus-strand RNA intermediate complex (R_IDS) serves as a template for (+)RNA synthesis with constant rate k_4p and subsequently dissociates into dsRNA (R_DS) and viral protein (P_RC). The newly synthesized (+)RNA (R_P) can then either be transported out of the RC into the cytoplasm with rate constant k_Pout, it can be used for a further round of replication with rate constant k₃, or it is used to assemble new virions, which are then released from the cell (v_p). We assume that all species in the RC are protected from active degradation, and decay together with the membrane vesicles with a common rate constant μ_RC.

Since the RCs might represent a protective environment for DV replication by shielding DV RNA from the host's immune response recognition (Scutigliani and Kikkert, 2017), we introduced a cytosolic dsRNA species (R_R, Equation 14a). Therefore, for the dsRNA species within the RC, R_DS and R_IDS, μ_RC represents a transfer rate into the cytoplasm and leads to the accumulation of cytosolic dsRNA that is detectable by the innate immune sensor (Chazal et al., 2018), while the RNA species within the RC remain protected. In order to account for a slow transfer rate without introducing another model parameter, we use the RNA degradation rate inside the replication compartment, μ_RC, for this purpose. Similar to the single stranded RNA species within the cytoplasm, the cytosolic dsRNA degrades with rate μ_RV.

Finally, Equations (17) to (19) model the assembly and release of new virus particles. To produce one infectious virion, the newly synthesized (+)RNA (R_P) assembles together with structural proteins: 180 C proteins, 180 E, and 180 prM/M proteins (Kuhn et al., 2002). Moreover, it has been shown that non-structural proteins support DV particle production as well, e.g., DV NS1 (Scaturro et al., 2015). However, we observed during model fitting that non-structural proteins were not rate limiting for virus assembly and therefore we neglected P_N in the assembly process, while we focused on structural proteins and host factors required or participating in virus assembly and release (for more information see Supplementary Material).

We model DV assembly and release (v_p) using a Michaelis-Menten type equation, as

$\begin{array}{ll}{v}_p=k_p{R}_P{\displaystyle \prod _j}\frac{P_j}{K_D·N_{P_j}+P_j},& (17)\end{array}$

with j ∈ {P_S, HF_PP}, N_Pj the number of each protein P_j, and cell line-specific virion release rate constant k_p, compare (Heldt et al., 2012). We require here that sufficient proteins per virion need to be available in order to reach the half-maximal virion release rate K_D. Furthermore, we introduced a second host factor HF_PP for particle production, with a cell line-specific basal production rate k_HFPP:

$\begin{array}{ll}\frac{dH{F}_{PP}}{dt}=k_{H{F}_{PP}}-{N_{H{F}_{PP}}v}_p.& (18)\end{array}$

The infectious virions (V_INF) are released from the cell and are then able to infect naïve cells with a reinfection rate constant k_re (Equation 19); they furthermore are assumed to degrade with rate constant μ_V, thus

$\begin{array}{ll}\frac{d{V}_{INF}}{dt}=v_p-k_{re}{V}_{INF}-{\mu }_V{V}_{INF}.& (19)\end{array}$

Host Cell Immune Response

We coupled the replication model with a simple model containing the central steps of the cell's IIR to infection. This HIR sub-model (Figure 2) comprises four main processes: (1) The recognition of viral RNA by cellular PRRs, leading to the initiation of a signaling cascade that causes (2) the production and release of IFN. (3) Subsequently, secreted IFN triggers the transcription of mRNAs of hundreds of ISG, which are then (4) translated into ISG proteins that act upon multiple processes in the DV lifecycle. To keep the model simple, we subsume the different ISGs by a single representative species, included in the model by its mRNA (I_R) and protein (I_P) species.

As mentioned above and to keep the model simple, we include only dsRNA recognition via RIG-I/MDA-5 into the model. As soon as dsRNA is detectable in the cytoplasm, it activates the HIR. We therefore modified the equation for R_R (Equation 14a) as follows (changes in bold), where k_rig describes the rate at which the RIG-I pathway is activated and IFN is produced when cytoplasmic R_R is bound by the receptor:

$\begin{array}{ll}\frac{d{R}_R}{dt}={\mu }_{RC}(R_{DS}+R_{IDS})-{\mu }_{R_V}{R}_R-k_{rig}{R}_R.(14\text{b})& (\mathrm{14b})\end{array}$

The dynamics of key components of the HIR, namely intracellular IFN (F_IN), secreted IFN (F_EX), ISG mRNA (I_R), and ISG protein (I_P) are given by the following ODEs:

$\begin{array}{ll}\frac{d{F}_{IN}}{dt}=k_{rig}{R}_R-k_s{F}_{IN}-{\mu }_F{F}_{IN},& (20)\end{array}$

$\begin{array}{ll}\frac{d{F}_{EX}}{dt}=k_s{F}_{IN}-k_{jak}{F}_{EX}-{\mu }_F{F}_{EX},& (21)\end{array}$

$\begin{array}{l}\frac{d{I}_R}{dt}=k_{jak}{F}_{ex}-k_{IC}{I}_R(Rib{o_{HIR}}_{tot}-IC)+k_{t}IC\\ -{\mu }_{I_R}{I}_R,& (22)\end{array}$

$\begin{array}{ll}\frac{dIC}{dt}=k_{IC}{I}_R(Rib{o_{HIR}}_{tot}-IC)-k_{t}IC-{\mu }_{IC}IC,& (23)\end{array}$

$\begin{array}{ll}\frac{d{I}_P}{dt}=k_{t}IC-{\mu }_{I_P}{I}_P.& (24)\end{array}$

Here, upon recognition of cytoplasmic dsRNA (R_R, Equation 14b), the cell produces IFN (F_IN, Equation 20) via the RIG-I/MDA-5 pathway with rate constant k_rig. IFN either degrades with rate constant μ_F or is secreted from the cell (F_EX, Equation 21) with rate constant k_s and then degrades extracellularly with rate constant μ_F. Extracellularly, F_EX binds to receptors that activate the JAK/STAT pathway, we assume this to happen with rate constant k_jak. Activation of the JAK/STAT pathway then triggers the production of ISG mRNAs (I_R, Equation 22), which we assume to degrade with rate constant μ_IR. Ribosomes (Ribo_HIR) bind to I_R to form a translation initiation complex (IC, Equation 23) with rate constant k_IC, which in turn degrades with rate constant μ_IC. The subsequent translation of the translation initiation complex (IC) with rate constant k_t leads to the production of ISG proteins (I_P, Equation 24). We assume that I_P degrades with rate constant μ_IP. Similar to TC, IC dissociates into Ribo_HIR and I_R (I_R + Ribo_HIR → IC → I_P + I_R + Ribo_HIR). Note that in order to prevent ribosomal competition in the model, we discriminate between ribosomes used for DV RNA translation (Ribo_DV) and ribosomes used for the HIR (Ribo_HIR).

The ISG proteins (I_P) affect numerous processes in the viral lifecycle. Here we focus on effects it has on the formation of the translation initiation complex (k₁) and the degradation of viral RNA in the cytoplasm (μ_RV) (Diamond, 2014; Castillo Ramirez and Urcuqui-Inchima, 2015). We furthermore assume that the ISGs cannot reach species inside of the replication vesicles, which thus provides a protected environment for viral replication (see Supplementary Material for details). We include these effects into the model by decreasing the corresponding reaction rate constant k₁ to

$\begin{array}{ll}\widehat{k_1}=\frac{k_1}{1+{\epsilon }_{k_1}{I}_P},& (25)\end{array}$

and increasing the degradation rate μ_RV to

$\begin{array}{ll}\widehat{{\mu }_{R_V}}={\mu }_{R_V}(1+{\epsilon }_{{\mu }_{RV}}{I}_P).& (26)\end{array}$

Furthermore, we take into account that IFN released from infected cells can protect naïve cells from infection by bringing them into an antiviral state, this has been integrated into the model by decreasing the corresponding reaction rate constant k_re to

$\begin{array}{ll}\widehat{k_{re}}=\frac{k_{re}}{1+{\epsilon }_{k_{re}}{F}_{EX}}.& (27)\end{array}$

The efficacy constants ε_* measures the efficacy of the inhibition on a range from 0 (no effect) to 1 (full inhibition).

Viral Countermeasures

DV is not only subjected to the HIR, but viral proteins in turn also impair the host's immune response, thus constituting a negative feedback loop of mutual inhibition. Several viral proteins have been described inhibiting HIR pathway activation. For example, DV NS3, NS4A, and NS5 inhibit the RIG-I pathway activation by the methylation of the DV RNA (DV NS5) or by blocking the RIG-I/MAVS interaction (DV NS4A) (Chazal et al., 2018). Additionally, by promoting the degradation of STAT2, DV NS5 impairs activation of the JAK/STAT pathway and thus ultimately inhibits ISG production upon exposure of the cell to exogenous IFN (Ashour et al., 2009; Mazzon et al., 2009). Therefore, we incorporated the ability of DV to circumvent the HIR in two ways ($\widehat{c_x}$): (i) by reducing the reaction rate of the RIG-I pathway activation that may lead to a decreased IFN production (k_rig), and (ii) by decreasing the reaction rate for the JAK/STAT pathway activation that may result in a decreased ISG expression (k_jak). Similarly, to the antiviral HIR effect, we incorporated these viral countermeasure effects into our model

$\begin{array}{ll}\widehat{c_x}=\frac{c_x}{1+{\epsilon }_x{P}_N},& (28)\end{array}$

with c_x ∈ {k_rig, k_jak} and its efficiency constant ${\epsilon }_x\in \left[10^{-5},1\right]$, dependent on DV non-structural protein concentration. Hence, we replaced the rate constants k_rig and k_jak in equations (14b), (20), (21), and (22) by the terms $\widehat{k_{rig}}$ and $\widehat{k_{jak}}$ as defined above (see Supplementary Material for details).

Parameter Estimation

We implemented the mathematical model in Matlab Release 2016b (The Mathworks). Twenty out of the total 56 model parameters were fixed based on evidence from literature, direct calculations or observations made during the optimization process, see Tables S6, S7. In brief, since infection experiments were carried out at a multiplicity of infection (MOI) of 10 and assuming that the fraction of infected cells follows a Poisson distribution (Flint et al., 2009; Wulff et al., 2012), we computed an initial virus concentration of V₀ = 10 infectious virus particles per ml per cell. Washing of cells to remove unbound virus was performed thrice after 1 h for a duration of 6 min, we therefore set ω_t = 1 h, ω_d = 0.1 h and assume a washing strength of ω_s = 100. Model parameters for translation and replication rates were estimated based on the DV genome length of approximately 10,700 nucleotides and DV polyprotein length of 3,400 amino acids. During the fitting process, we observed no difference whether the ribosomes bind to viral RNA (R_V) or host cell mRNA (I_R) and set k₁ = k_IC. Assuming a translation velocity of 3 to 8 AA/s per polysome and assuming 10 ribosomes bound to each 2,000 AA (Dahari et al., 2009), we obtain $k_2=100h^{-1}$. The translation rate (k_t) for the ISG proteins (I_P) was calculated accordingly as $k_t=120h^{-1}$, using the IFIT1 protein as a representative ISG with a length of 478 AA (Safran et al., 2010). RNA synthesis rate constants were calculated as $k_{4m}=k_{4p}=1.01h^{-1}$, using a transcription rate of 180 nt per min (Dahari et al., 2009). Degradation rates for extracellular virus and IFN were set to ${\mu }_V=0.4h^{-1}$ and ${\mu }_F=0.15h^{-1}$ (Schmid et al., 2015). Note that for simplicity, we assumed that the intracellular IFN degradation equals the extracellular degradation, μ_F. We observed a higher stability of virus within endosomes (μ_VE) compared to extracellular virus (μ_V) and fixed the degradation rate of virus within endosomes to μ_VE = 0.5·μ_V. The degradation rate for luciferase ${\mu }_L=0.35h^{-1}$ as well as the polyprotein cleavage rate $k_c=1h^{-1}$ were taken from our HCV replication model (Binder et al., 2013). The translation initial complexes TC and IC are assumed to be more stable than free cytosolic DV RNA genome (R_V) and ISG mRNA (I_R) due to the bound ribosomes. Therefore, the degradation rates μ_TC and μ_IC are assumed to be lower than the degradation rates μ_RV and μ_IR. The degradation rate for ISG proteins was fixed to ${\mu }_{I_P}=0.03h^{-1}$, corresponding to a half-life of t_1/2 = 24 h based on literature showing half-lives for ISG proteins in the range of 12 h and 2.3 days (Ronni et al., 1993; Martensen and Justesen, 2004; Haller et al., 2007; Bogunovic et al., 2013). The half-maximal virion release rate (K_D) needed for the virus assembly and release term (v_p) was approximated based on the experimental measurements of extracellular virus particles. Here, we calculated that in Huh7 cells, K_D = 1.8 virions/ml per measurement time point were produced, while in A549 cells, K_D = 0.7 virions/ml per measurement time point were produced. The number of structural proteins required to produce one virion has been taken from literature with N_PS = 180 molecules/virion (Kuhn et al., 2002). During the fitting process, we observed a 10 times higher basal production rate for the host factor involved in assembly/release in Huh7 cells than in A459 cells. We therefore set $k_{H{F}_{PP}}^{Huh7}=10·k_{H{F}_{PP}}^{A549}$. Furthermore, we observed that the initial concentration of the cell line specific host factor involved in virus assembly and release was fitted to the same value and thus set $H{F}_{P{P}_0}^{Huh7}=H{F}_{P{P}_0}^{A549}=H{F}_{P{P}_0}$. However, we assume that the initial host factor (HF_RC0, HF_PP0) and ribosome (Ribo_DV0, Ribo_HIR0) concentrations, as well as the number of consumed host factors in the virus assembly and release process (N_HFPP) are ≥1 molecules/virion. The antiviral HIR and DV countermeasure efficiency constants were estimated within ${\epsilon }_{*}\in \left[10^{-5},1\right]$, while the remaining model parameters have been estimated within the range [10⁻⁵, 10³]. Initial specie concentrations were V_A0 = V_E0 = V_INF0 = 0 virions/ml for virus species, R_V0 = TC₀ = R_RC0 = R_DS0 = R_IDS0 = R_R0 = R_P = P_P0 = P_S0 = P_N0 = L₀ = P_RC0 = 0 molecules/ml for viral RNA, protein and luciferase species and F_IN0 = F_EX0 = I_R0 = IC₀ = I_P0 = 0 molecules/ml for the IFN and ISG species, while the initial concentrations of host factors (HF_RC0 ≠ HF_PP0 ≠ 0 molecules/ml) and ribosomes (Ribo_DV0 ≠ Ribo_HIR0 ≠ 0 molecules/ml) have been estimated (for more details, see Supplementary Material).

To fit the model to the experimental data, we computed $R_P^{tot}$ = (V_E + R_V + TC + R_RC + R_DS + R_IDS + R_R + R_P), $V^{tot}=(V+V_A+V_{INF})$, and $I_R^{tot}=(I_R+IC)$ and introduced four scaling factors fScale_L, fScale_R, fScale_F, and fScale_IR to rescale experimental measurements acquired in relative values (Luciferase, DV RNA, and ISG mRNA) and pg/ml (IFN). Remaining free model parameters were then estimated from the experimental data. Parameter estimation was performed using the Data2Dynamics Matlab toolbox (Raue et al., 2015), using a deterministic trust region algorithm (lsqnonlin) with Latin hyper cube multi-start, minimizing the log likelihood function (Raue et al., 2013) (for more details see Supplementary Material). Parameter estimation was performed simultaneously for the Huh7 and A549 cell lines, where only the DV replication sub-model was used in the Huh7 cells and the full model including the immune response sub-model in the A549 cell line. The only other parameters that were allowed to vary between the two cell lines were the initial concentrations of the host factor for the formation of the minus-strand synthesis complex (HF_RC0) as well as the basal production (k_HFPP) of the host factor for particle production (HF_PP0). It is likely that more processes are cell line specific, however, here we summarized all model parameters that did not show any impact on the model fit and focused mainly on cell line specific host factor availability and HIR effects. Tables S6, S7 summarize the final, resulting model parameters used after model fitting.

Model Analysis

Simulation of Antiviral Intervention

We used the model to study potential antiviral strategies. For this purpose, we extended the model by effects of direct-acting antiviral drugs (DAAs). A drug efficacy parameter 0 ≤ ε ≤ 1 was introduced to simulate drug effects on viral attachment (k_a), viral entry (k_e), formation of the translation initiation complex (k₁), translation (k₂), polyprotein cleavage (k_c), replication complex formation (k_Pin), minus- and (+)RNA synthesis (k_4m and k_4p), virus particle production and release (v_p), and infection of naïve cells (k_re), by multiplying the corresponding parameter with (1 − ε). We then calculated the time-averaged infectious virus particle concentration released from the cell upon drug administration (ε > 0) within a time interval of 5 days (Δt = 120 h), normalized to the time-averaged infectious virus concentration without drug treatment (ε = 0) as

$\begin{array}{ll}\psi =\frac{\langle V_{INF}{(t)}_{T_{\epsilon >0}}\rangle }{\langle V_{INF}{(t)}_{T_{\epsilon =0}}\rangle },& (29)\end{array}$

with

$\begin{array}{ll}\langle V_{INF}{(t)}_T\rangle =\frac{1}{\Delta t}{\displaystyle {\int }_t^{T+\Delta t}}dt{V}_{INF}(t),& (30)\end{array}$

where T refers to the time point of drug administration (T ≤ t ≤ T + Δt).

Identifiability and Sensitivity Analysis

We assessed model identifiability using the profile likelihood method, which analyzes both structural and practical identifiability. The profile likelihood method evaluates the change in the likelihood function after modification of one individual model parameter by re-optimizing the remaining model parameters (Raue et al., 2009; Kreutz et al., 2013; Maiwald et al., 2016), thus assessing if changes in a given parameter can be compensated by modifications in other model parameters. Based on the profile likelihood, we calculated 95% confidence intervals on model parameters, which imply parameter identifiability if the confidence interval is finite (Raue et al., 2009, 2015). Local and global sensitivity analysis were carried out in Matlab using the SensSB toolbox (Rodriguez-Fernandez and Banga, 2010) and the extended Fourier Amplitude Sensitivity Test (eFAST) (Marino et al., 2008). Sensitivities with regard to polyprotein (P_P), total (+)RNA ($R_P^{tot}$) and total Virus (V^tot) concentrations were calculated for the two time points t = 4 hpi (early during infection) and t = 72 hpi (at steady state).

Results

In order to in silico analyze the full DV lifecycle in the absence and presence of HIR mechanisms, we developed a detailed model of the intracellular DV lifecycle and coupled this model to an HIR model, taking into account the antiviral effect of an active immune response on the DV lifecycle as well as DV's ability to in return attenuate the HIR (Figures 1, 2). Model calibration was performed by estimating model parameters simultaneously on experimental data measured in two different cell lines, Huh7 and A549 cells (for details see Materials and Methods). For this purpose, we measured viral polyprotein (luciferase readout), (+)RNA, extracellular virus, and IFN in both cell lines, while ISG mRNA was measured only in A549 cells, as the Huh7 cell do not show activation of the interferon system after DV infection. Polyprotein (luciferase) showed a 1-log₁₀ higher translation activity in Huh7 cells compared to A549 cells (Figure 3A). Similarly, Huh7 cells showed a higher extracellular infectious virus concentration compared to A549 cells (Figure 3C). However, in both cell lines, the extracellular virus concentration drops after reaching a peak (~32 hpi in Huh7 and 36 hpi in A549 cells). Nevertheless, against our expectations, DV (+)RNA measurements showed the opposite trend with a faster RNA production in A549 cells (Figure 3B). Additionally, IFN has been measured in both cell lines and showed an increase in secreted IFN in A549 cells (which is followed by ISG mRNA) and a baseline IFN level in Huh7 cells (Figures 3D,E).

FIGURE 3

Figure 3. Best model fit of the DV replication model (Huh7) and the coupled model (A549) compared to experimental data measured in Huh7 and A549 cells (for parameter values see Tables S6, S7). Experimental measurements are represented as mean μ± standard deviation σ. The DV replication model and the coupled model were fitted simultaneously to the Huh7 and A549 data sets with cell line specific differences mediated by host factors, the antiviral HIR effect, and DV countermeasures (Equations 1 to 28). (A) shows the model fit of luciferase compared to the luciferase measurements (L = Luc), (B) model fit of total (+)RNA to the (+)RNA measurements ($R_P^{tot}$ = (+)RNA), (C) model fit of extracellular virus to its measurements (V^tot= Virus), (D) model fit compared to measurements of the HIR ($I_R^{tot}=$ ISG mRNA and F_EX = extracellular IFN), while (E) shows the model prediction of the coupled model with cell line specific Huh7 parameter values and the knocked-out RIG-I pathway activation $k_{rig}=0h^{-1}$ compared to IFN measured in Huh7 cells.

Our coupled model, developed based on best biological knowledge, showed high agreement with our experimental cell-line specific measurements after fitting model parameters to the data (Figure 3). Due to the high degree of freedom of the model and in order to prevent overfitting, we analyzed structural and practical identifiability of model parameters. Results are shown in Figure S12; as can be seen in the figure, due to the high model complexity, not all model parameters are identifiable. In particular the parameters for replication within the RC (k₃) and DV RNA export out of the RC (k_Pout) are non-identifiable. Both parameters concern the use of newly synthesized plus-strand RNA and reflect the allocation of such newly synthesized RNA to either further rounds of RNA replication (k₃) or to export from the replication compartment and use for protein translation. The fact that these two parameters are non-identifiable is surprising at first, as allocation of newly synthesized RNA between these processes should significantly affect viral replication dynamics. However, this can be explained by other processes that are rate-limiting. In fact, we observed a similar behavior in our HCV replication model (Binder et al., 2013), where in high permissive cell lines the HCV RNA-dependent RNA polymerase became in fact rate limiting for RNA replication inside of the replication vesicles, and led to a similar insensitivity of the model to the parameter k₃. This is reflected in DV as well, with only limited impact of parameters k₃ and k_Pout on the replication dynamics in both cell lines.

Host Dependency and Restriction Factors

The first question we addressed was cell line specificity. In contrast to our expectations, we observed a faster onset and more efficient viral replication in the HIR-competent A549 cells. Here, our model was not able to describe the DV RNA dynamics in A549 cells that seemed unaffected by the HIR and showed a faster increase and an overall 2.7-fold higher amount of DV RNA (time-averaged concentration) compared to Huh7 cells. Viruses strongly depend on their living hosts and hijack host cell membranes, proteins, and lipids for their own replication. We thus speculated that other host processes explain this difference between the two cell lines. We tested different such potential host factors by including corresponding cell-line specific parameters into the model, compare Table S3, and discriminated between these models using model selection based on Aikaike's Information Criterion (AIC). In line with our previous findings with our published HCV model (Binder et al., 2013), we introduced an unspecific cell line specific host factor, HF_RC, that participates in the assembly of the replication complex and RC formation. In model fitting, this host factor showed a higher availability leading to a faster onset of DV replication in A549 cells compared to Huh7 cells (Table S6). We furthermore tested cell-line specific effects on different host factor supported processes such as virus entry, that improved the overall model fit without explaining the cell line specific DV RNA dynamics (see Table S3).

We additionally observed a cell line specific variation in the extracellular DV dynamics, resulting in a 2.8-fold lower extracellular virus concentration (time-averaged) in A549 cells, that could not be described by the HIR alone. Thus, we introduced another unspecific host factor, HF_PP, involved in virus particle production and release, with a cell line specific basal production, k_HFPP, and a cell line specific virus assembly and release rate, k_p. During model parameter estimation, we observed a faster virus assembly and release and an around 10 times faster basal production of the host factor involved in DV assembly and release in Huh7 cells compared to A549 cells. This basal host factor production was the key parameter for the lower virus concentration in A549 cells in steady state. Furthermore, this host factor represented a limiting species for particle production and release, since the drop in the extracellular DV concentration following the peak was associated with a drop in the host factor concentration. Interestingly, the availability of structural proteins had no effect on the drop in released virus (see Figure S1).

The Antiviral HIR Effect and Viral Countermeasures

During an acute infection, the IIR represents the first line of defense against an invading pathogen. The IIR is mounted by the production of IFN and subsequent ISG translation; the ISGs in turn subsequently inhibit multiple steps in the viral lifecycle. Having developed a detailed model coupling the DV lifecycle with key players of the HIR, we studied the antiviral modes of action in detail and introduced three possible antiviral HIR effects on (i) the translation initiation (k₁), (ii) the degradation of free cytosolic DV RNA (μ_RV), and (iii) the reinfection of naïve cells (k_re) into the model. Selection of these three main mechanisms was based on model selection using the least squared error with the AIC to account for model complexity. For details we refer to the Supplementary Material. By comparing the model fits and its AICs, we observed the best model fit and lowest AIC for a model in which the HIR inhibits the translation initiation complex formation (k₁), followed by a model, that increases the degradation of free cytosolic viral RNA (μ_RV). However, the model considering only the increase in the cytosolic RNA degradation (μ_RV) resulted in a very high cytosolic DV RNA degradation rate of $\widehat{{\mu }_{R_V}}=987h^{-1}$. The model that led to an antiviral state by inhibiting reinfection of naïve cell infection led to the third best model. Since we are interested in combinatory effects, we chose all three antiviral ISG and IFN dependent effects as our working model.

In the combined HIR effect model, the inhibition of the translation initiation (k₁) and the reinfection of susceptible cells (k_re) by the HIR showed the highest efficacy constants with ε_k1 = ε_kre = 1 (Table S7). Comparing the inhibitory effect on the effective rate constants over time, the rate constant for k₁ dropped by 99.983% from its initial value, while k_re showed a negligible 1.6% decrease (Figure 4 and Table 1). The cytoplasmic RNA degradation rate was increased by 58.7%. DV has the ability to evade the HIR by decreasing or inhibiting its own recognition, correspondingly, the rate constants for RIG-I pathway activation and dsRNA recognition was reduced by 93.6%, while the JAK/STAT pathway activation was reduced by 88.6% (Table 1).

FIGURE 4

Figure 4. Percentage increase and decrease of model parameters of the antiviral ($\widehat{k_1},\widehat{k_{re}},\widehat{{\mu }_{R_V}}$) and countermeasure ($\widehat{k_{rig}},\widehat{k_{jak}}$) effect of the HIR and DV on model parameters as a function of time (Equations 25 to 28; Table 1). The antiviral HIR mediated inhibition of the (A) translation initiation complex formation k₁, (B) naïve cell infection/reinfection k_re, (C) cytosolic DV RNA degradation. The DV mediated countermeasure of (D) RIG-I pathway k_rig and (E) JAK/STAT pathway k_jak, as well as (F) the ISG protein concentration over time. Note that $\widehat{k_1}$ and $\widehat{{\mu }_{R_V}}$ are ISG dependent, while $\widehat{k_{re}}$ is IFN dependent.

TABLE 1

Table 1. Effect of the immune response on DV replication parameters and of Dengue on parameters of the immune pathways—change in parameter values over 100 h for HIR affected processes in the DV lifecycle and HIR pathways that are targeted by DV: Translation initiation complex formation (k₁), naïve cell infection/reinfection (k_re), cytosolic RNA degradation (μ_RV), RIG-I pathway (k_rig), and JAK/STAT pathway (k_jak).

In order to further mathematically analyze the interplay of antiviral effects (k_rig and k_jak) and the viral ability to attenuate the HIR (ε_krig and ε_kjak), we performed a bifurcation analysis at time 72 hpi. Here, we compared the (+)RNA concentration to various effect combinations: (i) the recognition of DV dsRNA (k_rig) vs. DV's ability to attenuate its own recognition (ε_krig) and (ii) the activation of the JAK/STAT pathway (k_jak) vs. DV countermeasures targeting the JAK/STAT pathway (ε_kjak) (Figure 5). In the first scenario, we observed a clear separation: with increasing k_rig the HIR wins and the infection is effectively cleared with only minimal residual (+)RNA, while with increasing ε_krig the virus wins and the infection is ongoing. In contrast, in the second scenario we found that increasing the activation of the JAK/STAT signaling pathway (k_jak) did not lead to significant decreases in viral RNA levels.

FIGURE 5

Figure 5. Plus-strand RNA concentration for various model parameter combinations for: (A) the antiviral RIG-I pathway activation (k_rig) vs. the viral countermeasure targeting the RIG-I pathway (ε_krig) and (B) the JAK/STAT pathway activation (k_jak) vs. its viral counteract (ε_kjak). The black lines represent the model parameter combinations that have been estimated from the data (Table S7).

Antiviral Drug Intervention

We were further interested in using the mathematical model to identify processes with a high impact on the DV lifecycle, as those might constitute attractive targets for antiviral drug development. For this purpose, we performed a global sensitivity analysis to analyze the effect of all model parameters on viral polyprotein, total DV (+)RNA and extracellular virus concentrations at two distinct time points, 4 and 72 hpi (Figures 6, 7).

FIGURE 6

Figure 6. Global Sensitivity analysis performed for the DV replication model for polyprotein (A,B), (+)RNA (C,D), and extracellular virus (E,F), as well as for two different time points: 4 hpi (A,C,E) and 72 hpi (B,D,F). The red line is a negative control used for the sensitivity analysis that is not part of the model indicating that sensitivities below the red line are negligible. Significant differences to the negative control have been calculated by performing a t-test (p-values: *** ≤ 0.001, ** ≤ 0.01, * ≤ 0.05).

FIGURE 7

Figure 7. Global Sensitivity analysis performed for the DV replication model coupled to the HIR model for polyprotein (A,B), (+)RNA (C,D), and extracellular virus (E,F), as well as for two different time points: 4 hpi (A,C,E) and 72 hpi (B,D,F). The red line is a negative control used for the sensitivity analysis that is not part of the model indicating that sensitivities below the red line are negligible. Significant differences to the negative control have been calculated by performing a t-test (p-values: *** ≤ 0.001, ** ≤ 0.01, * ≤ 0.05).

Both cell lines showed a comparable sensitivity profile for polyprotein, DV (+)RNA, and the extracellular virus and showed high sensitivities to processes associated with cell infection, polyprotein translation and processing, and DV RNA synthesis. For all the three species, early processes in the viral lifecycle were associated with highly significant sensitivities at 4 hpi, such as virus attachment (k_a), entry (k_e), and fusion (k_f), as well as polyprotein translation (k₂). Later in infection (72 hpi), ongoing polyprotein translation as well as processes within the RC dominated in their sensitivities for the three studied species. Especially polyprotein cleavage (k_c) became the dominant process with the highest impact of all steps involved in viral protein translation. Of the processes occurring inside of the RC, the most sensitive rate constants were associated with RNA synthesis (k_4m, k_4p) and the host factor involved in the formation of the RC (HF_RC) for both cell lines. For the extracellular virus, the number of host factors (N_HFPP) involved in assembly and release showed a higher sensitivity compared to the number of viral structural proteins (N_PS). Amongst the HIR model parameters, the RIG-I pathway activation (k_rig) showed a slightly higher, significant sensitivity on the polyprotein species. Furthermore, the HIR efficacy constant decreasing the rate constant of the naïve cell infection (ε_μRV) showed the highest sensitivity of all antiviral HIR constants, albeit not reaching statistical significance.

As a next step, we were interested in the question whether the highly sensitive processes identified in the previous analysis might represent potent drug targets. We therefore performed a theoretical antiviral intervention by simulating a possible drug administration. In this simulation, we monitored the release of infectious virus for 5 days following drug administration. Several processes in the DV lifecycle were inhibited by simulated drug administration at 0 hpi, 24 hpi, and 72 hpi (Figure 8). An early drug administration at 0 hpi led to an efficient viral clearance in both cell lines, using a hypothetical drug acting on any process in the DV lifecycle except for putative drugs targeting reinfection. With the support from the HIR, the overall drug efficacy constants necessary to eradicate the virus in A549 cells were lower. In particular drugs targeting translation initiation and the DV RNA synthesis were able to induce viral clearance even with low drug efficacy constants, and administering a drug targeting the DV RNA synthesis process showed a viral clearance with the lowest drug efficacy constant (ε ≈ 0.5) in A549 cells. For drugs targeting any one of the remaining processes, drug efficacy constants higher than ε ≥ 0.9 were needed to clear the viral infection. Administering a hypothetical drug at 24 and 72 hpi led to comparable viral clearance patterns, but with higher drug efficacy constants. Obviously, if a drug is administered late in the viral lifecycle and targets early processes of the viral lifecycle such as virus attachment, endocytosis and fusion as well as formation of the (membranous) replication compartment, leads to a loss of the drug effect and non-clearance of the DV infection in both cell lines. In both cell lines, the DV infection can still be cleared when blocking DV RNA synthesis and virus assembly/release with <3% DV left with the highest drug efficacy of 1 (thus completely shutting of RNA synthesis and assembly/release), an outcome which cannot be achieved by targeting any of the other processes in the DV lifecycle according to our model simulations.

FIGURE 8

Figure 8. Antiviral intervention study with a drug administration at three different time points (0, 24, and 72 hpi) in (A) A549 and (B) Huh7 cells. The fold change of infectious virus (ψ) with and without drug administration for the core processes in the DV lifecycle (Equations 29 and 30). A fold change of 1 means no difference between the model with and without drug administration, while a fold change of 0 shows viral eradication to a successful drug treatment.

Discussion

In the present study, we investigated the intracellular virus replication and HIR dynamics in two different cell lines: Huh7 cells with very low HIR-competence and highly HIR-competent A549 cells. Several cell population models have been developed to analyze DV dynamics under influence of the innate (Schmid et al., 2015) and the adaptive immune response (Ben-Shachar and Koelle, 2014; Clapham et al., 2014, 2016; Ben-Shachar et al., 2016). These models, however, do not take intracellular processes into account and thus lack molecular detail. In order to study the intracellular dynamics during the DV lifecycle, we developed the first mathematical model that reflects the initial dynamics of key components of the intracellular DV genome replication. Our detailed model is derived from previous mathematical models that have been used to describe the intracellular RNA replication of a HCV replicon system after RNA transfection (Dahari et al., 2007; Binder et al., 2013). We extended these models by including virus entry and release of infectious virus particles. Furthermore, we coupled the virus dynamics model to a model of the HIR activation and effector phases in order to study the modes of action of the HIR, and to analyze potential antiviral intervention strategies acting at the level of intracellular mechanisms.

Host Factors

Our experimental measurements were performed in two different cell lines: Huh7 cells which show no interferon response, and A549 cells showing a strong immune reaction. However, Huh7 cells are based on hepatoma (liver) cells, whereas A549 cells are of pulmonary epithelial origin, thus they likely differ in several other aspects as well. In fact, some characteristics of our experimental data cannot be explained by the lack of an interferon response in the Huh7 cells alone. Contrary to our expectations, we observed a faster onset and more efficient DV genome replication in the immuno-competent A549 cells. We therefore tested which other host factors may explain such cell line specific differences. We set up several different models for such host factors, fitted the corresponding models to the experimental data, and compared different models using AIC; details are given in the Supplementary Material. In our previous HCV study, we have shown that host factors involved in replication complex formation play a crucial role in cell permissiveness and viral replication efficiency (Binder et al., 2013). Similarly, for DV, such a host factor best explained differences in replication efficiency between the two cell lines. According to our model, the more efficient RNA replication (earlier increase and a higher steady state concentration of total (+)RNA in the A549 cell line) is directly associated with a higher concentration of this putative host factor in A549 cells, similar to our previous results considering HCV replication in different Huh7 cell clones (Binder et al., 2013).

Concerning the extracellular virus dynamics, our model was not able to explain the drop in infectious virus titers observed in the experimental data after ~40 h post infection—at different degrees in both cell lines—by a limitation in structural viral proteins. In fact, our simulation results show that structural proteins do not limit the process of particle production and release. Similar to our finding, Heldt et al. (2012) in a mathematical model of influenza A virus replication did not find a limitation in structural proteins and suggested that transport and budding processes might limit the viral production. Furthermore, the drop in virus titers that we observed in our data is qualitatively present in both cell lines, i.e., in the presence and in the absence of the HIR, it is therefore unlikely that it is due to effects of the HIR on virus assembly and release. Therefore, we integrated another unspecific host factor, HF_PP, that is involved in virion assembly, maturation, and release into the mathematical model with a cell line-specific basal host factor production and a cell line specific virus assembly and release rate. Fitting of this extended model resulted in a higher production rate of this assembly/release host factor in Huh7 cells, explaining the higher viral steady state level in these cells. Several host factors affecting DV assembly / release are known; we recently employed siRNA screening to identify such factors and described a mechanisms involving Fibroblast Growth Factor Receptor 4 (FGFR4), a host factor supporting DV RNA replication when FGFR4 concentration is high, but leading to increased assembly and maturation of virus particles when this host factor is depleted (Cortese et al., 2019). While FGFR4 is one potential mechanism, the exact identity and mechanisms of host factors differences between A459 and Huh7 based cell lines needs more exploration in the future.

DV and the HIR

We next employed our mathematical model to characterize the interplay between virus replication on the on hand and the HIR on the other. During DV infection, activation of the interferon system leads to the transcription of hundreds of antiviral ISG proteins at different time points, with effects on multiple steps in the viral lifecycle. In case of HCV, which is closely related to DV and one of the best studied (+)RNA viruses, ISG proteins have been identified to act on almost every step in the HCV replication cycle (Schoggins and Rice, 2011; Metz et al., 2013; Gokhale et al., 2014). Integrating such a multitude of mechanisms into a mathematical model is therefore a daunting task. To keep things simple, we tested different potential antiviral ISG mechanisms individually by including them into our mathematical model and retained the combination of mechanisms leading to the lowest AIC values in model comparison. As we assumed that the intracellular RC protects the newly synthesized viral RNA from detection by and effector mechanisms of the HIR, we did not include any ISG effects on species inside of the RC in our model (Welsch et al., 2009; Belov and van Kuppeveld, 2012; Romero-Brey et al., 2012; Cortese et al., 2017).

According to our single effect models, ISGs inhibiting RNA translation initiation and/or promoting the cytoplasmic RNA degradation led to best fits of the experimental data. However, this model resulted in a 98,600% increase in the degradation rate constant with $\widehat{{\mu }_{R_V}}=987h^{-1}$, corresponding to an unrealistic RNA half-life of t_1/2 ≈ 2 sec. A model including only IFN dependent inhibition of the reinfection of naïve cells (promoting an antiviral state in susceptible cells) was not able to reproduce the experimental data.

A combination of mechanisms based on model selection criteria described above resulted in a model including ISG effects on (1) translation initiation, (2) cytosolic RNA degradation, and (3) new infection of naïve cells. In this model, DV RNA degradation was increased by 59%, resulting in a degradation rate and half-life of $\widehat{{\mu }_{R_V}}=4.4h^{-1}$ and t_1/2 = 9 min, respectively. Concerning the reinfection of naïve cells, we observed an inhibition of about 2%, which was rather negligible. Since cells were infected with a high MOI in our experiments, i.e., virtually every cell is infected, viral spread and infection of naïve cells play only a minor role in our experimental data.

While DV is subject to ISG effects, it has also developed several strategies to evade the antiviral HIR by antagonizing and inhibiting the induction of the HIR and the antiviral state induced by it. Several DV NS proteins have been described as highly potent inhibitors of IFN signaling and production. For example, DV NS4B protein has been shown to inhibit STAT1 phosphorylation (Munoz-Jordan et al., 2003), while the DV NS5 protein is well-known to degrade STAT2 and thus result in an inhibition of type I IFN signaling (Ashour et al., 2009). According to our model simulations, during the course of infection, DV inhibits both phases of the HIR, the RLR-mediated induction of IFN by ~94%, as well as IFN signaling through the JAK/STAT signaling pathway by ~ 89%. However, in our model sensitivity analysis at 72 hpi, we found that inhibition of the JAK/STAT signaling pathway may be the more important viral defense mechanism: Increasing the efficiency of the JAK/STAT pathway activation in a sensitivity analysis did not lead to viral eradication, but still resulted in ongoing viral replication with a constant viral RNA concentration of 73%, indicating that DV efficiently counteracts activation of this pathway. In fact, DV's ability to efficiently counteract the JAK/STAT pathway has been confirmed experimentally (Muñoz-Jordán and Fredericksen, 2010). In contrast, we found that increasing the efficiency of DV recognition by the RIG-I pathway led to viral replication at a significantly lower level of only 11% remaining DV RNA. DV's ability to target the JAK/STAT pathway and thus prevent the establishing of an antiviral cellular state mediated by IFN therefore is an efficient and important viral survival mechanism.

Comparison to HCV

DV and HCV are both (+)RNA viruses of the family Flaviviridae and share key features in their lifecycles, but there are striking differences in their clinical manifestation. While a primary dengue infection is acute and occasionally associated with severe complications (DHF, DSS) but does not lead to chronic infection, the rather asymptomatic acute hepatitis C infection may develop into lifelong chronic hepatitis C with life threatening secondary manifestations, such as liver cirrhosis or hepatocellular carcinoma without successful treatment.

Comparing the dynamics of our DV model with our previous HCV model (Binder et al., 2013), we observed that the overall dynamics of luciferase and total (+)RNA in DV is comparable with the HCV dynamics with a highly dynamic initial phase that results in steady states for all the measured species. Most estimated model parameters involved in DV replication showed higher rate constants in DV compared to HCV (Table 2). Considering that DV is causing an acute infection, the faster DV replication seems reasonable, while HCV that may develop into chronicity is in comparison rather slow in its lifecycle.

TABLE 2

Table 2. Comparison of DV and HCV model parameters.

Hypothetical Drug Therapy Against DENV

The recent Zika outbreak in Brazil showed the potential health risks of (re-)emerging viruses. Therefore, a comprehensive understanding of the virus-host interaction is necessary in order to suggest antiviral treatment strategies. According to our global sensitivity analysis and simulated antiviral interventions, the most effective drug targets in the DV lifecycle are processes associated with viral entry, translation, and DV RNA synthesis. However, a drug administration earlier than 24 hpi is highly unrealistic, since dengue symptoms usually start 4 to 7 days following a mosquito bite and last for 3 to 10 days (CDC, 2014). However, targeting viral entry is suggested to prevent viral spread. Later in infection, processes associated with DV RNA synthesis and virus assembly and release still represented the most promising drug targets. The antiviral effect on post-translational and early RNA synthesis proposed by our antiviral drug intervention study might be achievable by drugs like Bromocriptine, which has shown antiviral effects against all DV-Serotypes (Kato et al., 2016). In combination with inhibitors of the DV RNA-dependent RNA polymerase, effective antiviral treatment strategies may be possible. Since all processes in the DV lifecycle depend on host factors, a future antiviral therapy may focus on host factor-targeting with the development of pan-serotype or even pan-viral antiviral drugs. As an example, the global sensitivity analysis of our model showed a high impact of the host factors involved in RC formation on DV RNA and assembly/release. To this end, the inclusion of further host factors in viral replication models might be an important challenge for future, in-silico based design of anti-DV treatment strategies.

Limitations and Outlook

In the current study, we developed the first detailed mathematical model of the intracellular DV lifecycle, coupling viral entry, protein translation, RNA replication and assembly and release with a model of the host cell immune response to infection. It has been shown in literature that stochastic effects play an important role in the activation of the IIR and individual cells in a population respond differently (Rand et al., 2012). Schmid et al. (2015) have shown that on a single cell level the IFN response to DV is stochastic and leaky with a fraction of remaining unprotected cells, in which DV replication is ongoing, emphasizing the complex nature of the IIR and virus-host interactions (Schmid et al., 2015). However, we here studied intracellular processes following DV infection in a single, “average” cell, and thus we do not take into account inter-cell differences. Since cells were infected with a high MOI, where virtually every cell is infected, viral spread is negligible and therefore, the impact of IFN released from infected cells to render non-infected, IFN-exposed cells non-permissive to DV infection is not relevant in our experimental data. We furthermore neglected cell proliferation in our model, which would require a multi-scale model combining effects at the cell population scale with a detailed intracellular model. Overall, we model an average response of an infected cell in order to study the DV lifecycle in absence and presence of the HIR, identifying HIR modes of action and sensitive processes, which might represent suitable targets for antiviral treatment.

In order to keep the HIR sub-model tractable, we simplified the activation of the HIR and took only key players of the HIR into account. Here, we model the recognition of dsRNA that is present in the cytoplasm, assuming that the replication vesicles represent a protective environment in which no RNA recognition occurs. We thus assume that DV replication intermediates are subject of detection, either when leaked into the cytosol through the replication vesicle pore or by replication vesicle decay. However, other cytosolic DV RNA species might be recognized as well, such as highly structured RNA regions in the single-stranded genome.

Furthermore, following the HIR activation, we subsume the different ISG proteins by a single species. This is a simplification that we make to keep the model simple. It is known that different ISGs are active at different time points (Metz et al., 2012), even after uniform IFN treatment (Schmid et al., 2015), hence, we here model an “average” effect. However, with our coupled model, we set the basis to study the DV-host interaction. Modeling the IIR in detail, possibly even coupling it to the adaptive immune response is needed in order to better understand and prevent severe dengue complications and to evaluate treatment strategies that suppress high-level viremia.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

LK, MB, and RB contributed conception and design of the study. BS and AR performed experiments. CZ and LK developed the model. CZ implemented and analyzed the model and data and wrote the first draft of the manuscript. CZ, AR, MB, RB, and LK wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

LK, CZ, MB, and RB received funding from the German BMBF, grant numbers 031A602A, B and C (EraSysApp/SysVirDrug). LK acknowledges funding from the BMBF, grant number 031L0032 (LiSyM). Parts of this work were done under the auspices of the U.S. Department of Energy under contract 89233218CNA000001 and supported by NIH grant R01-AI078881 awarded to AP.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.00725/full#supplementary-material

Abbreviations

AIC, Aikaike's information criterion; (+)RNA, Positive-sense RNA; DAA, Direct acting antiviral; DC, Dendritic cell; DF, Dengue fever; DHF, Dengue haemorrhagic fever; dsRNA, Double-stranded RNA; DSS, Dengue shock syndrome; DV, Dengue virus; rER, Rough Endoplasmic reticulum; HIR, Host cell immune response; HCV, Hepatitis C virus; HF, Host factor; hpi, hours post infection; IFN, Interferon; IIR, Innate immune response; ISG, Interferon stimulated gene; JAK, Janus kinase; MDA-5, Melanoma differentiation associated gene-5; MOI, Multiplicity of infection; NS, Non-structural protein; ODE, Ordinary differential equations; PRR, Pattern recognition receptor; RC, Replication compartment; RIG-I, Retinoic acid inducible gene-I; RLR, RIG-I like receptors; STAT, Signal transducer and activator of transcription factor; TLR, Toll-like receptors; WHO, World health organization.

References

Ashour, J., Laurent-Rolle, M., Shi, P-Y., and García-Sastre, A. (2009). NS5 of dengue virus mediates STAT2 binding and degradation. J. Virol. 83, 5408–5418. doi: 10.1128/JVI.02188-08

PubMed Abstract | CrossRef Full Text | Google Scholar

Aunins, T. R., Marsh, K. A., Subramanya, G., Uprichard, S. L., Perelson, A. S., and Chatterjee, A. (2018). Intracellular hepatitis C modeling predicts infection dynamics and viral protein mechanisms. J. Virol. 92:JVI.02098-7. doi: 10.1128/JVI.02098-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartenschlager, R., and Miller, S. (2008). Molecular aspects of dengue virus replication. Future Microbiol. 3, 155–165. doi: 10.2217/17460913.3.2.155

PubMed Abstract | CrossRef Full Text | Google Scholar

Belov, G. A., and van Kuppeveld, F. J. (2012). (+)RNA viruses rewire cellular pathways to build replication organelles. Curr. Opin. Virol. 2, 740–747. doi: 10.1016/j.coviro.2012.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-Shachar, R., and Koelle, K. (2014). Minimal within-host dengue models highlight the specific roles of the immune response in primary and secondary dengue infections. J. R. Soc. Interface 12:20140886. doi: 10.1098/rsif.2014.0886

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-Shachar, R., Schmidler, S., and Koelle, K. (2016). Drivers of inter-individual variation in dengue viral load dynamics. PLoS Comput. Biol. 12:e1005194. doi: 10.1371/journal.pcbi.1005194

PubMed Abstract | CrossRef Full Text | Google Scholar

Benzine, T., Brandt, R., Lovell, W. C., Yamane, D., Neddermann, P., De Francesco, R., et al. (2017). NS5A inhibitors unmask differences in functional replicase complex half-life between different hepatitis C virus strains. PLOS Pathog. 13:e1006343. doi: 10.1371/journal.ppat.1006343

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhatt, S., Gething, P. W., Brady, O. J., Messina, J. P., Farlow, A. W., Moyes, C. L., et al. (2013). The global distribution and burden of dengue. Nature 496, 504–507. doi: 10.1038/nature12060

PubMed Abstract | CrossRef Full Text | Google Scholar

Binder, M., Sulaimanov, N., Clausznitzer, D., Schulze, M., Hüber, C. M. C. M., Lenz, S. M., et al. (2013). Replication vesicles are load- and choke-points in the hepatitis C virus lifecycle. PLoS Pathog. 9:e1003561. doi: 10.1371/journal.ppat.1003561

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogunovic, D., Boisson-Dupuis, S., and Casanova, J-L. (2013). ISG15: leading a double life as a secreted molecule. Exp. Mol. Med. 45:e18. doi: 10.1038/emm.2013.36

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonhoeffer, S., May, R. M., Shaw, G. M., and Nowak, M. A. (1997). Virus dynamics and drug therapy. Proc. Natl. Acad. Sci. U.S A. 94, 6971–6976. doi: 10.1073/pnas.94.13.6971

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, L. P., Luther, C., Moo-Llanes, D., Ramsey, J. M., Danis-Lozano, R., and Peterson, A. T. (2015). Climate change influences on global distributions of dengue and chikungunya virus vectors. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 370:20140135. doi: 10.1098/rstb.2014.0135

PubMed Abstract | CrossRef Full Text | Google Scholar

Castillo Ramirez, J. A., and Urcuqui-Inchima, S. (2015). Dengue virus control of type I IFN responses: a history of manipulation and control. J. Interf. Cytokine Res. 35, 421–430. doi: 10.1089/jir.2014.0129

PubMed Abstract | CrossRef Full Text | Google Scholar

CDC (2014). Epidemiology Dengue. CDC.

Chazal, M., Beauclair, G., Golè Ne Gracias, S., Rie Najburg, V., Simon-Loriè, E., Dé, F., et al. (2018). RIG-I recognizes the 5′ region of dengue and zika virus genomes. Cell Rep. 24, 320–328. doi: 10.1016/j.celrep.2018.06.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Clapham, H. E., Quyen, T. H., Kien, D. T. H., Dorigatti, I., Simmons, C. P., Ferguson, N. M., et al. (2016). Modelling virus and antibody dynamics during dengue virus infection suggests a role for antibody in virus clearance. PLoS Comput. Biol. 12:e1004951. doi: 10.1371/journal.pcbi.1004951

PubMed Abstract | CrossRef Full Text | Google Scholar

Clapham, H. E., Tricou, V., Van Vinh Chau, N., Simmons, C. P., and Ferguson, N. M. (2014). Within-host viral dynamics of dengue serotype 1 infection. J. R. Soc. Interface 11, 504–507. doi: 10.1098/rsif.2014.0094

PubMed Abstract | CrossRef Full Text | Google Scholar

Clausznitzer, D., Harnisch, J., and Kaderali, L. (2015). Multi-scale model for hepatitis C viral load kinetics under treatment with direct acting antivirals. Virus Res. 218, 96–101. doi: 10.1016/j.virusres.2015.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cortese, M., Goellner, S., Acosta, E. G., Chatel-Chaix, L., Ruggieri, A., and Correspondence, R. B. (2017). Ultrastructural characterization of zika virus replication factories. Cell Rep. 18, 2113–2123. doi: 10.1016/j.celrep.2017.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Cortese, M., Kumar, A., Matula, P., Kaderali, L., Scaturro, P., Erfle, H., et al. (2019). Reciprocal effects of fibroblast growth factor receptor signaling on dengue virus replication and virion production. Cell Rep. 27, 2579–2592.e6. doi: 10.1016/j.celrep.2019.04.105

PubMed Abstract | CrossRef Full Text | Google Scholar

Dahari, H., Ribeiro, R. M., Rice, C. M., and Perelson, A. S. (2007). Mathematical modeling of subgenomic hepatitis C virus replication in Huh-7 cells. J. Virol. 81, 750–760. doi: 10.1128/JVI.01304-06

PubMed Abstract | CrossRef Full Text | Google Scholar

Dahari, H., Sainz, B., Perelson, A. S., and Uprichard, S. L. (2009). Modeling subgenomic hepatitis C virus RNA kinetics during treatment with alpha interferon. J. Virol. 83, 6383–6390. doi: 10.1128/JVI.02612-08

PubMed Abstract | CrossRef Full Text | Google Scholar

Dalrymple, N. A., Cimica, V., and Mackow, E. R. (2015). Dengue virus NS proteins inhibit RIG-I/MAVS signaling by blocking TBK1/IRF3 phosphorylation: dengue virus serotype 1 NS4A is a unique interferon-regulating virulence determinant. MBio 6, e00553–e00515. doi: 10.1128/mBio.00553-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Diamond, M. S. (2014). IFIT1: a dual sensor and effector molecule that detects non-2′-O methylated viral RNA and inhibits its translation. Cytokine Growth Factor Rev. 25, 543–550. doi: 10.1016/j.cytogfr.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Diamond, M. S., and Pierson, T. C. (2015). Molecular insight into dengue virus pathogenesis and its implications for disease control. Cell 162, 488–492. doi: 10.1016/j.cell.2015.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Flint, S. J., Enquist, L. W., Racaniello, V. R., and Skalka, A. M. (2009). Principles of Virology, 3rd Edn. Washington, DC: ASM Press.

Google Scholar

Gokhale, N. S., Vazquez, C., and Horner, S. M. (2014). Hepatitis C virus. Strategies to evade antiviral responses. Future Virol. 9, 1061–1075. doi: 10.2217/fvl.14.89

CrossRef Full Text | Google Scholar

Guedj, J., Dahari, H., Rong, L., Sansone, N. D., Nettles, R. E., Cotler, S. J., et al. (2013). Modeling shows that the NS5A inhibitor daclatasvir has two modes of action and yields a shorter estimate of the hepatitis C virus half-life. Proc. Natl. Acad. Sci. U.S.A. 110, 3991–3996. doi: 10.1073/pnas.1203110110

PubMed Abstract | CrossRef Full Text | Google Scholar

Haller, O., Staeheli, P., and Kochs, G. (2007). Interferon-induced Mx proteins in antiviral host defense. Biochimie 89, 812–818. doi: 10.1016/j.biochi.2007.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Heldt, F. S., Frensing, T., Pflugmacher, A., Gröpler, R., Peschel, B., and Reichl, U. (2013). Multiscale modeling of influenza A virus infection supports the development of direct-acting antivirals. PLoS Comput. Biol. 9:e1003372. doi: 10.1371/journal.pcbi.1003372

PubMed Abstract | CrossRef Full Text | Google Scholar

Heldt, F. S., Frensing, T., and Reichl, U. (2012). Modeling the intracellular dynamics of influenza virus replication to understand the control of viral RNA synthesis. J. Virol. 86, 7806–7817. doi: 10.1128/JVI.00080-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho, D. D., Neumann, A. U., Perelson, A. S., Chen, W., Leonard, J. M., and Markowitz, M. (1995). Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature 373, 123–126. doi: 10.1038/373123a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kato, F., Ishida, Y., Oishi, S., Fujii, N., Watanabe, S., Vasudevan, S. G., et al. (2016). Novel antiviral activity of bromocriptine against dengue virus replication. Antiviral Res. 131, 141–147. doi: 10.1016/j.antiviral.2016.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreutz, C., Raue, A., Kaschek, D., and Timmer, J. (2013). Profile likelihood in systems biology. FEBS J. 280, 2564–2571. doi: 10.1111/febs.12276

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhn, R. J., Zhang, W., Rossmann, M. G., Pletnev, S. V., Corver, J., Lenches, E., et al. (2002). Structure of dengue virus: implications for flavivirus organization, maturation, and fusion. Cell 108, 717–725. doi: 10.1016/S0092-8674(02)00660-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindenbach, BD, Evans, MJ, Syder, AJ, Wölk, B, Tellinghuisen, TL, Liu, CC, et al (2005). Complete Replication of Hepatitis C Virus in Cell Culture. Science. 309, 623–626. doi: 10.1126/science.1114016

PubMed Abstract | CrossRef Full Text | Google Scholar

Laske, T., Heldt, F. S., Hoffmann, H., Frensing, T., and Reichl, U. (2016). Modeling the intracellular replication of influenza A virus in the presence of defective interfering RNAs. Virus Res. 213, 90–99. doi: 10.1016/j.virusres.2015.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Maiwald, T., Hass, H., Steiert, B., Vanlier, J., Engesser, R., Raue, A., et al. (2016). Driving the model to its limit: profile likelihood based model reduction. PLoS ONE 11:e0162366. doi: 10.1371/journal.pone.0162366

PubMed Abstract | CrossRef Full Text | Google Scholar

Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. (2008). A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254, 178–196. doi: 10.1016/j.jtbi.2008.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Martensen, P. M., and Justesen, J. (2004). Small ISGs coming forward. J. Interferon Cytokine Res. 24, 1–19. doi: 10.1089/107999004772719864

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazzon, M., Jones, M., Davidson, A., Chain, B., and Jacobs, M. (2009). Dengue virus NS5 inhibits interferon-α signaling by blocking signal transducer and activator of transcription 2 phosphorylation. J. Infect. Dis. 200, 1261–1270. doi: 10.1086/605847

PubMed Abstract | CrossRef Full Text | Google Scholar

Metz, P., Dazert, E., Ruggieri, A., Mazur, J., Kaderali, L., Kaul, A., et al. (2012). Identification of type I and type II interferon-induced effectors controlling hepatitis C virus replication. Hepatology 56, 2082–2093. doi: 10.1002/hep.25908

PubMed Abstract | CrossRef Full Text | Google Scholar

Metz, P., Reuter, A., Bender, S., and Bartenschlager, R. (2013). Interferon-stimulated genes and their role in controlling hepatitis C virus. J. Hepatol. 59, 1331–1341. doi: 10.1016/j.jhep.2013.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrison, J., Aguirre, S., and Fernandez-Sesma, A. (2012). Innate immunity evasion by dengue virus. Viruses 4, 397–413. doi: 10.3390/v4030397

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Jordán, J. L., and Fredericksen, B. L. (2010). How flaviviruses activate and suppress the interferon response. Viruses 2, 676–691. doi: 10.3390/v2020676

PubMed Abstract | CrossRef Full Text | Google Scholar

Munoz-Jordan, J. L., Sanchez-Burgos, G. G., Laurent-Rolle, M., and Garcia-Sastre, A. (2003). Inhibition of interferon signaling by dengue virus. Proc. Natl. Acad. Sci. U.S.A. 100, 14333–14338. doi: 10.1073/pnas.2335168100

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, N. E. A., Quam, M. B., and Wilder-Smith, A. (2013). Epidemiology of dengue: past, present and future prospects. Clin. Epidemiol. 5, 299–309. doi: 10.2147/CLEP.S34440

PubMed Abstract | CrossRef Full Text | Google Scholar

Nasirudeen, A. M. A., Wong, H. H., Thien, P., Xu, S., Lam, K. P., and Liu, D. X. (2011). RIG-I, MDA5 and TLR3 synergistically play an important role in restriction of dengue virus infection. PLoS Negl. Trop. Dis. 5:e926. doi: 10.1371/journal.pntd.0000926

PubMed Abstract | CrossRef Full Text | Google Scholar

Neufeldt, C. J., Cortese, M., Acosta, E. G., and Bartenschlager, R. (2018). Rewiring cellular networks by members of the flaviviridae family. Nat. Rev. Microbiol. 16, 125–142. doi: 10.1038/nrmicro.2017.170

PubMed Abstract | CrossRef Full Text | Google Scholar

Neumann, A. U. (1998). Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-alpha therapy. Science 282, 103–107. doi: 10.1126/science.282.5386.103

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowak, M. A., and Bangham, C. R. (1996). Population dynamics of immune responses to persistent viruses. Science 272, 74–79. doi: 10.1126/science.272.5258.74

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowak, M. A., Bonhoeffer, S., Hill, A. M., Boehme, R., Thomas, H. C., and McDade, H. (1996). Viral dynamics in hepatitis B virus infection. Proc. Natl. Acad. Sci. U.S.A. 93, 4398–4402. doi: 10.1073/pnas.93.9.4398

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A. S. (2002). Modelling viral and immune system dynamics. Nat. Rev. Immunol. 2, 28–36. doi: 10.1038/nri700

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A. S., Essunger, P., Cao, Y., Vesanen, M., Hurley, A., Saksela, K., et al. (1997). Decay characteristics of HIV-1-infected compartments during combination therapy. Nature 387, 188–191. doi: 10.1038/387188a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A. S., and Guedj, J. (2015). Modelling hepatitis C therapy-predicting effects of treatment. Nat. Rev. Gastroenterol. Hepatol. 12, 437–445. doi: 10.1038/nrgastro.2015.97

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A. S., Neumann, A. U., Markowitz, M., Leonard, J. M., and Ho, D. D. (1996). HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science 271, 1582–1586. doi: 10.1126/science.271.5255.1582

PubMed Abstract | CrossRef Full Text | Google Scholar

Perelson, A. S., and Ribeiro, R. M. (2013). Modeling the within-host dynamics of HIV infection. BMC Biol. 11:96. doi: 10.1186/1741-7007-11-96

PubMed Abstract | CrossRef Full Text | Google Scholar

Perry, S. T., Prestwood, T. R., Lada, S. M., Benedict, C. A., and Shresta, S. (2009). Cardif-mediated signaling controls the initial innate response to dengue virus in vivo. J. Virol. 83, 8276–8281. doi: 10.1128/JVI.00365-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Rand, U., Rinas, M., Schwerk, J., Nöhren, G., Linnes, M., Kröger, A., et al. (2012). Multi-layered stochasticity and paracrine signal propagation shape the type-I interferon response. Mol. Syst. Biol. 8:584. doi: 10.1038/msb.2012.17

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranjit, S., and Kissoon, N. (2011). Dengue hemorrhagic fever and shock syndromes. Pediatr. Crit. Care Med. 12, 90–100. doi: 10.1097/PCC.0b013e3181e911a7

PubMed Abstract | CrossRef Full Text | Google Scholar

Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmüller, U., et al. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics 25, 1923–1929. doi: 10.1093/bioinformatics/btp358

PubMed Abstract | CrossRef Full Text | Google Scholar

Raue, A., Schilling, M., Bachmann, J., Matteson, A., Schelker, M., Schelke, M., et al. (2013). Lessons learned from quantitative dynamical modeling in systems biology. PLoS ONE 8:e74335. doi: 10.1371/annotation/ea0193d8-1f7f-492a-b0b7-d877629fdcee

PubMed Abstract | CrossRef Full Text | Google Scholar

Raue, A., Steiert, B., Schelker, M., Kreutz, C., Maiwald, T., Hass, H., et al. (2015). Data2Dynamics: a modeling environment tailored to parameter estimation in dynamical systems. Bioinformatics 31, 3558–3560. doi: 10.1093/bioinformatics/btv405

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodenhuis-Zybert, I. A., Wilschut, J., and Smit, J. M. (2010). Dengue virus life cycle: viral and host factors modulating infectivity. Cell. Mol. Life Sci. 67, 2773–2786. doi: 10.1007/s00018-010-0357-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Fernandez, M., and Banga, J. R. (2010). SensSB: a software toolbox for the development and sensitivity analysis of systems biology models. Bioinformatics 26, 1675–1676. doi: 10.1093/bioinformatics/btq242

PubMed Abstract | CrossRef Full Text | Google Scholar

Romero-Brey, I., Merz, A., Chiramel, A., Lee, J. Y., Chlanda, P., Haselman, U., et al. (2012). Three-dimensional architecture and biogenesis of membrane structures associated with hepatitis C virus replication. PLoS Pathog. 8:e1003056. doi: 10.1371/journal.ppat.1003056

PubMed Abstract | CrossRef Full Text | Google Scholar

Rong, L., and Perelson, A. S. (2009). Modeling HIV persistence, the latent reservoir, and viral blips. J. Theor. Biol. 260, 308–331. doi: 10.1016/j.jtbi.2009.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronni, T., Melén, K., Malygin, A., and Julkunen, I. (1993). Control of IFN-inducible MxA gene expression in human cells. J. Immunol. 150, 1715–1726.

PubMed Abstract | Google Scholar

Safran, M., Dalah, I., Alexander, J., Rosen, N., Iny Stein, T., Shmoish, M., et al. (2010). GeneCards version 3: the human gene integrator. Database 2010:baq020. doi: 10.1093/database/baq020

PubMed Abstract | CrossRef Full Text | Google Scholar

Scaturro, P., Cortese, M., Chatel-Chaix, L., Fischl, W., and Bartenschlager, R. (2015). Dengue virus non-structural protein 1 modulates infectious particle production via interaction with the structural proteins. PLoS Pathog. 11:e1005277. doi: 10.1371/journal.ppat.1005277

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmid, B., Rinas, M., Ruggieri, A., Acosta, E. G., Bartenschlager, M., Reuter, A., et al. (2015). Live cell analysis and mathematical modeling identify determinants of attenuation of dengue virus 2'-O-methylation mutant. PLoS Pathog. 11:e1005345. doi: 10.1371/journal.ppat.1005345

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoggins, J. W., and Rice, C. M. (2011). Interferon-stimulated genes and their antiviral effector functions. Curr. Opin. Virol. 1, 519–525. doi: 10.1016/j.coviro.2011.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Screaton, G., Mongkolsapaya, J., Yacoub, S., and Roberts, C. (2015). New insights into the immunopathology and control of dengue virus infection. Nat. Rev. Immunol. 15, 745–759. doi: 10.1038/nri3916

PubMed Abstract | CrossRef Full Text | Google Scholar

Scutigliani, E. M., and Kikkert, M. (2017). Interaction of the innate immune system with positive-strand RNA virus replication organelles. Cytokine Growth Factor Rev. 37, 17–27. doi: 10.1016/j.cytogfr.2017.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Shresta, S., Kyle, J. L., Snider, H. M., Basavapatna, M., Beatty, P. R., and Harris, E. (2004). Interferon-dependent immunity is essential for resistance to primary dengue virus infection in mice, whereas T- and B-cell-dependent immunity are less critical. J. Virol. 78, 2701–2710. doi: 10.1128/JVI.78.6.2701-2710.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Stafford, M. A., Corey, L., Cao, Y., Daar, E. S., Ho, D. D., and Perelson, A. S. (2000). Modeling plasma virus concentration during primary HIV infection. J. Theor. Biol. 203, 285–301. doi: 10.1006/jtbi.2000.1076

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuiskunen Bäck, A., and Lundkvist, Å. (2013). Dengue viruses – an overview. Infect. Ecol. Epidemiol. 3:19839. doi: 10.3402/iee.v3i0.19839

CrossRef Full Text | Google Scholar

Wei, X., Ghosh, S. K., Taylor, M. E., Johnson, V. A., Emini, E. A., Deutsch, P., et al. (1995). Viral dynamics in human immunodeficiency virus type 1 infection. Nature 373, 117–122. doi: 10.1038/373117a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Welsch, S., Miller, S., Romero-Brey, I., Merz, A., Bleck, C. K. E., Walther, P., et al. (2009). Composition and three-dimensional architecture of the dengue virus replication and assembly sites. Cell Host Microbe 5, 365–375. doi: 10.1016/j.chom.2009.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

WHO (2012). Global Strategy for Dengue Prevention and Control 2012–2020. World Health Organization.

Google Scholar

World Health Organization (2016). Dengue Vaccine: WHO Position Paper – July 2016. Weekly Epidemiological Record (WER). World Health Organization, 2016, 349–364.

Google Scholar

Wulff, N. H., Tzatzaris, M., and Young, P. J. (2012). Monte carlo simulation of the spearman-kaerber TCID50. J. Clin. Bioinformatics 2:5. doi: 10.1186/2043-9113-2-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Xagorari, A., and Chlichlia, K. (2008). Toll-like receptors and viruses: induction of innate antiviral immune responses. Open Microbiol. J. 2, 49–59. doi: 10.2174/1874285800802010049

PubMed Abstract | CrossRef Full Text | Google Scholar

Zitzmann, C., and Kaderali, L. (2018). Mathematical analysis of viral replication dynamics and antiviral treatment strategies: from basic models to age-based multi-scale modeling. Front. Microbiol. 9:1546. doi: 10.3389/fmicb.2018.01546

PubMed Abstract | CrossRef Full Text | Google Scholar