Live Cell Analysis and Mathematical Modeling Identify Determinants of Attenuation of Dengue Virus 2'-O-Methylation Mutant.

Dengue virus (DENV) is the most common mosquito-transmitted virus infecting ~390 million people worldwide. In spite of this high medical relevance, neither a vaccine nor antiviral therapy is currently available. DENV elicits a strong interferon (IFN) response in infected cells, but at the same time actively counteracts IFN production and signaling. Although the kinetics of activation of this innate antiviral defense and the timing of viral counteraction critically determine the magnitude of infection and thus disease, quantitative and kinetic analyses are lacking and it remains poorly understood how DENV spreads in IFN-competent cell systems. To dissect the dynamics of replication versus antiviral defense at the single cell level, we generated a fully viable reporter DENV and host cells with authentic reporters for IFN-stimulated antiviral genes. We find that IFN controls DENV infection in a kinetically determined manner that at the single cell level is highly heterogeneous and stochastic. Even at high-dose, IFN does not fully protect all cells in the culture and, therefore, viral spread occurs even in the face of antiviral protection of naïve cells by IFN. By contrast, a vaccine candidate DENV mutant, which lacks 2'-O-methylation of viral RNA is profoundly attenuated in IFN-competent cells. Through mathematical modeling of time-resolved data and validation experiments we show that the primary determinant for attenuation is the accelerated kinetics of IFN production. This rapid induction triggered by mutant DENV precedes establishment of IFN-resistance in infected cells, thus causing a massive reduction of virus production rate. In contrast, accelerated protection of naïve cells by paracrine IFN action has negligible impact. In conclusion, these results show that attenuation of the 2'-O-methylation DENV mutant is primarily determined by kinetics of autocrine IFN action on infected cells.

The time-resolved data after wildtype DENV infection (cf. Fig 6B, C) suggest that the timing of viral spread as well as IFN release play an important role in the competition between virus and the IFN system. Regarding the high diffusion coefficients of DENV (2.6×10 4 µm 2 /h, [S1]) and IFN (1.1×10 5 µm 2 /h, [S2]), spread occurs over many cell diameters on the relevant time scale of hours and spatial gradients of DENV as well as IFN are negligible. Therefore, we developed a population-based delay-differential equation (DDE) model with uniform distribution of secreted IFN and extracellular virus.
The observed single-step growth curves of virus replication, virus production and IFN secretion after high dose infection with DENV (cf. Fig 4C, D) illustrate that these processes are heterogeneously distributed among single cells. As an approximation of the temporal distribution of these processes, we take into account the delay of virus replication, virus production and IFN secretion in form of constant time delays.
In our population-based model (Fig 6A) we initially consider a number of naïve cells S 0 which are susceptible to virus infection, as well as an initial extracellular viral load V 0 . Susceptible cells S can become infected by extracellular virus V with infection rate r V (cf. blue arrow in Fig 6A). Depending on the time elapsed since viruses have entered the host cells, infected cells I can acquire up to three different features in parallel: -After a time delay τ R , infected but not yet virus replicating cells I R turn into virus replicating cells I R . The time delay of virus replication τ R is related to the mean expression onset of virus replication in individual cells monitored after infection with DENV-faR (cf. Fig 5 and Fig 6B).
-After a time delay τ V , infected but not yet virus releasing cells I V become virus producing cells I V and release new generated infectious virus particles with the virus production rate v V . The time delay of virus production τ V represents the average time required to replicate, translate viral proteins, assemble and release new virus particles (cf. Fig 4C).
-After a time delay τ F , infected but not yet IFN expressing cells I F turn into IFN secreting cells I F and produce IFN F with the IFN secretion rate v F . The time delay of IFN secretion τ F corresponds to the mean onset of IFN expression. It includes the necessary time of viral recognition, downstream signal transduction and IFN transcription, translation and secretion (cf. Fig 4D).
Since the time delay τ V also comprises the duration of virus replication, the condition τ V > τ R must be fulfilled, while we do not require any temporal link between τ F and both other time delays. Whether the production of new virus particles takes longer than IFN secretion or vice versa will be determined by the parameter estimation.
According to our experimental study, IFN has two antiviral effects (cf. Fig 1A, B and S1 Fig): -Pre-stimulation with IFN demonstrated that secreted IFN can protect naïve cells against DENV infection. To integrate this paracrine IFN response in our model, we take into account an IFN-induced switch from susceptible cells to protected cells P with the protection rate r F (solid green arrow in Fig 6A).
-Stimulation of infected cells with IFN in an early phase after viral entry can activate antiviral defense mechanisms, which reduce virus replication and consequently virus production. To keep the model simple, so that its parameters can be identified from experimental data, we neglect a direct autocrine impact of IFN in the first instance (dashed green arrow in Fig 6A). After the model parameterization, we will examine the effect of secreted IFN on infected cells directly with an extended version of our model, which is described in detail in section 4.
In order to include the removal of extracellular virus and IFN by cellular uptake and, to a lesser extent, by extracellular degradation, we assume a decrease of virus and IFN with the rate constant of virus degradation d V and the rate constant of IFN degradation d F , respectively.
Moreover, we observed that infected cells proliferate less and die faster than noninfected cells. To compensate this propagation difference and, additionally, to keep the number of model parameters low, we consider an inhibited proliferation in infected cells, whereas susceptible and protected cells proliferate with the rates p S and p P .
To establish the DDE system of our model, particular attention must be paid to the formulation of the three features of the infected cells comprising virus replication, virus production and IFN expression, which begin at time τ R , τ V and τ F , respectively. As the time delays τ m for m ∈ {R, V, F} are free parameters and only bound to the condition τ V > τ R , the properties of the infected cells can arise at different times, overlap each other and take place in parallel. Therefore we consider each feature of the infected cells separately and calculate the number of inactive infected cells I m (t) and active infected cells I m (t) for m ∈ {R, V, F} depending on time t ∈ R. For this purpose we first require the total number of infected cells I at time t I(t) = I m (t) + I m (t) ∀m ∈ {R, V, F} (1.1) and the corresponding differential equatioṅ (1. 2) The number of inactive infected cells I m at time t consists of all cells which become infected in the time period [t − τ m , t], since after τ m inactive cells I m progress to active cells I m . Thus, I m is given by To obtain the differential equationİ m , we differentiate both sides of (1.3) with respect to time t and apply the fundamental theorem of calculus: Using the preceding equations, the differential equationİ m results froṁ Taking all considerations together, our model is described by the following DDE system: (1.6) The corresponding initial value problem of our model is defined by the DDE system (1.6) for time t ≥ 0 together with initial conditions, which we specify next. At the starting time t 0 = 0 only two initial values are unequal to zero, namely the number of susceptible cells S(0) = S 0 as well as the extracellular viral load V (0) = V 0 . Additionally, we must provide history functions for those variables in our system (1.6), which have retarded arguments and refer to the past [t 0 − τ m , t 0 ] for m ∈ {R, V, F}. Since we will simulate experiments that start with the initial infection at time t 0 = 0, we must set V (t) = 0 for t < 0. This implies, however, that the choice of the history function for the number of susceptible cells S is irrelevant for the solution of the DDE system (1.6), as the delayed variable S(t − τ m ) with τ m > 0 is always multiplied by V (t − τ m ) and V (t − τ m ) = 0 for t < τ m . For simplicity, we take S(t) = 0 for t < 0.
To solve the DDE system (1.6) we utilize the freely available RADAR5 solver written in ANSI Fortran-90 [S3]. The RADAR5 framework comprises an algorithm to calculate numerically the solution of stiff delay-differential equations based on an adapted 3stage Radau IIA collocation method [S4, S5, S6]. This algorithm corresponds to a certain implicit Runge-Kutta method of order 5 in which the Runge-Kutta equations are solved by means of a suitable Newton process [S7]. The Newton iterations require two Jacobian matrices, first the standard Jacobian matrix J and, second, the Jacobian matrix with respect to the retarded variables J τ . These Jacobian matrices are either computed internally or can be provided by the user. In order to increase the accuracy of the numerical solution, we implemented both time-dependent Jacobian matrices of our DDE system (1.6) within the RADAR5 framework.
In addition to the RADAR5 solver we also tested Matlab's dde23 solver for delaydifferential equations with constant delays [S8, S9]. Although the implementation in Matlab is more convenient and modifications of the model can be achieved with less effort than within the RADAR5 framework, the significantly faster calculation in Fortran makes the RADAR5 framework more suitable for the purpose of parameter estimation, which will be detailed next.

Parameterization of the model based on wildtype DENV data
For the parameterization of the model (1.6), we focus on the time-resolved flow cytometry data set along with ELISA quantification of IFN-λ after infection of the IFIT1deGFP reporter cells with DENV-faR (cf. Fig 6B, C). This data set was obtained by measuring a specified proportion of the cell suspensions and thus provides absolute numbers of virus replicating and IFN responding cells with respect to the measured volume.
At the beginning of the experiments, the IFIT1deGFP reporter cells were infected at a MOI of 0.1 TCID 50 /cell. Based on the definition of MOI given by and the initial number of naïve cells S 0 ∼15000 in the relevant experiments, we set the initial extracellular viral load V 0 (2.7) = 0.1×15000 = 1500 arbitrary units/ml. This speci-fication implies, that the viral load related infection rate r V can be determined up to a scaling factor.
In addition, we take advantage of an independent virus stability experiment to obtain an estimate of the virus degradation rate constant d V (S7C Fig). The DENV-faR stability analysis was performed by incubating IFN-incompetent BHK-21 cells with a high number of DENV-faR particles i.e., the BHK-21 cells were infecting and the cell culture medium was not changed. The extracellular amount of virus was identified by TCID 50 assay over time. Since the production of new infectious virus particles in host cells takes approximately 16 h (cf. Fig 4C and S7B Fig), the measured extracellular amount of virus decreased continuously between 0 and 12 h after incubation. According to the shape of the observed extracellular virus kinetic V(t), we suppose an exponential virus decay and consider as objective function To fit the parameters of the objective function (2.8) to the virus stability data, we utilized the trust-region-reflective least-squares algorithm of Matlab's optimization toolbox [S10]. After applying the optimization method with 10 4 different random initial values, the best fit yielded d V = 0.4/h and provides a satisfying match with the data (cf. S7C Fig, blue curve).
Moreover, we found that fitting of the IFN secretion rate v F and the rate constant of IFN degradation d F at once would lead to a correlation of both parameters, and thus we set To determine the remaining model parameters, we use the time-resolved data set consisting of two independent flow cytometry measurements along with ELISA quantifications post DENV-faR infection (cf. Fig 6B, C) by comparing the following observed kinetics with the related model readouts (ROs): (RO1) The total number of cells A serves as a control and is described in our model by at time t ∈ R and m ∈ {R, V, F}.
(RO2) The number of experimentally identified IFIT1deGFP -DENV-faR double-negative cells are represented byS(t) defined as (2.10) S(t) comprises both susceptible cells and already infected but not yet virus replicating cells I R , which implies that the DENV-faR fluorescence protein is not yet visible in the host cell.
(RO3) The observed DENV-faR positive cells correspond to the calculated number of virus replicating cells I R (t).
(RO5) Since IFN-λ is the most prominent type of IFN in our cell system, we equate the measured level of secreted IFN-λ with the simulated amount of released IFN F (t).
To exploit both independently performed experiments, the model dynamics listed in (RO2)-(RO5) are simultaneously fitted to the related time course data of experiment 1 and 2 by allowing only the initial number of susceptible cells S 0 to vary between experiments, while all other parameters are identical.
As optimization method we use the trust-region-reflective least-squares algorithm of Matlab's optimization toolbox [S10] by calling the Fortran program of our DDE system with a binary Matlab executable (mex) subroutine. Since we compare several components of the model with their respective data, the chi-squares statistic χ 2 of the model parameters θ k , k = 1, . . . , N p , is given by ( 2.11) where N p stands for the number of model parameters and N c denotes the number of model components y j (t i θ 1 , . . . , θ Np ) that describe the associated data set (t i , d j (t i )) with the measurement error σ j (t i ) for i = 1, . . . , N t and j = 1, . . . , N c at N t observation times. To account for the fact, that the data used for fitting were obtained via different measurement techniques and comprise small as well as large values, we assume for the measurement errors a 10% deviation band relative to the measured kinetic data: To find a reliable best fit parameter set of our optimization problem min {θ k k=1,...,N p } χ 2 (θ 1 , . . . , θ N p ) , (2.12) the least-squares minimization has to be repeated for a large number of different random initial values. It turned out that the calculation of the DDE system (1.6) is considerably more time consuming than solving comparable large ODE systems. To speed up calculation, we utilized Matlab's parallel computing toolbox to run several optimization processes simultaneously on a computer cluster. Due to parallelization we were able to repeat the optimization procedure for more than 10 4 different random initial values. The resulting best fit correctly reproduces the observed kinetic data of both experiments (Fig 6C).
In order to assess the uncertainties of the best fit parameter values, we calculate likelihood-based confidence intervals for the estimated parameters by applying the profile-likelihood method [S12]. Likelihood-based confidence regions dependent on a threshold in the likelihoods Q χ 2 (1 − a, DF ), which represents the (1 − a) quantile of the χ 2 -distribution with DF degrees of freedom. As we aim to compute the confidence bound of each individual parameter, the degree of freedom in this case is equal to 1.
The application of the profile-likelihood method for our model regarding the parameters that were estimated based on the kinetic wildtype DENV data set shown in Fig 6C, reveals that all these parameters are identifiable within narrow confidence bounds ( Fig  7A and S1 Table). Since the model is constrained by experimental wildtype DENV data, we can utilize it to make quantitative predictions about the competition between spreading DENV and IFN-induced antiviral protection.
To validate the parameterization, we examined whether key parameter values agree with the results of independent measurements. As described in the main text, we found that the estimated time delays of virus replication (cf. Fig 5) and virus production (cf. Fig 7B) comply with experimental data that were not used for model fitting.
In addition to the time delays, we also investigated the parameterized IFN response in the model using an IFN stimulation experiment, in which A549-IFIT1deGFP reporter cells were treated with IFN-λ (cf. S3B Fig). To imitate the initial conditions of this experiment, we consider a number of susceptible cells S 0 along with the applied extracellular IFN concentration F 0 = 10 ng/ml. When simulating the IFN response, we assume that the recognition of extracellular IFN F turns susceptible cells S into protected cells P with the protection rate r F . The removal of IFN by cellular uptake and, to a lesser extent, by extracellular degradation is taken into account through the IFN degradation rate d F . The dynamics of this small submodel at time t ∈ R ⩾0 is defined by the following ODE system:Ṡ (2.16) We fix S 0 , r F and d F to the fitted values given in S1 Table. The ODE system (2.16) was then solved using a standard ODE solver of MATLAB based on explicit Runge-Kutta method [S13, S14]. The direct comparison shows a good agreement between the predicted fraction of protected cells and the observed IFN-responding cells (Fig 7C). This test indicates that the parameterized model describes the IFN response correctly and, additionally, underscores the accuracy of the estimated protection rate r F .
Taken together, the validity of the model is corroborated by the parameterization within narrow confidence intervals from wildtype DENV data and, additionally, by the consistency with key results of independent experiments that were not considered for model fitting.

Detection of DENV mutant-specific parameters by utilizing the knowledge from wildtype fitting
Driven by the question, which antiviral factors have the greatest influence on viral fitness, we compared wildtype DENV with the E217A DENV mutant, which is unable to modify its RNA genome by 2'-O-methylation [S15]. Our studies show that the E217A mutant elicits faster IFN production resulting in an earlier onset of the IFN response, and barely spreads in IFN-competent cells (Fig 8B, C). In order to analyze the differences between wildtype and attenuated E217A mutant infections on a quantitative level, we utilize the population-based model parameterized exclusively with wildtype DENV data (cf. section 2) to identify E217A mutant-specific parameters.
Motivated by findings from the literature, we argue that four parameters might vary between wildtype DENV and E217A mutant infections. Since studies demonstrate that IFIT1 inhibits translation of mutant RNA [S16, S17], the generation of new infectious virus particles could take longer or occur at a lower rate in the E217A mutant. Moreover, as other reports show that 2'-O-unmethylated RNA is detected more readily by intracellular pattern recognition receptors [S18], infected cells might also produce IFN faster or to a larger extent.
For the parameter optimization regarding the E217A mutant, we therefore take advantage of the established wildtype DENV parameter set (cf. S1 Table) and estimate only the four just mentioned potentially E217A mutant-specific parameters (MP) in line with their respective constrains: (MP1) Delay of mutant virus production τ V mut > τ R (MP2) Mutant virus production rate v V mut ≤ v V (MP3) Delay of IFN secretion after mutant virus infection τ F mut (MP4) IFN secretion rate after mutant virus infection v F mut ≥ v F .
In the same way as in the wildtype DENV case, we determine the E217A mutantspecific parameters (MP1)-(MP4), by fitting the model dynamics listed in (RO2)-(RO5) (cf. section 2 page 5) simultaneously to the related time-resolved data set consisting of two independent flow cytometry measurements along with ELISA quantifications post DENV-faR E217A mutant infection (cf. Fig 8B, C).
We repeated the optimization procedure for more than 5×10 3 different random initial values in the same way as for the data concerning wildtype DENV. The obtained best fit provides a good match with the observed kinetics after DENV-faR E217A mutant infection of both experiments (cf. Fig 8C).
Again, applying the profile-likelihood method, we calculate the 95% confidence intervals of the estimated E217A mutant-specific parameters (MP1)-(MP4) according to (2.13). The resulting narrow confidence bounds confirm the identifiability of all four parameters. More importantly, the comparison of the confidence regions reveals that only two of the four E217A mutant-specific parameters differ strongly from their wildtype DENV values: The time delay of IFN secretion is ∼24 h shorter and the virus production rate is ∼8-fold lower after E217A mutant infection in contrast to wildtype DENV (S2 Table and Fig 8D).

Extension of the population-based model to elucidate the paracrine and autocrine effects of IFN
The prediction of the data-based quantitative model, that a decreased virus production and an accelerated IFN induction make the difference between wildtype DENV and the attenuated E217A mutant, raises the questions how much impact either factor individually has on the outcome of infection, or even if the two factors are possibly related to each other. The latter question is prompted by our observation that stimulation of infected cells with IFN in an early phase after infection causes a reduction of virus replication and consequently virus production (cf. Fig 1A, B and S1 Fig). However, this autocrine effect of IFN must be temporally limited, since stimulation with IFN has no further influence after a certain time period post infection.
To analyze the impact of IFN on infected cells, we expand our established model (1.6) by assuming that the recognition of IFN in a certain time window of IFN responsiveness τ P after viral entry inhibits the production of new infectious virus particles (Fig 9A,  [S19, S20, S21].
At first, we consider a non-negative, continuous random variable T , which represents the waiting time until an infected cell recognizes IFN. The probability P to receive no IFN stimulus within the time spant ∈ R ≥0 is defined as where h denotes the so-called hazard function [S20]. In our case, the hazard function is equal to the time-dependent probability rate of recognizing IFN given by with the rate of IFN-induced inhibition of virus production r P and the amount of extracellular IFN F (t) (Fig 9A, solid green inhibition link).
The distinction whether or not an infected cell is stimulated by IFN within the time window τ P after the time point of infection t I ≥ 0 (Mathematical Modeling Fig 1) is relevant at the time Therefore, we compute the probability that IFN is not recognized up to the time under the condition that there was no stimulation through IFN before the time point of viral entry. The condition is necessary, since sensing of IFN before infection leads to antiviral protection in our model. The probability of an infected cell to produce virus without antiviral IFN action is thus calculated as where F denotes the antiderivative of F .
The probability (3.21) enables us to discriminate between infected cells I VP , which receive no IFN stimulus and thus produce virus with the rate constant v V after τ V , and infected, non-virus producing cells I VP , which recognize IFN within the time window τ P (cf. Mathematical Modeling Fig 1). Accordingly, we include I VP , I VP and F in our DDE model (1.6) and obtain the following DDE system with four time delays: (3.22) The associated initial value problem of the extended model is defined by the DDE system (3.22) for time t ≥ 0 together with initial conditions, which we specify in the following. At the starting time t 0 = 0, only the number of susceptible cells S(0) = S 0 and the extracellular viral load V (0) = V 0 are non-zero. For those variables, that have retarded arguments, we use as history functions the constant zero function, that is, S(t) = 0, V (t) = 0 as well as F (t) = 0 for t < 0.
The full delay-differential equation model (3.22) is numerically solved by applying the RADAR5 solver written in ANSI Fortran-90 [S3]. To obtain a precise solution, we implemented the time-dependent standard Jacobian matrix and the Jacobian matrix with respect to the retarded variables within the RADAR5 framework.
The results of the IFN-stimulation experiment depicted in Fig 1B indicate an IFN responsive time window τ P of approximately 6 h post infection with wildtype DENV at a MOI of 10 TCID 50 /cell. The infection with such a high viral dose causes a synchronized infection of the entire cell population. Since the delay of virus replication after a non-synchronized infection is expected to be delayed, we simulate low dose infections with a time window of IFN responsiveness of τ P = 8 h. Additionally, we fix the rate of IFN-induced inhibition of virus production r P = 10 4 r F , as we suspect that restriction of virus replication is achieved faster than antiviral protection.
According to findings from the literature, IFIT1 is one of those ISGs that can be induced directly after viral recognition in an IFN-independent manner [S22, S23]. An early IFNindependent IFIT1 induction could thus also be contributing to the predicted decreased virus production in IFIT1 reporter cells after E217A mutant infection. To incorporate IFN-independent antiviral effects in the model, we assume that virus producing cells release new infectious virus particles with a reduced rate v V mut = 0.5 v V post DENV mutant infection, which represents half of the virus production rate v V after infection with wildtype DENV.
The full model (3.22) with the mentioned modifications listed in S3 Table matches the observed kinetics of both wildtype and mutant DENV replication and IFN response equally well as the original model (1.6). Moreover, the full model allows us to study the competition between spreading DENV and the antiviral immune response induced in an autocrine as well as paracrine manner by secreted IFN. In addition, we can separately analyze the effect of IFN on infected or naïve cells by setting either r P or r F to zero and subsequently compare the relative importance of both antiviral effects on a quantitative level (cf. main text and Fig 9B).