Ebola virus infection modeling and identifiability problems

The recent outbreaks of Ebola virus (EBOV) infections have underlined the impact of the virus as a major threat for human health. Due to the high biosafety classification of EBOV (level 4), basic research is very limited. Therefore, the development of new avenues of thinking to advance quantitative comprehension of the virus and its interaction with the host cells is urgently needed to tackle this lethal disease. Mathematical modeling of the EBOV dynamics can be instrumental to interpret Ebola infection kinetics on quantitative grounds. To the best of our knowledge, a mathematical modeling approach to unravel the interaction between EBOV and the host cells is still missing. In this paper, a mathematical model based on differential equations is used to represent the basic interactions between EBOV and wild-type Vero cells in vitro. Parameter sets that represent infectivity of pathogens are estimated for EBOV infection and compared with influenza virus infection kinetics. The average infecting time of wild-type Vero cells by EBOV is slower than in influenza infection. Simulation results suggest that the slow infecting time of EBOV could be compensated by its efficient replication. This study reveals several identifiability problems and what kind of experiments are necessary to advance the quantification of EBOV infection. A first mathematical approach of EBOV dynamics and the estimation of standard parameters in viral infections kinetics is the key contribution of this work, paving the way for future modeling works on EBOV infection.

1. Introduction

Ebola was characterized for the first time in 1976 close to the Ebola River located in the Democratic Republic of the Congo (WHO, 1978). Since then, outbreaks of EBOV among humans have appeared sporadically causing lethal diseases in several African countries, mainly in Gabon, South Sudan, Ivory Coast, Uganda, and South Africa (CDC, 2014). Among the most severe symptoms of the EBOV disease are fever, muscle pain, diarrhea, vomiting, abdominal pain and the unexplained hemorrhagic fever (Calain et al., 1999). Fatalities are predominantly associated with uncontrolled viremia and lack of an effective immune response. However, the pathogenesis of the disease is still poorly understood (Peters and Peters, 1999; Feldmann et al., 2003).

Ebola virus belongs to the family of Filoviridae, from Latin filum which means thread (Carter and Saunders, 2013). Ebola virus is classified in Tai Forest, Sudan, Zaire, Reston, and Bundibugyo. The human Ebola epidemics have been mainly related to infection by the Zaire and Sudan strains. Filovirus virions possess several shapes, a property called pleomorphism (Feldmann et al., 2003). These shapes are appearing as either U-shaped, 6-shaped, or other configurations, e.g., Figure 1.

FIGURE 1

Figure 1. Ebola virus molecular structure. The Ebola genome is composed of 3 leader, nucleoprotein (NP), virion protein 35 (VP35), VP40, glycoprotein (GP), VP30, VP24, polymerase (L) protein and 5 trailer (adapted from SIB SWISS Institute of Bioinformatics, 2014).

The natural hosts of EBOV still remain unsettled, but it is tenable that EBOV persists in animals which transmit the virus to non-human primates and humans (Knipe et al., 2001). It has been reported that fruit bats are capable of supporting EBOV replication without becoming ill and may serve as a major reservoir (Swanepoel et al., 1996; Knipe et al., 2001; Leroy et al., 2009; Formenty, 2014). EBOV can spread from an infected person to others through direct contact with blood or body fluids (e.g., saliva, sweat, feces, breast milk, and semen), objects (i.e., needles) that have been contaminated with the virus and infected fruit bats or primates (Peters and Peters, 1999; Feldmann et al., 2003; CDC, 2014). The 2014 Ebola epidemic is the largest ever reported in history, affecting multiple countries in West Africa and being imported to other countries: one infection case was reported in Spain while in the United States one death and two locally acquired cases in healthcare were reported (CDC, 2014).

EBOV can infect a wide variety of cell types including monocytes, macrophages, dendritic cells, endothelial cells, fibroblasts, hepatocytes, adrenal cortical cells, and several types of epithelial cells, all supporting EBOV replication. Monocytes, macrophages, and dendritic cells are early and preferred replication sites of the virus (Knipe et al., 2001). Furthermore, murine studies have revealed that EBOV can infect cells in different compartments, showing high viral titers in liver, spleen, kidney and serum (Mahanty et al., 2003).

Due to its high infectivity and fatality, the virus is classified as a biosafety level-4 agent, restricting basic research for Ebola disease (Halfmann et al., 2008). Infection parameters and quantification of the interactions between the virus and its target cells remain largely unknown. Therefore, the development of new avenues of thinking to bring forward quantitative comprehension of the relationship between the virus and the host is urgently needed. To this end, mathematical models can help to interpret experimental results on quantitative grounds. Model simulations can infer predictions to initiate further and conclusive experiments that may solve relevant scientific questions and advance knowledge of EBOV infection.

Recently, mathematical models have played a central role to capture the dynamics of different virus infections (Nowak and May, 2000). Among the most popular are HIV (Kirschner, 1996; Wu et al., 1998; Duffin and Tullis, 2002; Perelson, 2002; Hernandez-Vargas et al., 2010; Hernandez-Vargas and Middleton, 2013; Jaafoura et al., 2014), hepatitis virus (Ribeiro et al., 2002; Reluga et al., 2009; Guedj et al., 2013) and influenza virus infection models (Baccam et al., 2006; Handel et al., 2010; Smith and Perelson, 2011; Pawelek et al., 2012; Hernandez-Vargas et al., 2014). These models have been instrumental to study the mechanisms that control viral kinetics in order to provide a quantitative understanding and to formulate recommendations for treatments. Similarities of parameter values for EBOV infection to other viral infections that promote outbreaks, e.g., influenza virus infection, could be expected. Nevertheless, to the best of our knowledge, there has not been any mathematical approach until now to describe EBOV dynamics. This and the interaction of EBOV virus with non-human primate epithelial cells is the key contribution of this work.

2. Materials and Methods

2.1. Mathematical Model

The mathematical model proposed here to represent EBOV dynamics is based on the well established target cell-limited model (Nowak and May, 2000), see Figure 2. This has served to model several viral diseases, among them HIV infection (Wu et al., 1998; Perelson, 2002), hepatitis virus infection (Ribeiro et al., 2002) and influenza virus infection (Baccam et al., 2006; Hernandez-Vargas et al., 2014). A detailed reference for modeling of viral dynamics can be found in Nowak and May (2000).

FIGURE 2

Figure 2. Schematic representation of the model for EBOV infection. Target cells (U) are replenished with rate λ and die with rate ρ. Virus (V) infects target cells (U) with rate β. Infected cells are cleared with rate δ. Once cells are productively infected (I), they release virus at rate p and virus particles are cleared with rate c.

Using ordinary differential equations (ODEs), the EBOV infection model is considered as follows:

$\begin{array}{cc}\frac{dU}{dt}=\lambda -\rho U-\beta UV& (1)\end{array}$

$\begin{array}{cc}\frac{dI}{dt}=\beta UV-\delta I& (2)\end{array}$

$\begin{array}{cc}\frac{dV}{dt}=pI-cV& (3)\end{array}$

EBOV target cells can be either in a susceptible (U) or an infected state (I). Cells are replenished with a constant rate λ and die with rate ρ. Note that the condition λ = U₀ρ should be satisfied to guarantee homeostasis in the absence of viral infection, such that only ρ is a parameter to be determined. Virus (V) infects susceptible cells with rate constant β. Infected cells are cleared with rate δ. Once cells are productively infected, they release virus at rate p and virus particles are cleared with rate c.

The initial number of susceptible cells (U₀) can be taken from the experiment in Halfmann et al. (2008) as 5 × 10⁵. The initial value for infected cells (I₀) is set to zero. The viral titer in Halfmann et al. (2008) is measured in foci forming units per milliliter (ffu/ml). The initial viral load (V₀) is estimated from the data using the fractional polynomial model of second order (Royston and Altman, 1994). The best model based on the Akaike Information Criterion (AIC) is presented in Figure 3, providing an estimate of 9 ffu/ml for V₀. The parameter ρ is fixed from literature as 0.001 day⁻¹ (Moehler et al., 2005). The effect of fixing this value on the model output is evaluated with a sensitivity analysis.

FIGURE 3

Figure 3. Data preparation. Fitted statistical model for the wild-type Vero cells infected with EBOV at a low multiplicity of infection (MOI) (Halfmann et al., 2008)

2.2. Experimental Data

As described in the previous section, this paper is mainly focused on the interaction between the virus and the target cells. A safe way to study the virus life cycle was proposed in Halfmann et al. (2008). The disease pathogenesis of EBOV in non-human primates is known to be more faithful in portraying the human condition than in rodents (Knipe et al., 2001). Replication kinetics of EBOV are studied in Vero cells, a cell line derived from kidney epithelial cells of African green monkeys (Halfmann et al., 2008). This non-human primate is a known source of Filoviridae virus infection, e.g., the European Marburg outbreak from 1967 (Knipe et al., 2001). Wild-type Vero cells and a Vero cell line expressing VP30 were tested to reveal their ability to confine EBOV to its complete replication cycle. In this study, viral kinetics for wild-type Vero cells infected with EBOV at different multiplicities of infection (MOI) were considered (Halfmann et al., 2008). The viral growth data is presented in Figure 3. Further details on the data, methods and experiments can be found in Halfmann et al. (2008).

2.3. Parameter Estimation

Parameter fitting is performed minimizing the root mean square (RMS) difference on log scale between the model output, ŷ_i, and the experimental measurement, y_i:

$\begin{array}{cc}\text{RMS}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{({\log }_{10}{y}_i-{\log }_{10}{\widehat{y}}_i)}^2}}& (4)\end{array}$

where n = 5 (Halfmann et al., 2008) is the number of measurements. Differential equations are solved by R 3.1.2 (R Core Team, 2014) using the deSolve package (Soetaert et al., 2010). The minimization of RMS is performed using the Differential Evolution (DE) algorithm employing the DEoptim package (Storn, 1997; Mullen et al., 2011). The DE global optimization algorithm does not rely on initial parameter guesses and converged faster than the other tested methods, including genetic algorithms and the quasi-Newton (BFGS, L-BFGS-B) algorithms.

2.4. Parameter Uncertainty

Viral load variability is very large for several viral infectious diseases (Mahanty et al., 2003; Baccam et al., 2006; Toapanta and Ross, 2009; Groseth et al., 2012). In order to consider the large variability of biological problems, a bootstrap method is applied to the data series presented in Halfmann et al. (2008). Bootstrapping is a statistic method for assigning measures of accuracy to estimates (Davison and Hinkley, 1997; Xue et al., 2010). The nonparametric bootstrap requires data to be independent and identically distributed while the parametric bootstrap requires to impose on the data a distribution assumption which is usually unknown. For the data in Halfmann et al. (2008), three bootstrap approaches were considered: (i) the conventional parametric approach assumes a log-normal distribution of the measurement, (ii) the nonparametric approach assumes uniform distribution in the measurement range, and (iii) the weighted bootstrap assigns to the cost function a vector of random weights from exponential distribution with mean one and variance one (Ma and Kosorok, 2005; Xue et al., 2010).

For each repetition, the model parameters are refitted to obtain the corresponding parameter distribution. The 95% confidence interval of parameter estimates is computed using the outcome of the bootstrap method (Xue et al., 2010). For each parameter, the 2.5 and 97.5% quantiles of the estimates are used to form the 95% confidence interval.

2.5. Parameter Identifiability and Sensitivity

A critical obstacle to overcome in mathematical modeling is how to verify whether model parameters are identifiable based on the measurements of output variables (Xia, 2003; Xia and Moog, 2003; Wu et al., 2008; Miao et al., 2011). A system that is algebraically identifiable may still be practically non-identifiable if the amount and quality of the measurements is insufficient and the data shows large deviations. The novel approach proposed in Raue et al. (2009) exploits the profile likelihood to determine identifiability and is considered here. This method is able to detect both structurally and practically non-identifiable parameters.

Identifiability properties are studied for the model Equations (1–3) and the data set in Halfmann et al. (2008). The idea behind this approach is to explore the parameter space for each parameter θ_i by re-optimizing the RMS with respect to all other parameters θ_{j ≠ i}. In particular, for each parameter θ_i, a wide range of values centered at the optimized value is generated in an adaptive manner. Re-optimization of RMS with respect to the other parameters is done for each value of parameter θ_i. The main task is to detect directions where the likelihood flattens out (Raue et al., 2009). The resulting profiles are plotted vs. each parameter range to assess the parameter identifiability visually.

In model fitting, some parameters may have little effect on the model outcome, while other parameters are so closely related that simultaneous fitting could be a difficult task. For this aspect, the scatter plots using pairs of parameters over different bootstrap replicates will be reported. Furthermore, sensitivity analysis of the estimated parameters is performed (Brun et al., 2001; Soetaert, 2014). For each data point the derivative of the corresponding modeled variable value with respect to the selected parameter is computed. The normalized sensitivity function reads as

$\begin{array}{cc}\frac{\partial y_i}{\partial {\Theta }_j}·\frac{w_{{\Theta }_j}}{w_{y_i}}& (5)\end{array}$

where y_i denotes the model variables, Θ_j is the parameter of interest, and the ratio w_Θj/w_yi is the normalized factor corresponding to its nominal value (Soetaert and Petzoldt, 2010). Summary statistics of the sensitivity functions can be used to qualify the impact of the parameter on the output variables, i.e., the higher the absolute value of the sensitivity summary statistics, the more important the parameter (Brun et al., 2001). For the model in Equations (1–3), the sensitivity functions will be plotted vs. time to illustrate the parameters' role on the model output. The parameters that have little effect do not need to be fine-tuned extensively in model fitting.

2.6. Cross-Validation

It is important to prove how the model predictions will generalize to an independent data set, revealing how accurately the predictive value of a model is in practice. In this paper, the parameter set obtained from the data of wild-type Vero cells infected at low MOI is used to predict the replication kinetics of the data at high MOI presented in Halfmann et al. (2008).

3. Results

Although significant progress has been made to the identification and characterization of EBOV, human data is very limited due to the long asymptomatic periods of the virus and its high mortality. Animal models are pivotal to shed light on this lethal disease. Due to the very close similarities with the human immune system, non-human primates are the preferred animal model for several viral infections e.g., HIV). Moreover, EBOV infection has been adapted to guinea pigs and mice (Feldmann et al., 2003), serving as a flexible model in comparison to human and non-human primates. In this work, we focus on the interaction between the virus and the host cells. In vitro data can be very convenient due to the important simplification of the in vivo complexity of biological problems. Thus, for parameter fitting procedures, we consider the experimental data from Halfmann et al. (2008), which investigates EBOV kinetics in a Vero cell line.

Before rigorous optimization methods can be applied to estimate the model parameters using experimental data, the verification of parameter identifiability is required. The omission of identifiability analyses may result in incorrect fits and consequently incorrect interpretations. The identifiability analysis in the model Equations (1–3) has been broadly studied (Xia, 2003; Xia and Moog, 2003; Wu et al., 2008; Miao et al., 2011; Hernandez-Vargas et al., 2014). All parameters in the model Equations (1–3) were shown to be algebraically identifiable given measurements of viral load and initial conditions (U₀, I₀, and V₀) (Wu et al., 2008). However, the difference between structural identifiability and practical identifiability in the presence of measurement error requires further identifiability studies. To address practical identifiability, the approach proposed by Raue et al. (2009) is considered here for the data presented in Figure 3.

The resulting RMS profiles in Figure 4 for β, p and c show a convex shape of which the optimization routine can reach their minimum. Note that the profile of δ is flat in one tail, suggesting that parameter δ can be chosen arbitrarily small without affecting the fit quality (Raue et al., 2009). In spite of this, the lower bound of this parameter has a clear biological constraint. To be precise, the half-life of an infected cell cannot be longer than that of an uninfected cell. There is experimental evidence that the half-life of epithelium cells in lung is 17–18 months in average (Rawlins and Hogan, 2008). In view of this, the infected cell death rate (δ) is fixed at 10⁻³.

FIGURE 4

Figure 4. Parameter Identifiability. RMS profile of model parameters. Each parameter is varied in a wide range around the optimized value. Subsequently, the DE algorithm is used to refit the remaining parameters to the data set of Halfmann et al. (2008). The vertical dashed lines indicate the value obtained from the optimization for all four parameters collectively.

Bootstrapping can provide more insights into the distribution of parameter values based on experimental data in Halfmann et al. (2008). For the sake of clarity, we present only the weighted bootstrap (Xue et al., 2010) in the results, the other two methods can be found in the supplementary material. Distributions of the model parameters are shown in Figure 5. Bootstrap estimates for the viral clearance (median c = 1.05 day⁻¹) is slightly below other viral infection results (Table 1). For example, clearance of influenza virus varied from 2.6 to 15 day⁻¹ in (Baccam et al., 2006; Miao et al., 2011; Pawelek et al., 2012; Hernandez-Vargas et al., 2014). This may be attributed to the fact that the viral clearance is computed for in vitro experiments.

FIGURE 5

Figure 5. Weighted bootstrap results. Top row: Distributions from 1000 sample estimates are presented for the three parameters: β, p and c. Bottom row: Scatter plot between bootstrap parameters. The parameter ρ is fixed during the bootstrapping at 0.001 (Moehler et al., 2005). Numerical values for the model Equations (1–3) are presented in the Table 1.

TABLE 1

Table 1. Estimates of infection parameters*.

EBOV is known to replicate at an unusually high rate that overwhelms the protein synthesis of infected cells (Sanchez, 2001). Consistent with this observation, bootstrap estimates revealed a very high rate of viral replication, p = 62 (95% CI: 31−580) (Table 1). Although the scatter plot in Figure 5 shows that the estimate of p can be decreased given a higher effective infection rate (β), a replication rate of at least 31.8 ffu/ml cell⁻¹ day⁻¹ is still needed to achieve a good fit of the viral replication kinetics in Figure 3.

Scatter plots are a graphical sensitivity analysis method, and a simple but useful tool to test the robustness of the results. The estimated parameters are plotted against each other. Scatter plots for the parameters in Figure 5 provide visual evidence that these parameters strongly depend on one another such that their individual values can not be independently determined. That is, increasing the values of p increases the estimations of c. Decreasing the estimations of β increases the estimation of both c and p. However, the green curves in Figure 5 provide the most likely region where the parameters values can be found. In order to verify this intuition, we fix the viral clearance rate (c) at 4.2 (Miao et al., 2010) and then estimate the others two parameters (β and p). The results of 1000 bootstrap replicates reveal that fixing the parameter c improves the fitting with a narrow confidence interval (see Supplementary Materials 1.3).

The sensitivity study for the mathematical model Equations (1–3) is performed in a similar fashion to Brun et al. (2001); Soetaert (2014). Figures 6A–E show the effect on the viral load when varying the respective parameter by 10, 20 and 50% around its nominal value. It can be seen that the healthy cell death rate (ρ), which in the virus-free steady state represents the cell turnover, has little effect on the viral load kinetics. This can be attributed to the fact that the experiment was performed in vitro and within a short period. Similarly, the effect of the infected cell death rate (δ) can also be neglected. This could be explained by the fact that the observed Ebola viral load was not decreasing (Figure 3), contrary to observations in other viral infections, e.g., influenza virus (Baccam et al., 2006). The remaining three parameters (β, p, and c) are sensitive, in the sense that a small change in parameter value can lead to a large difference in viral kinetics. Figure 6F summarizes in detail the parameter sensitivity functions. It is clear that the three parameters β, p, and c govern the infection kinetics while the effect of the two parameters ρ and δ can be neglected for this data set. Therefore, fixing both ρ and δ is adequate for the presented problem.

FIGURE 6

Figure 6. Sensitivity of parameters. (A–E) Plotting of viral titer variation vs. time. The dashed line is the viral kinetics obtain from nominal parameter values. Three color shades in each figure represent the viral load variation range when varying the corresponding parameter by a percentage denoted in the legend. (F) Parameters sensitivity function over time, the values in y-axis are calculated using Equation (5).

Moreover, both β and p can be seen as consistently increasing the viral load because their respective sensitivity functions are always positive, in contrast to the parameter c. Note that the absolute magnitude of change in the sensitivity functions of these three parameters is approximately equal over time (Figure 6F). The strong similarity in the sensitivity functions indicates that the corresponding parameters have equivalent effect on the viral titer. For instance, the sensitivity functions of β and p are very similar so that almost the same output of viral titer will be generated by increasing β if p is decreased correspondingly. A similar statement can also be made about the relationship between c and β.

Computational simulations for the best fitting of the proposed mathematical model Equations (1–3) plotted in Figure 7B show that the virus grows exponentially from day 1 to 5 post infection. This is consistent with the mathematical analysis developed in Nowak et al. (1996), which deduced that the virus initially grows exponentially and can be better modeled as exp(r₀t) while the susceptible cell population remains relatively constant, where r₀ is the leading eigenvalue which solves the equation r²₀ + (c + δ)r₀ − (βpU₀ − cδ) = 0.

FIGURE 7

Figure 7. Model fitting for EBOV kinetics. Viral titer data with low MOI from Halfmann et al. (2008) and simulations from the best fit shown in Table 1 are in panel (A) for the host cells and (B) for the viral titer.

Viral titer peaks at high levels, more than 10⁷ ffu/ml, which in general is 10 fold higher than those reported in influenza virus infection (Toapanta and Ross, 2009; Hernandez-Vargas et al., 2014). In addition, the viral titer reaches a plateau at day 6 and may remain at those levels (Figure 7B). No depletion of infected cells is observed in the period of observation. This could be a combined effect attributed to either high infection rate or high replication rate, and to the slow clearance of infected cells. To achieve virus titer levels as reported in Halfmann et al. (2008), either a high infection rate (β) of susceptible cells, or a high replication rate is required (Figure 5). Note that even though these estimations were performed in vitro, in vivo murine studies for EBOV infection (Mahanty et al., 2003) showed similar kinetics and time scales as those presented in Figure 7B.

3.1. Transmission Measures

Infectivity is a critical parameter to assess the ability of a pathogen to establish an infection (Diekmann et al., 1990). To determine infectivity, we compute the reproductive number (R₀), which is defined as the expected number of secondary infections produced by an infected cell in its lifetime (Diekmann et al., 1990; Heffernan et al., 2005). On the one hand, if R₀ is less than one, each infected individual produces on average less than one infected individual, and therefore the infection will be cleared from the population. On the other hand, if R₀ is greater than one, the pathogen is able to invade the susceptible population. This epidemiological concept can be applied to the model Equations (1–3) and computed as follows (Nowak et al., 1996):

$\begin{array}{cc}{R}_0=\frac{\lambda p\beta }{c\rho \delta }& (6)\end{array}$

As expected, the estimated reproductive number in EBOV infection is very high, see Figure 8A and numerical results in Table 1. These results can be attributed to the fact that no depletion of virus was observed and to a slow clearance of infected cells. Thus, both parameters δ and c increase the value of R₀. Note that very high estimates of the reproductive number in highly viremic influenza virus strains from in vitro experiments have also been reported, with an average of 1.7 × 10³ (Pinilla et al., 2012). It is worth to mention that fitting the model to in vitro data in Halfmann et al. (2008) could lead to small estimates for c and δ in comparison to an in vivo situation. Nevertheless, estimates of the epithelial cell half-life were 6 months in the trachea and 17 months in the lungs in average (Bowden, 1983; Rawlins and Hogan, 2008), which corresponds to a δ equal to 0.003 and 0.001, respectively. As mentioned previously, the δ was fixed at 0.001 in the computation of R₀. Therefore, the estimated values of R₀ interval are very likely to be positioned in a biologically plausible range, especially the upper bound. Notwithstanding, the estimate of R₀ presented here should be interpreted with care within the limits of the data used.

FIGURE 8

Figure 8. Transmission measures. Bootstrap estimate of (A) reproductive number and (B) infecting time in hours. Numerical values can be found in Table 1.

Recent viral modeling works (Holder et al., 2011; Pinilla et al., 2012) have also introduced the term infecting time, which represents the amount of time required for a single infectious cell to cause the infection of one more cell within a completely susceptible population. Strains with a shorter infecting time have a higher infectivity (Holder et al., 2011; Pinilla et al., 2012). From model Equations (1–3), this measure can be computed as follows:

$\begin{array}{cc}{t}_{\text{inf}}=\sqrt{\frac{2}{p\beta U_0}}& (7)\end{array}$

Bootstrap results showed that EBOV possesses an average infecting times of 9.49 h (Table 1) which is approximately 7 times slower than the infecting time of influenza virus (Holder et al., 2011). This number provides a reasonable explanation for the kinetics of susceptible cells which slowly decrease from day 1 to day 4 (Figure 7A), and quickly deplete within the last 2 days. This number could also explain the absence of viral replication within the first 5.6 h after infection. This period corresponds to the short decreasing period observed in Figure 7B. The initial decrease of viral load thus can be attributed to self-clearance of the virus when some viruses have infected cells but are not yet able to replicate.

The infectivity parameters in Figure 8 characterize the EBOV infection kinetics in the data in Halfmann et al. (2008). The slow infection time of EBOV is compensated by its efficient replication. As a result, a short delay is followed by a massive amount of virus. The above infectivity parameters contributed an explanation for the high levels of viral load even when the susceptible cells were already depleted at the end of the experiment.

The best set of estimated parameters is challenged to validate the data at high multiplicities of infection (MOI) in Halfmann et al. (2008). The initial viral load is estimated using the fractional polynomial model of second order providing V₀ at 460 ffu/ml. Figure 9 shows that the parameters derived from data at low multiplicity of infection are still consistent with data generated at high multiplicity of infection. The predicted kinetics follows the experimental data closely when changing the initial condition of the viral titer to 50 folds higher.

FIGURE 9

Figure 9. Cross-validation. Test of estimated parameters on an independent set of data. The viral replication kinetics in wild-type Vero cells infected with EBOV at a high multiplicities of infection (MOI) in Halfmann et al. (2008) are modeled starting from a higher initial viral load of V₀ = 460 ffu/ml. The (Mean) indicates the predicted kinetics using parameters obtained from bootstrap while (Best) refers to the predicted kinetics using the parameters resulting from the optimization.

4. Discussion

Ebola virus (EBOV) is highly pathogenic for humans, being nowadays one of the most lethal pathogens worldwide. Ebola fatalities are predominantly associated with uncontrolled viremia and lack of an effective immune response (i.e., low levels of antibodies and no cellular infiltrates at sites of infection) (Feldmann et al., 2003).

The work presented here focused on the interaction between EBOV and the host cells, i.e., epithelial cells of green monkey. Experimental data on the Vero cell line from non-human primates could help to better understand the virus infection dynamics in humans (Knipe et al., 2001). However, the in vitro studies must be translated carefully to avoid over-interpretation to the in vivo context, which can sometimes lead to erroneous conclusions. Especially, the EBOV infection has been known to have abnormal behavior in vivo where different cells types and the immune system are involved (Knipe et al., 2001). Additionally, given the fact that EBOV exhibits an asymptomatic period in humans (Leroy et al., 2000), the viral dynamics model in vivo should take the eclipse phase into consideration. This feature can be modeled by adding an appropriate eclipse phase term as has been done previously (Moehler et al., 2005; Baccam et al., 2006). Nevertheless, given the problem of parameter identifiability exposed in the results, a complex model would not bring any better understanding. Once more data would become available, future work could attempt to address this issue, especially in the in vivo context.

The exposed identifiability issues in the results reveal the problematic of parameter estimation using solemnly the viral load measurements. Here, our efforts to cope thoroughly with the identifiability issues spotted the current restrictions on the estimated parameters. These restrictions cannot be resolved without the progress of new experiments, more measurements are necessary to sort out the identifiability problems presented here, e.g., measurements of infected and non-infected cells. Another possible experiment is to determine the EBOV clearance rate in the absence of target cells. For instance, Pinilla et al. (2012) employed an experiment in a similar fashion to determine the viral infectivity loss (c). Known influenza virus titers were incubated without target cells and followed up to determine the remaining infectious titers (Pinilla et al., 2012). In this way the approximate values of the viral clearance rate could be determined and provide a more accurate estimates for the whole set of kinetics parameters, as shown in the Supplemental Material 1.3.

The high EBOV replication reported here is in agreement with recent findings by Misasi and Sullivan (2014) as well as documented in Knipe et al. (2001), reporting that early and coordinated disruptions by Ebola genes and proteins (VP24, VP30, and VP35) lead to elevated levels of virus replication. The bootstrap results suggested that the EBOV average infecting time is approximately 9.5 h, at least 5 fold slower than estimations from influenza virus infection (Pinilla et al., 2012). These simulations outline the EBOV kinetics in the data from Halfmann et al. (2008), suggesting that a slow infecting time of EBOV is compensated by its efficient replication.

The model results suggested that the saturation of viral growth as observed in the data is induced by the loss of susceptible cells. This result has to be re-evaluated with a more complete data set, as the present data set would also be appropriately described by a logistic-growth model (data not shown) with an unspecific limitation of resources. However, a logistic model can explain only the growth behavior of the virus. As pointed out before (Wu et al., 2008), a higher resolution of the data and later time points which exhibit the long-term behavior of the viral load are required for a full determination of the mechanisms at work.

EBOV infection from in vitro and even murine systems may differ considerably from humans. The latency phase in human is much longer than in animals and EBOV symptoms in humans may appear from 2 to 21 days after exposure to the virus, having an average time of 8–10 days (Peters and Peters, 1999). Remarkably, mice infected by intra-peritoneal injection develop symptomatic infection where EBOV will increase rapidly at day 4 and continue to increase until day 6, with death occurring at day 6–7 post-infection (Mahanty et al., 2003). These experimental observations are compatible with our simulation results, suggesting that the growth of infected cells starts at day 3 post infection (Figure 7) while almost the whole susceptible cell pool is depleted at day 6 post infection. It is worth to mention that EBOV kinetics were similar in different tissue compartments (Mahanty et al., 2003): liver, spleen, kidney and serum. Consequently, further modeling approaches should address the EBOV kinetics in different compartments of the infected host.

The in vitro system may mimick a human context where the immune response against EBOV is not working adequately. The onset of a CD8+ T cell response as well as of the antibody response (Gupta et al., 2003) rely on early regulation of cytokines in the asymptomatic phase of the disease (Mahanty et al., 2003; Ebihara et al., 2006; García-Sastre and Biron, 2006). Human EBOV infection revealed that patients infected by the Sudan strain had lower levels of tumor necrosis factor TNF-α and interferon IFN-γ compared to those found in patients with fatal Zaire strain infection (Hutchinson and Rollin, 2007). Additionally, the levels of IFN-α were found significantly higher in surviving patients with Sudan strain infection (Hutchinson and Rollin, 2007), whereas the levels of IL-6, IL-8, IL-10, and macrophage inflammatory proteins were higher in patients with fatal infections (Hutchinson and Rollin, 2007). Therefore, modeling the effects of IFN-I would limit the number of infected cells by the introduction of a resistant state with a possible impact on the value of the viral replication rate (p). Future modeling studies need to quantify the situation in vivo where the effect of the immune system is taken into account.

The modeling work developed in this paper paves the way for future mathematical models and experiments to shed light on the reasons for less efficient control of Ebola virus infections. Determining empirically the EBOV clearance rate in the absence of target cells would fulfill the picture of EBOV kinetics in vitro. In addition, due to the critical relevance of the cytokine effects in EBOV pathogenesis, future modeling attempts should be directed to establish a more detailed model of interactions between the relevant cytokines and EBOV. Further insights into immunology and pathogenesis of EBOV will help to improve the outcome of this lethal disease.

Conflict of Interest Statement

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.

Acknowledgments

This work was supported by iMed—the Helmholtz Initiative on Personalized Medicine. In addition, we thank the support provided by the Measures for the Establishment of Systems Medicine (e:Med) projects in Systems Immunology and Image Mining in Translational Bio- marker Research (SYSIMIT) and in identification of predictive response and resistance factors to targeted therapy in gastric cancer using a systems medicine approach (SYS-Stomach) by the Federal Ministry of Education and Research (BMBF), Germany. VKN has been supported by the President's Initiative and Networking Funds of the Helmholtz Association of German Research Centres (HGF) under contract number VH-GS-202.

Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fmicb.2015.00257/abstract

References

Baccam, P., Beauchemin, C., Macken, C. A., Hayden, F. G., and Perelson, A. S. (2006). Kinetics of influenza A virus infection in humans. J. Virol. 80, 7590–7599. doi: 10.1128/JVI.01623-05

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Bowden, D. (1983). Cell turnover in the lung. Am. Rev. Res. Dis. 128(2 Pt 2), S46–S48.

PubMed Abstract | Full Text | Google Scholar

Brun, R., Reichert, P., and Kuensch, H. R. (2001). Practical identifiability analysis of large environmental simulation models. Water Res. Res. 37, 1015–1030. doi: 10.1029/2000WR900350

CrossRef Full Text | Google Scholar

Calain, P., Bwaka, M. A., Colebunders, R., Roo, A. D., Guimard, Y., Katwiki, K. R., et al (1999). Ebola hemorrhagic fever in kikwit, democratic republic of the congo : clinical Observations in 103 Patients. J. Infect. Dis. 179(Suppl. 1), 1–7. doi: 10.1086/514308

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Carter, J., and Saunders, V. (2013). Virology: Principle and Applications. Chichester, UK: John Wiley.

Google Scholar

CDC. (2014). CDC Report to Ebola Virus Disease 2014. Technical report.

Davison, A. C., and Hinkley, D. V. (1997). Bootstrap Methods and Their Application. Cambridge, UK: Cambridge University.

Google Scholar

Diekmann, O., Heesterbeek, J., and Metz, J. (1990). On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases. J. Math. Biol. 28, 365–382. doi: 10.1007/BF00178324

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Duffin, R. P., and Tullis, R. H. (2002). Mathematical models of the complete course of HIV infection and AIDS. J. Theor. Med. 4, 215–221. doi: 10.1080/1027366021000051772

CrossRef Full Text | Google Scholar

Ebihara, H., Takada, A., Kobasa, D., Jones, S., Neumann, G., Theriault, S., et al. (2006). Molecular determinants of Ebola virus virulence in mice. PLoS Pathog. 2:e73. doi: 10.1371/journal.ppat.0020073

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Feldmann, H., Jones, S., Klenk, H.-d., and Schnittler, H.-j. (2003). Ebola virus: from discovery to vaccine. Nat. Rev. 3, 677–685. doi: 10.1038/nri1154

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Formenty, P. (2014). “Chapter 9 - ebola virus disease,” in Emerging Infectious Diseases, eds N. Ergnl, F. Can, L. Madoff, and M. Akova (Amsterdam: Academic Press), 121–134.

García-Sastre, A., and Biron, C. a. (2006). Type 1 interferons and the virus-host relationship: a lesson in détente. Science 312, 879–882. doi: 10.1126/science.1125676

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Groseth, A., Marzi, A., Hoenen, T., Herwig, A., Gardner, D., Becker, S., et al. (2012). The Ebola virus glycoprotein contributes to but is not sufficient for virulence in vivo. PLoS Pathog. 8:e1002847. doi: 10.1371/journal.ppat.1002847

PubMed Abstract | Full Text | 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, 3–8. doi: 10.1073/pnas.1203110110

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Gupta, M., Mahanty, S., Greer, P., Towner, J. S., Shieh, W.-J., Zaki, S. R., et al. (2003). Persistent Infection with Ebola Virus under Conditions of Partial Immunity. J. Virol. 78, 958–967. doi: 10.1128/JVI.78.2.958-967.2004

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Halfmann, P., Kim, J., Ebihara, H., Noda, T., Neumann, G., Feldmann, H., et al. (2008). Generation of biologically contained Ebola viruses. Proc. Natl. Acad. Sci. U.S.A. 105, 1129–1133. doi: 10.1073/pnas.0708057105

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Handel, A., Longini, I. M., and Antia, R. (2010). Towards a quantitative understanding of the within-host dynamics of influenza A infections. J. R. Soc. Interf. 7, 35–47. doi: 10.1098/rsif.2009.0067

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Heffernan, J. M., Smith, R. J., and Wahl, L. M. (2005). Perspectives on the basic reproductive ratio. J. R. Soc. Interf. 2, 281–293. doi: 10.1098/rsif.2005.0042

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hernandez-Vargas, E. A., Colaneri, P., Middleton, R. H., and Blanchini, F. (2010). Discrete-time control for switched positive systems with application to mitigating viral escape. Int. J. Rob. Nonlin. Control 21, 1093–1111. doi: 10.1002/rnc.1628

CrossRef Full Text | Google Scholar

Hernandez-Vargas, E. A., and Middleton, R. H. (2013). Modeling the three stages in HIV infection. J. Theor. Biol. 320, 33–40. doi: 10.1016/j.jtbi.2012.11.028

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hernandez-Vargas, E. a., Wilk, E., Canini, L., Toapanta, F. R., Binder, S. C., Uvarovskii, A., et al. (2014). Effects of aging on influenza virus infection dynamics. J. Virol. 88, 4123–4131. doi: 10.1128/JVI.03644-13

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Holder, B. P., Simon, P., Liao, L. E., Abed, Y., Bouhy, X., Beauchemin, C. A. A., et al. (2011). Assessing the in vitro fitness of an oseltamivir-resistant seasonal A/H1N1 influenza strain using a mathematical model. PLoS ONE 6:e14767. doi: 10.1371/journal.pone.0014767

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hutchinson, K. L., and Rollin, P. E. (2007). Cytokine and chemokine expression in humans infected with Sudan Ebola virus. J. Infect. Dis. 196(Suppl.), S357–S363. doi: 10.1086/520611

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Jaafoura, S., de Goër de Herve, M. G., Hernandez-Vargas, E. A., Hendel-Chavez, H., Abdoh, M., Mateo, M. C., et al. (2014). Progressive contraction of the latent HIV reservoir around a core of less-differentiated CD4(+) memory T Cells. Nat. Commun. 5, 1–8. doi: 10.1038/ncomms6407

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kirschner, D. (1996). Using mathematics to understand HIV immune dynamics. AMS Notices 43, 191–202.

Google Scholar

Knipe, D. M., Howley, P. M., Griffin, D. E., Lamb, R. A., Martin, M. A., Roizman, B., et al. (2001). Field Virology. Philadelphia, PA: Lippincott Williams and Wilkins.

Leroy, E. M., Baize, S., Volchkov, V. E., Fisher-Hoch, S. P., Georges-Courbot, M. C., et al. (2000). Human asymptomatic Ebola infection and strong inflammatory response. Lancet 355, 2210–2215. doi: 10.1016/S0140-6736(00)02405-3

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Leroy, E. M., Epelboin, A., Mondonge, V., Pourrut, X., Gonzalez, J. P., Muyembe-Tamfum, J. J., et al. (2009). Human Ebola outbreak resulting from direct exposure to fruit bats in Luebo, Democratic Republic of Congo, 2007. Vector Borne Zoonotic Dis. 9, 723–728. doi: 10.1089/vbz.2008.0167

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ma, S., and Kosorok, M. R. (2005). Robust semiparametric m-estimation and the weighted bootstrap. J. Multivar. Anal. 96, 190–217. doi: 10.1016/j.jmva.2004.09.008

CrossRef Full Text | Google Scholar

Mahanty, S., Gupta, M., Paragas, J., and Bray, M. (2003). Protection from lethal infection is determined by innate immune responses in a mouse model of Ebola virus infection. Virology 312, 415–424. doi: 10.1016/S0042-6822(03)00233-2

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Miao, H., Hollenbaugh, J. A., Zand, M. S., Holden-Wiltse, J., Mosmann, T. R., Perelson, A. S., et al. (2010). Quantifying the early immune response and adaptive immune response kinetics in mice infected with influenza A virus. J. Virol. 84, 6687–6698. doi: 10.1128/JVI.00266-10

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Miao, H., Xia, X., Perelson, A. S., and Wu, H. (2011). on Identifiability of Nonlinear Ode Models and Applications in Viral Dynamics. SIAM Rev. Soc. Industr. Appl. Math. 53, 3–39. doi: 10.1137/090757009

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Misasi, J., and Sullivan, N. J. (2014). Camouflage and misdirection: The full-on assault of ebola virus disease. Cell 159, 477–486. doi: 10.1016/j.cell.2014.10.006

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Moehler, L., Flockerzi, D., Sann, H., and Reichl, U. (2005). Mathematical model of influenza a virus production in large-scale microcarrier culture. Biotechnol. Bioeng. 90, 46–58. doi: 10.1002/bit.20363

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mullen, K., Ardia, D., Gil, D., Windover, D., and Cline, J. (2011). DEoptim: An R package for global optimization by differential evolution. J. Stat. Softw. 40, 1–26. Available online at: http://www.jstatsoft.org/v40/i06/bibtex

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 | Full Text | CrossRef Full Text | Google Scholar

Nowak, M. A., and May, R. (2000). Virus Dynamics: Mathematical Principles of Immunology and Virology, Vol. 291. Oxford, UK: Oxford University Press.

Pawelek, K. A., Huynh, G. T., Quinlivan, M., Cullinane, A., Rong, L., and Perelson, A. S. (2012). Modeling within-host dynamics of influenza virus infection including immune responses. PLoS Comp. Biol. 8:e1002588. doi: 10.1371/journal.pcbi.1002588

PubMed Abstract | Full Text | 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 | Full Text | CrossRef Full Text | Google Scholar

Peters, C., and Peters, J. (1999). An introduction to Ebola: the virus and the disease. J. Infect. Dis. 179, ix–xvi. doi: 10.1086/514322

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Pinilla, L. T., Holder, B. P., Abed, Y., Boivin, G., and Beauchemin, C. A. (2012). The H275Y neuraminidase mutation of the pandemic A/H1N1 influenza virus lengthens the eclipse phase and reduces viral output of infected cells, potentially compromising fitness in ferrets. J. Virol. 86, 10651–10660. doi: 10.1128/JVI.07244-11

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

R Core Team (2014). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

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 | Full Text | CrossRef Full Text | Google Scholar

Rawlins, E. L., and Hogan, B. L. M. (2008). Ciliated epithelial cell lifespan in the mouse trachea and lung. Am. J. Physiol. Lung Cell. Mol. Physiol. 295, L231–L234. doi: 10.1152/ajplung.90209.2008

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Reluga, T., Dahari, H., and Perelson, A. (2009). Analysis of hepatitis C virus infection models with hepatocyte homeostasis. SIAM J. Appl. Math. 69, 999–1023. doi: 10.1137/080714579

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ribeiro, R. M., Lo, A., and Perelson, A. S. (2002). Dynamics of hepatitis B virus infection. Microbes Infect. 4, 829–835. doi: 10.1016/S1286-4579(02)01603-9

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Royston, P., (Imperial College School of Medicine, London) Altman, D. G., and (Imperial Cancer Research Fund, London). (1994). Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling. Appl. Statist. 43, 429–467.

Sanchez, A. (2001). Filoviridae: Marburg and Ebola Viruses. Philadelphia, PA: John Wiley & Sons, Ltd.

SIB SWISS Institute of Bioinformatics (2014). Ebolavirus. Genève.

Smith, A. M., and Perelson, A. S. (2011). Influenza A virus infection kinetics: quantitative data and models. Syst. Biol. Med. 3, 429–445. doi: 10.1002/wsbm.129

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Soetaert, K. (2014). Package rootsolve: roots, gradients and steady-states in r.

Google Scholar

Soetaert, K., and Petzoldt, T. (2010). Inverse modelling, sensitivity and monte carlo analysis in R using package FME. J. Stat. Softw. 33, 1–28. Available online at: http://www.jstatsoft.org/v33/i03/bibtex

Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010). Solving differential equations in r: Package desolve. J. Stat. Softw. 33, 1–25. Available online at: http://www.jstatsoft.org/v33/i09/bibtex

Storn, R. (1997). Differential evolution a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optimizat. 11, 341–359. doi: 10.1023/A:1008202821328

CrossRef Full Text | Google Scholar

Swanepoel, R., Leman, P. A., Burt, F. J., Zachariades, N. A., Braack, L. E., Ksiazek, T. G., et al. (1996). Experimental inoculation of plants and animals with Ebola virus. Emerg. Infect. Dis. 2, 321–325. doi: 10.3201/eid0204.960407

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Toapanta, F. R., and Ross, T. M. (2009). Impaired immune responses in the lungs of aged mice following influenza infection. Res. Res. 10, 1–19. doi: 10.1186/1465-9921-10-112

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

WHO (1978). Ebola Haemorrhagic Fever in Sudan, 1976: Report of a World Health Organization International Study Team. Technical report.

Google Scholar

Wu, H., Ding, A. A., and De Gruttola, V. (1998). Estimation of HIV dynamic parameters. Stat. Med. 17, 2463–2485. doi: 10.1002/(SICI)1097-0258(19981115)17:21<2463::AID-SIM939>3.0.CO;2-A

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wu, H., Zhu, H., Miao, H., and Perelson, A. S. (2008). Parameter identifiability and estimation of HIV/AIDS dynamic models. Bull. Math Biol. 70, 785–799. doi: 10.1007/s11538-007-9279-9

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Xia, X. (2003). Estimation of HIV/AIDS parameters. Automatica 39, 1983–1988. doi: 10.1016/S0005-1098(03)00220-6

CrossRef Full Text | Google Scholar

Xia, X., and Moog, C. (2003). Identifiability of nonlinear systems with application to HIV/AIDS models. IEEE Trans. Automat. Control 48, 330–336. doi: 10.1109/TAC.2002.808494

CrossRef Full Text | Google Scholar

Xue, H., Miao, H., and Wu, H. (2010). Sieve estimation of constant and time-varying coefficients in nonlinear ordinary differential equation models by considering both numerical error and measurement error. Ann. Statist. 38, 2351–2387. doi: 10.1214/09-AOS784

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar