We developed a mathematical model to represent the intracellular replication of HCV shown schematically in Figure 1. The model allows the study of aspects such as translation of positive-strand HCV RNA after cell entry, transfer of the positive strand to the membranous web where it is used for replication, production of negative- and positive-strand HCV RNA within replication complexes and secretion of positive-strand RNA as virions. The replication of HCV RNA has been studied in detail c.f. (Chatel-Chaix and Bartenschlager, 2014; Li et al., 2015).

The system of ordinary differential equations used to represent the dynamics of intracellular infection over time is

where R_t represents positive strand HCV RNA used for translation, R_c represents positive strands within replication complexes used for replication, R_m represents minus (or negative) strand RNA and a represents the time a cell has been infected. Positive strand HCV RNA forms the viral genome. After cell entry, cellular machinery translates this positive strand RNA into a polyprotein in the cytoplasm (Shi and Lai, 2006). However, after polyproteins are made the positive strand must also be used for replication and must be copied into minus stand RNA. We assume that the positive-strand HCV RNA used for translation (R_t in Equation 1) moves from the cytoplasm into what is called the membranous web and interacts with the proteins needed for replication to become a species we call R_c at rate σ per strand. We also assume the positive strand in the cytoplasm, R_t has a natural decay rate of μ_t per strand. Lastly, positive strands need to be assembled into virions which are then exported from the infected cell. Virion assembly occurs in association with cytosolic lipid droplets (Chatel-Chaix and Bartenschlager, 2014). As it is not clear whether the positive strand RNA in the membranous web needs first to be transported into the cytosol for viral assembly, we assume both R_t and R_c can be assembled into virions and exported at rate ρ(a). The time-dependence of ρ will be discussed below. Further, we assume positive-strand HCV RNA in the replication complex (R_c) can move out of the replication complex and membranous web and back into the cytoplasm to become R_t at rate θ. More detailed models can be developed that separate virion assembly from secretion and that include a separate compartment of positive strand RNA used for virion assembly (cf. Benzine et al., 2017), but here for simplicity we have combined these steps.

Minus-strand HCV RNA (R_m) is formed in the replication complex by copying the positive strand R_c at maximum rate r. As in Guedj and Neumann (2010), it is assumed that host factors limit the replication of negative-strand RNAs, so that as the maximum number R_max is reached replication slows according to a logistic growth law. The positive strands in replication complexes, R_c, are copied from the negative strand template at rate α per template. We consider that both R_c and R_m are in the replication complex and decay at the same per capita rate μ_c.

In order to have a positive equilibrium when the model represents an established infection, the parameters need to satisfy the relations: ϕ2>σθϕ1 and αr>(ϕ2-σθϕ1)μc in which ϕ₁ = θ + ρ + μ_t and ϕ₂ = σ + ρ + μ_c.

We also analyze in vivo data in which the effects of antiviral treatment on kinetics of HCV RNA levels in plasma are measured. To fit this data we introduced a new multiscale model depicted in Figure 2.

The intracellular portion of the multiscale model with treatment is represented by the following partial differential equations (PDEs) in which t represents clock time and a the age of an infected cell:

We have assumed that intracellular infection is initiated by the introduction of R_t0 positive HCV RNA strands into the cytoplasm of a cell. Typically, we shall assume that infection is the result of a single virion, carrying a single positive-strand HCV RNA, entering a cell, so that R_t0 = 1. Further, we shall assume that the individual's being treated with antivirals are chronically infected and have reached steady state in which R¯t(a), R¯c(a) and R¯m(a) are the steady state distributions of positive-strand HCV RNA in translation and in replication complexes and negative-strand HCV RNA in replication complexes, respectively, in the absence of treatment and are given by the steady state solutions of the ODEs in Equation (1). Further, we let ϵ_α be the effectiveness of therapy in decreasing or blocking positive-strand RNA replication, ϵ_r the effectiveness of therapy in decreasing or blocking negative-strand RNA replication, and ϵ_s the effectiveness of therapy in decreasing or blocking secretion of positive-strand RNA, where for each of the ϵ's, ϵ = 1 corresponds to a 100% effective drug and ϵ = 0 corresponds to a completely ineffective or absent drug. Further κ_t is a factor by which therapy changes the degradation rate of positive-strand RNA used for translation and κ_c is the factor by which therapy changes the degradation rate of both positive and negative strand RNA in replication complexes.

To complete the multiscale model, we then coupled the intracellular model to an established HCV cellular infection model (Equation 5) (Neumann et al., 1998; Canini and Perelson, 2014).

in which, T are target cells, I, infected cells and V the HCV RNA concentration in plasma. Target cells become infected at rate β, have a constant source rate s and a natural per capita decay rate d. The parameter δ(a) represents the death rate of an infected cell of age “a” and the effects of therapy on the virus export are given by ϵ_s. Here for simplicity we shall only analyze the case in which δ(a) is a constant, δ. Virus in the plasma is cleared from the circulation at per capita rate c. Here we have assumed that at t = 0, the time therapy starts, the system is in steady state, where Ī(a) is the steady-state distribution of infected cells, which can be shown to be Ī(a) = βV0T0e-δa. T0=cβN, where N is the steady state total amount of virus secreted by an infected cell over its lifetime, N=ρ∫0∞(R¯t(a)+R¯c(a)e−δada, and V0=s-dT0βT0. See Rong et al. (2013). The coupling between the intracellular and extracellular models occurs through the equation for V, the virus in the plasma. The amount of plasma virus depends on the number of infected cells and number of virions being packaged and exported per infected cell. This coupling has been used before (Guedj et al., 2013; Rong et al., 2013; Rong and Perelson, 2013).

The model equations were discretized in space, i.e., age, and integrated in time using the method of lines (MOL) approach (Sadiku and Obiozor, 2000; Shakeri and Dehghan, 2008) where the partial derivatives in age were approximated by finite-differences and the solution at the grid points was integrated along lines in time. We integrated the equations using the Matlab® ordinary differential equation Runge-Kutta solver ode45.

The domain was discretized on a uniform grid of 201 mesh points between ages 0 and 50 days, as it is unlikely for an infected cell to live longer than this. The boundary of domain at a = 50 was defined as a simple outflow boundary condition and was incorporated into the numerical solution by linearly extrapolating the solution to two buffer grid points outside the domain. The partial derivatives in age were approximated with fourth-order centered finite differences.

We verified the convergence of the numerical solution to an accuracy of 10⁻³ by varying the number of spatial grid points and the time integration error tolerance.

The simulation time was varied according to the length of time that virus was detected in plasma after therapy initiation in the data we analyzed. The computer run time were typically a few seconds for a single simulation using a laptop computer.

We used the Matlab® nonlinear optimization program fmincon to fit the solutions of the model to the experimental data by minimizing the L₂ norm of the residual difference between the model solution and the data. This routine was chosen due to the possibility of specifying lower and upper bounds for the parameters we wanted to estimate. The algorithm we used was “interior-point” as it satisfies the bounds at all iterations.

The data we fit to validate the model was obtained from different sources. We extracted in vitro data from Keum et al. (2012) and Binder et al. (2013) using the on-line tool WebPlotDigitizer (Rohatgi, 2016). We also fit clinical trial data from Guedj et al. (2013) that we had access to.

Because our models have a large number of parameters we numerically approximated the Hessian of the objective function at the optimal parameter values. At a minimum, the gradient of the objective function is zero. If an eigenvalue of the Hessian is zero at the minimum, then the gradient remains zero along the direction of the associated eigenvector. That is, the solution is not unique (identifiable) (Beck, 2014). Here, for each of the data sets we fit, at the optimum, all of the eigenvalues of the Hessian were positive, and the condition number was below 10⁴, indicating that the parameters were locally identifiable.

To validate the intracellular mathematical model, we first compared the results of Equation (1) to transfection experiments performed by Binder et al. (2013). In that paper the authors used two distinct cell lines to assess HCV RNA replication over 72 h: (a) Huh7-Lunet cells which are highly permissive to HCV RNA replication and (b) Huh7 cells (Huh7-lp) which presents lower levels of HCV RNA replication. They measured positive-strand and negative-strand RNA by strand specific quantitative Northern blotting. Binder et al. (2013) developed a complex mathematical model that included 13 molecular species with 16 parameters in two compartments: the cytoplasm and a replication compartment.

Using the three equation mathematical model, Equation (1), we were able to fit the dynamics of both positive and negative strand HCV RNA in both the high and low permissive cell lines (Figure 3). Our model was able to replicate the initial decay seen after transfection with both types of cells and the plateau during the 72 h measured (Figures 3A,B).

In fitting the data, the parameters used to describe the age-dependent virion export rate, ρ(a), were fixed at ρ = 0.1 d⁻¹, τ = 0.5 d⁻¹ and k = 0.8 d⁻¹. We set τ at 0.5 d⁻¹ based on the fact that Keum et al. (2012) could not detect any extracellular virus until 12 h post-infection. We further tested different values of τ and k and chose the values that gave the best fits to the data in both the Binder and Keum experiments. Choosing the export rate as a time-dependent function rather than a constant allowed us to have an initial delay followed by a smooth transition to the maximum export ρ. Regarding the maximum export rate, ρ, we at first chose the value estimated by Guedj et al. (2013) based on fitting in vivo data. However, using this value did not give good fits to the in vitro data. We then scanned through different values and chose the one giving the best fit. HCV uses the host cell export machinery and thus it is not surprising that these parameters differ between in vitro and in vivo systems.

The initial number of HCV positive strands introduced into these cells to initiate HCV replication in this in vitro system was R_t0 = 4, 000 molecules cell⁻¹. Other parameters of the model were estimated using the fmincon routine in Matlab and are shown in Table 1.

Another form of validation we performed was testing the model predictions by comparing to positive-strand measurements using a replication deficient replicon (Binder et al., 2013). By setting the rate at which positive-strand RNA goes from use in translation to use in replication (σ) to zero we could compare the results obtained with the model to the measurements reported by Binder et al. (2013) Without replication, the initial amount of transfected HCV RNA decays exponentially and no negative-strand is formed. Further as Binder et al. show the decay of positive strand RNA is similar in both the high permissive and low permissive cell lines. We simulated the intracellular model with the parameters that were estimated for the highly-permissive cell line (Table 1) and the results are shown in Figure 4. The results using the parameters for the low permissive cell line are the same.

We also compared the intracellular model to experiments in vitro performed by Keum et al. (2012) in which a high multiplicity of infection was used (MOI = 5 or 6) so that only one round of infection occurred. Theoretically, with an MOI of 5, 99.3% of cells should be infected with a least one infectious virion (Keum et al., 2012). A cell culture adapted HCV, JFH-m4, was incubated with Huh7.5.1 cells for 3 h to initiate infection. At subsequent times cells and supernatant were harvested to measure HCV RNA levels intra-cellularly and the amount secreted into the medium. Keum et al. quantified the number of positive and negative HCV RNA strands using real-time RT-PCR. As shown in Figure 6 the number of cell-associated positive strands initially decreased reaching a minimum of about 1 positive strand per cell at 6 h post-infection (pi). Intracellular negative strand, which serves as a template for making new positive strands, was first detected at 6 h pi. Our model was able to reproduce the observed intracellular HCV RNA dynamics (Figure 6) as well as the dynamics of positive strand HCV RNA secreted into the media (Figure 7). As before we fixed the export rate with ρ = 0.1 d⁻¹, τ = 0.5 d⁻¹, and k = 0.8 d⁻¹. The initial time t₀ = 0, R_t0 = 12.8 and no therapy was given (Figure 7). Other parameters were estimated and were found to be α = 30 d⁻¹, μ_t = 24 d⁻¹, r = 3.18 d⁻¹, μ_c = 1.05 d⁻¹, R_max = 100 molecules, σ = 0.1 d⁻¹ and θ = 1.2 d⁻¹. As both the cell line and virus used in these experiments are different than the ones used by Binder et al. (2013), it is surprising that resulting parameters do not differ very much from those we estimated in the previous section for high and low-permissive cells.

We validated the coupled multiscale model by fitting Equations (4) and (5) to data obtained from patients treated with one dose of 10 or 100 mg of daclatasvir (DCV) (Guedj et al., 2013). DCV inhibits the action of the HCV NS5A protein, which has been shown to play an important role in HCV RNA replication and secretion (Lee, 2013; Scheel and Rice, 2013). This data was previously analyzed by Guedj et al. (2013) using a much simpler multiscale model that only considered HCV positive strand RNA dynamics.

We assumed that the parameters that represent in vivo infection dynamics are different from those we estimated for in vitro infection as both the virus and target cells are different. We also assumed that there was no superinfection, so that only one virus infects each cell. Using the same approach as for the intracellular model, we performed a sensitivity analysis of the coupled model parameters in order to determine how sensitive the predicted viral load is to each parameter. We chose to vary each parameter one at a time and compared how they affected the predicted viral load at day 2 on therapy.

The sensitivity index was calculated using Equation (6) and the results are shown in Figure 8. Intracellular parameters such as the replication and decay rates of HCV RNA, α, r, μ_c are the ones which the viral load is most sensitive to. The parameters that represents the export rate, ρ, and infected cell decay rate δ, are also important to define the viral load.

A baseline in vivo set of parameters was fixed based on the literature: α = 30 d⁻¹, ρ = 8.18 d⁻¹, δ = 0.14 d⁻¹, and c = 22.3 d⁻¹ were taken from Rong et al. (2013) and ϵ_s = 0.998 was taken from Guedj et al. (2013). The remaining parameters were estimated and their values are shown in Table 2.

Figures 9, 10 depict the results obtained with the multiscale model for each patient. We fixed the replication rate of positive strand HCV RNA α = 30 d⁻¹ and considered no enhancement in HCV RNA decay with therapy, κ_t = κ_c = 1. Our model predicted that initiation of therapy affects the replication of both positive and negative strands and that initially there is a slightly increase in the number of positive strand HCV RNAs used for translation (Figure 10). This increase is most likely due to the fact that DCV effectively blocks secretion of positive strands thus allowing them to accumulate in the cytoplasm. Therapy also blocks the appearance of new replication complexes, which only decrease in the presence of the drug (Figure 10).

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.