The standard model that has been used and modified for studying hepatitis C viral dynamics is the Neumann et al. model [20]. The three variables this model keeps track of are the target cells T, in Equation (1a), the infected cells I in Equation (1b) and the extracellular virus V in Equation (1c). The target cells T are produced at constant rate s, and decreased by the number of cells infected by virus in blood V at constant rate β and their death rate d. The infected cells, I, increase with the new infections at rate βV(t)T(t) and die at constant rate δ. The virus V is produced at rate p by each infected cell and is cleared at constant rate c. The ϵ term denotes the effectiveness of the anti-viral treatment that decreases viral production from p to (1 − ϵ)p. The previously established ensemble of ODEs for this model is:

From the mathematical perspective, the model is simple and can be solved analytically by appropriate assumptions.

The multiscale model of HCV dynamics has been formulated in recent years [25, 41, 43] in order to study HCV dynamics in patients and decide among various treatment options. Figure 1 depicts a systematic overview of the model, which is described by the following partial differential equations (PDEs):

subject to the boundary conditions I(0, t) = βV(t)T(t), I(a, 0) = Ī(a), R(0, t) = 1, and R(a,0)=R-(a).

The variable I(a, t) for infected cells, that existed in the standard model as simply I(t), and the newly introduced variable R(a, t) for intracellular viral RNA (vRNA) both depend on the age of infection a and the time duration from therapy initiation t. Hence, a and t are two different times, with the use of partial derivatives in Equations (2b) and (2d). Model parameters of T, Equation (2a), and I, Equation (2b), are similar to the standard model. The quantity of vRNA R, in Equation (2d), depends on its synthesis α and its degradation μ and secretion from the cell as virus particles ρ. The quantity of extracellular virus V shown in Equation (2c) depends on the number of assembled and released virions and their clearance rate c. A schematic description of the multiscale model is shown in Figure 1.

An important consideration in this model is that the treatment starts after the infection has reached its steady state. The steady states of the different functions are R-(a,t), Ī(a, t), V- and T-. Given N, the total number of virions produced by a cell in its life-span, it was shown in Rong et al. [41] that those values are:

Unlike the standard model, three different antiviral effects of therapy can be simulated in the multiscale model. The decrease in viral RNA synthesis is represented by ε_α, the reduction in secretion by ε_s and the increase in viral degradation by κ ≥ 1.

Through the method of characteristics, as was derived in Rong et al. [41], an analytical solution was found for the variable R(a, t). The same method was applied to derive a solution for I(a, t). The ensemble of Equations (2) represents the full model. The analytical solutions for R(a, t) and I(a, t) were described as follows:

From the system of Equations (2) it can be noticed that computing V(t) necessitates an integral. If a < t, in other words the cell age is younger than the time of treatment, i.e., infection occurs after initiation of treatment, the term I(a, t) of the integral in Equation (2) depends itself on V(t) and T(t) by consideration of Equation (5). As was shown in Guedj et al., Rong et al., Rong and Perelson [25, 41, 43], this makes the analytical solution for V(t) approximative.

In this section we wish to approximate the function y(t) for solving the differential equation dydt=f(t,y), where f(t, y) is known. Finite-difference methods of higher order are used. Given a known value of y at time step n, y_n, compute the value at the next time step y_n+1 by means of f(t_n, y_n).

Since serum HCV RNA is stable in chronic HCV subjects, we assume that the model is at steady state before the start of treatment, as shown in Equation (3) and discussed in Rong et al. [41]. While individuals with recent HCV infection might not be at steady state at initiation of DAA therapy [3], it may be feasible to assume a quasi-steady state with lower pre-treatment HCV RNA levels, as observed in some subjects, which could be adjusted by changing model parameters.

We implemented the multiscale model with the Rosenbrock method of order S = 4. We define:

where yn=(TnVn) and

We note that the term (1-ϵs)∫0∞ρ(a)R(a,t)I(a,t)da is constant at each iteration and is computed between each one once. We can now explicitly write the values of the constants and describe the whole scheme. The Equations (13) and (14) become:

Note that the computations of g₃ in Equation (19) and g₄ in Equation (20) use the same value of f, which implies that it needs to be computed only once. The error term also disregards the value of g₃.

The scheme for one iteration is shown in Algorithm 1. An important observation is how the error term e is computed. It derives only from the error induced by the Rosenbrock iteration, not the computation of the integral term (1-ϵs)∫0∞ρ(a)R(a,t)I(a,t)da.

The system (I/hγ-Jn)gi=f(yn+∑j=1i-1αijgj)+1h∑j=1i-1γijgj needs to be solved for four different values of i. Since the left hand matrix is constant, it is decomposed into its LU decomposition once. The Jama package [53] is used to perform the LU decomposition for solving the system in a highly efficient manner.

The short-term approximation [41] is shown in Figure 2 in blue. The results are shown for only 2 days since the behavior is smooth and conserved on longer time scales in agreement with previous results (Figure 2A in [41]). It is clear that after 12 h the value converges far from the PDE solution. This is expected since the effect of the treatment on the infection rate is not taken into account. In practice, most simulations are valuable for more than half a day and are focused on the long-term approximation.

In Rong and Perelson [43] it is hypothesized that the long-term approximation is an underestimate of the PDE model solution since some infection events are being ignored. However, with realistic parameters characteristic of potent therapy, the difference between them is very small. Figure 3 depicts the result of our scheme, measuring the difference between the long-term approximation and the solution of the multiscale PDE model (Figure 2B in [41]). We show that the long-term approximation is converging and we are indeed obtaining that the long-term approximation is an underestimate of the PDE model solution (as the difference between the long term approximation and the solution of the PDE model is negative). The difference of less than 0.01 that we are obtaining is very small compared to the y-axis units shown in Figure 1 of Rong and Perelson [43]. Therefore, it is possible to view in Figure 3 this small effect that is an improvement over what was inferred in Rong and Perelson [43] regarding the long-term approximation. This finding is reiterated here for completeness and is referred to below in the sub-section about the Rosenbrock method.

An important distinction of the age-based PDE HCV model equations is their stiffness, or how they can be numerically unstable even with a small stepsize. We have verified that with small time-steps the results are similar for all methods. But as we increase the size of the steps, the instability of the equations can be observed over 20 days. All explicit methods such as the standard version of Runge-Kutta quickly diverge from the solution with small time-steps. On the other hand an implicit Runge-Kutta method tends to the correct solution and the Rosenbrock method does so efficiently as reported below.

As previously mentioned, the integral term (1-ϵs)∫0∞ρ(a)R(a,t)I(a,t)da impedes the use of adaptive stepsize when a canned solver is used. Such methods are essential to allow fast computations of complex models like ours. We present in Table 5 the number of iterations computed using the Rosenbrock method, shown in the previous section, in comparison with the default canned method for three different time step values of h and h_a. With the smallest time step, the Rosenbrock method is more than five times as efficient as the default canned method, giving results of similar accuracy. This is crucial since the methods with fixed time step can take up to tens of seconds per day of simulation with h = 0.001 and h_a = 0.01. For most research needs, simulations are expected to be performed numerous times for various instances and parameter values. Under those conditions a five-fold increase in speed, through the decrease in the number of iterations, provides a much needed advantage when testing large number of parameters over long time periods. Interestingly with the highest value of h, the Rosenbrock method takes additional iterations. This is due to the greediness of the time step adjustment, which when increased too quickly and thereby inducing a large error backtracks and starts again with a smaller value.

Simpler methods than Rosenbrock exist with adaptive stepsize, such as Dormand–Prince [48], Fehlberg (RKF) [55] and Cash-Karp (RKCK) [56]. These are explicit methods that can be easily implemented but are vulnerable to the stiffness of the equations, while Rosenbrock is an implicit method that is indeed found to be more stable. In the following comparisons we demonstrate that the Rosenbrock method is more efficient than the Runge-Kutta 4th order method and at the same time it is more stable than both the Runge-Kutta 4th order and the Dormand-Prince methods.

First, a comparison of the computing time between the Rosenbrock method and the other methods mentioned above is reported in Table 6 with starting time step values of h = 0.001 and h_a = 0.01 for 14 days of simulation. Each run with a certain method was repeated 10 times, after which the mean and standard deviation was taken. The Rosenbrock method is over 8 times faster than the standard Runge-Kutta 4th order method, which is a non-adaptive time step method. While Dormand–Prince is an adaptive time step method that is even faster than Rosenbrock, as will be shown below, it is an explicit method that does not always converge to the correct solution. The fastest time of execution for the Dormand-Prince method is not surprising since an explicit method is easier to evaluate computationally than an implicit method, in particular because there is no need for matrix inversions.

Second, a comparison of the methods behavior over 14 days of simulation is presented in Figure 4. All the methods were started with time steps of h = 0.1 and h_a = 0.1. The baseline (the default implementation in the SciPy library of an ODE solver) was computed for h = 0.001 and h_a = 0.01 since those time step values are small enough to ensure that all methods converge to the same solution. The advantage of the Rosenbrock method in terms of stability and convergence to the correct solution is clearly noticed.

Finally, we also report in Table 7 the sensitivity of the Rosenbrock method. We compute the changed value of the result after 2 days of simulation given a variation of ±10% on 10 parameters. We notice that most of the perturbations are symmetrical. We omit ε_α and ε_s since they effectively determine the fraction of α and ρ used in the model. The results (Table 7) indicate that the Rosenbrock method is sufficiently robust for our needs.

We have developed a user-friendly simulator with a GUI for the multiscale model (Figure 5) that is freely available at http://www.cs.bgu.ac.il/~dbarash/HCVsimulator. The parameters file input is shown in Figure 6. For illustration, two example simulations are provided, one over 2 days (Figure 7) and one over 28 days (Figure 8). In Figure 7 we show how different parameter ensembles can be displayed simultaneously. In that case we chose ε_α and ε_s to be 0 or 0.99, resulting in three curves. We can observe that increasing the ε_s parameter, which decreases the assembly/secretion of intracellular virions, has a strong short term effect that tempers quickly. In contrast the reduction in the synthesis of vRNA, modeled by an increase of the parameter ε_α, shows a slower but more efficient effect that exhibits stronger results after half a day. Combining both factors shows a multiphasic viral decline which was observed under treatment with HCV NS5A inhibitors [25]. In Figure 8 we present the difference between the long-term approximation and the Rosenbrock method for those three sets of parameters over longer (than 2 days) treatment durations (e.g., 28 days). Although the highest error is present when ε_s = 0, the error remains below a tenth of the differences of the log values. Thus, in all cases the trends from the first 2 days continue in the next 4 weeks and the error grows minimally.

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.