simode: R Package for statistical inference of ordinary differential equations using separable integral-matching

In this paper we describe simode: Separable Integral Matching for Ordinary Differential Equations. The statistical methodologies applied in the package focus on several minimization procedures of an integral-matching criterion function, taking advantage of the mathematical structure of the differential equations like separability of parameters from equations. Application of integral based methods to parameter estimation of ordinary differential equations was shown to yield more accurate and stable results comparing to derivative based ones. Linear features such as separability were shown to ease optimization and inference. We demonstrate the functionalities of the package using various systems of ordinary differential equations.


Background
This paper presents the simode R package [1] aimed for conducting statistical inference on systems of ordinary differential equations (ODEs). Systems of ODEs are commonly used for the mathematical modeling of the rate of change of dynamic processes such as in mathematical biology [2], biochemistry [3] and compartmental models in epidemiology [4], to mention a few areas. Inference of ODEs involves the 'standard' statistical problems such as studying the identifiability of a model, estimating model parameters, predicting future states of the system, testing hypotheses, and choosing the 'best' model. However, dynamic systems are typically very complex: nonlinear, high dimensional and only partly measured. Moreover, data may be sparse and noisy. Thus, statistical learning (inference, prediction) of dynamical systems is not a trivial task in practice. In particular, numerical application of standard estimators, like the maximum likelihood or the least squares, may be difficult or computationally costly. Therefore, special computational platforms that allow for performing statistical inference for ODEs were recently developed. We first briefly mention some relevant packages we are aware of and then point out the main focus of simode.
Existing software implementations that are most relevant to this work are the following. CollocInfer R package of [5] implements the profiling methodology of [6] and some extensions (there exist also a Matlab version). In the area of systems biology, [7] present Data2Dynamics, a modeling environment for Matlab that can be used for constructing dynamical models of biochemical reaction networks for large datasets and complex experimental conditions, and to perform efficient and reliable parameter estimation for model fitting. [8] developed the episode R package that implements adaptive integral-matching (AIM) algorithm for learning polynomial or rational ODEs with a sparse network structure. Other software libraries not directly related to our methodological framework are described or used in [9], [10], [11], and [12] which focus on stochastic modeling or on more specific domains. [3] uses PLAS (Power Law Analysis and Simulation; Copyright 1996-2012 by António Ferreira) a software suitable to analyze power-law differential equations. Also developed by António Ferreira is S-timator, a Python library dedicated for analyzing ODE-based models.
The R package simode is substantially different from all the above tools in a sense that will be now explained and made clear.

The focus of simode
The statistical methodologies applied in the package are based on recent publications that study theoretical and applied aspects of smoothing methods in the context of ordinary differential equations ( [13], [14], [15], [16], [17]). In that sense simode is closer in spirit to CollocInfer R package of [5], and episode R package of [8]. Unlike CollocInfer we do not consider penalized estimation which balances between data and model. Further, we focus on integral-matching criterion functions which were shown to be more robust than gradient based ones [13]. Using integral-matching criteria takes us closer in spirit to the episode R package of [8]. However, we focus on several minimization procedures of an integral-matching criterion function, taking advantage of the mathematical structure of the ODEs like separability of parameters from equations. Linear features such as separability were shown to ease optimization and inference ( [13], [14], [18], [15], [16], [19]). We demonstrate various functionalities of the package using different systems of ODEs. To be more specific, we demonstrate the ability of the package to implement a full estimation pipeline from point estimates to generating confidence intervals (using an example of S-system); deal with partially observed systems (using SIR example); user defined likelihood functions and system decoupling(using FitzHugh-Nagumo model); Monte-Carlo and multiple subjects expeirments (using Lotka-Volterra example) and models with an external input functions (using a seasonally forced Lotka-Volterra example).
The idea of separability of parameters and equations is now explained. Consider the following simple biochemical system taken from Chapter 2, Page 54 of [3], (1) In this system the production of x 2 depends on x 1 , x 3 , and x 4 which enter with different kinetic orders (power). Specifically, x 3 has a negative power which indicates an inhibiting effect since an increase in x 3 leads to reduced production of x 2 . The dynamics of the system for x 1 (0) = 2 and x 2 (0) = 0.1 is shown in Figure 1. This system is a special case of an S-system ( [3]) defined as Here, α j , β j are rate constants; g jk , h jk are kinetic orders that reflect the strength and directionality of the effect a variable has on a given influx or efflux. The above system is linear in α j , β j but nonlinear in g jk , h jk . In fact, one can view this system as a regression where the 'covariates' variables are x j (t), the solutions of the ODEs on the right hand side of the equations, while the 'response' variables are the derivatives x j (t) on the left hand side. Further, we can say that the system is linear in its rate constants but nonlinear in the kinetic orders (the powers). More generally, consider a system of ordinary differential equations given by where The simode R package is especially useful for handling ODE systems for which where θ = (θ N L , θ L ) , stands for the matrix transpose. Here θ N L , a vector of size p N L , stands for the 'nonlinear' parameters that can not be separated from the state variables x, while θ L , a vector of size p L , are the 'linear' parameters; note that p = p L + p N L . Setting  (2) is a special case of (3) with a vector field F as given in (4). In the special case of the biochemical example (1), fixing x 3 , x 4 , we have d = 2, so that the matrix g(x(t); θ N L ) is given by . Thus, the system can be written in the form where θ L = (α 1 , β 1 , α 2 , β 2 ) = (2, 2.4, 4, 2) , and θ N L = (g 11 , g 12 , h 11 , h 12 , g 21 , g 22 , h 21 , h 22 ) = (0, 1, 0.5, 0, 0.1, 0, 0, 1) . The vector field in the formulation above is separable in the linear parameter vector θ L , and therefore we refer to such systems as ODEs linear in the parameter θ L (as in a linear regression model). This linear property of the system turns out to be very useful for data fitting purposes where parameter estimation is required. We emphasize that in our view simode should NOT be considered as a competitor for the other tools and packages mentioned above but instead, a complementary tool. Indeed, real problems arising in the area of dynamic systems are complex and typically there is no one method that can handle all type of problems uniformly better than other methods. For instance, one may consider generating initial parameter estimates using simode and then derive final estimates using generalized profiling implemented in CollocInfer that enables additional flexibility in assuming that the ODEs description of the dynamics are only approximately correct.
The paper is organized as follows. In the next section we briefly present the statistical methodology implemented in simode. In Section 3 we describe in detail the use of simode for parameter estimation of ODEs. Section 4 deals with partial observed systems, while additional functionalities are demonstrated in Section 5. The last section includes a summary and some future directions.

Statistical methodology
Let x(t; θ, ξ), t ∈ [0, T ] be the solution of the initial value problem (3) given values of ξ and θ. We assume measurements of x are collected at discrete time points where the random variables ij are independent measurement errors (not necessarily Gaussian) with zero mean and finite variance. Consider the nonlinear least squares estimator of θ and ξ defined as a minimiser with respect to η and ζ of the least squares criterion function Typically, the ODEs (3) are nonlinear in x, and no analytic solution of the system exists, therefore numerical integration techniques are required in the estimation process. In fact, in the least squares criterion above the exact solution x will be approximated byx, a numerical solution (e.g., using Runge-Kutta) of the ODEs equation (3) for a given parameter and initial conditions. Thus, estimation methods such as nonlinear least squares or maximum likelihood require solving the system numerically for large set of potential parameters values, and then choosing an optimal parameter using some nonlinear optimisation technique. However, the combination of sparse and noisy data, nonlinear optimisation, and the need for numerical integration makes the parameter estimation a complex task (even for systems of low dimensions, e.g., [20]), and in many instances requires heavy computation. Therefore we adopt a 'smooth and match' approach for parameter estimation which includes two steps (i): bypassing numerical integration by using nonparametric smoothing of the data, and (ii): estimating the parameters by fitting the ODEs model to the estimated functions ( [21], [13], [14], [18], [15], [17], [22], [23], [24], [25], [26]; see also Chapter 8 of [27]). Then the resulting estimates are used as initial guess for optimization of the nonlinear least squares criterion function (7). In particular, in the first estimation stage which includes the two steps mentioned above, we consider integral-matching which is now described.

Integral matching
By integration, equation (3) yields the system of integral equations Here x(t) = x(t; θ, ξ) is the true solution of the ODE. Letx(t) stand for a nonparametric estimator (e.g., smoothing the data using splines or local polynomials) of the solution x of the ODEs equation (3) given observations (6). The criterion function of an integral-matching approach for a fully observed systems of ODEs takes the form where · denotes the standard Euclidean norm. The estimator of the parameter will be the minimiser of the criterion function (9), with respect to ζ and η. As its name suggests, integral-matching avoids the estimation of derivatives of the solution x as done in other smooth and match applications and hence is more stable ( [13]). While applying integral-matching leads to stable estimators, minimizing the criterion (9) is a complicated task in practice and good initial guess of parameter values is required for optimization. Therefore, the R package simode is designed to take advantage of separability of the ODE system, as demonstrated above in equation (5). This separability issue is now further explained.

Exploiting linear features of the ODEs
The simode package implements three separability scenarios corresponding to the following cases of equation (3): (a) ODEs linear in the parameters where F (x(t); θ) = g(x(t))θ; (b) ODEs semi-linear in the parameters where F (x(t); θ) = g(x(t); θ N L )θ L ; (c) ODEs nonlinear in the parameters where F (x(t); θ) has no separable form that can be exploited.
On top of the estimation stability the integral-matching criteria ensures, cases (a)-(b) above enable better optimization. Indeed, the above cases describe mathematical characteristics (separability) of the ODEs, and in what follows we use a smoothed version of separable nonlinear least squares ( [28]) as well as 'classical' nonlinear least squares. These two optimization methods can be applied in both cases (a) and (b). However, case (c) characterizes a model for which the only optimization method applicable is a nonlinear one since there is no separability of parameters.
Consider case (a) of ODEs linear in the parameters where F (x(t); θ) = g(x(t))θ. DenoteĜ Minimizing the integral criterion function (9) with respect to ζ and η results in the direct estimatorŝ where I d denotes the d × d identity matrix. Note that these estimators are well defined only if the inverse matrices in (10) and (11) exist. Necessary and sufficient conditions for √ n-consistency of the 'direct integral estimators' (10) and (11) are provided in [13]. Furthermore, the extensive simulation study in the aforementioned paper has demonstrated that using integrals as above instead of derivatives yields more accurate estimates. Indeed, it is well known (see e.g., [3] and [29]) that estimating derivatives from noisy and sparse data may be rather inaccurate. Additional application of the direct integral method to a variety of synthetic and real data was shown to yield accurate and stable results in [17] and [18]. Clearly, in this special case of ODEs linear in the parameter θ the complex task of nonlinear optimization reduces to the least squares solutions (10) and (11) which are easy to obtain and therefore a substantial computational improvement in optimization performance is achieved. Now consider case (b) above of ODEs semi-linear in the parameters for which in equation (3) the model can be written in the form F (x(t); θ) = g(x(t); θ N L )θ L . Then, for a given θ N L , minimizing the integral criterion function (9) yields least squares solutions similar to (10)-(11) which we denote byξ(θ N L ) andθ L (θ N L ) with the notation emphasizing the dependence of the linear solutions on the nonlinear parameters. Plugging backξ(θ N L ) and θ L (θ N L ) into the integral criterion function (9) results with where we have definedĜ(t; is minimized and a solutionθ N L is obtained, estimators for ξ and θ follow immediately and are given byξ(θ N L ) andθ L (θ N L ), respectively. This optimization procedure was considered in [15] and is a form of separable nonlinear least squares ( [28]). Note that the we apply nonlinear optimization only for estimating the nonlinear parameters θ N L , and hence the dimension of the optimization problem has been substantially reduced.
Finally, case (c) above requires nonlinear optimization for estimating ξ and θ and the dimension of the optimization problem can not be reduced.

Parameter estimation of ODEs using simode
In this section we demonstrate the main functionality of simode package (version 1.1.2), for parameter estimation of ODEs, via exploring the cases (a) ODEs linear in the parameters where F (x(t); θ) = g(x(t))θ; (c) ODEs nonlinear in the parameters where F (x(t); θ) has no separable form that can be exploited.
We continue with the biochemical example described in equation (1). Consider equation (1) where we drop the last two equations and set x 3 = 0.5, x 4 = 1 as constants. Further, assume that the zero kinetic parameters are known , so that the system is given by Thus, the linear parameters of the system are θ L = (α 1 , β 1 , α 2 , β 2 ) = (2, 2.4, 4, 2) , and the nonlinear parameters are θ N L = (g 12 , h 11 , g 21 , h 22 ) = (1, 0.5, 0.1, 1) . Here is the system given in equation (13) as it should be written to be used later by the package: The code above is the symbolic setup of the ODE system. The resulting objects of variables, parameters and initial conditions are Since we are working with symbolic objects we have created a function solve ode that uses the ode function of deSolve package ( [30]). The following code generates observations according to the statistical model defined in equation (6) where the distribution of the measurement error is Gaussian with standard deviation of 0.05. The resulting 'true' ODE solutions and the stochastic observations are presented in Figure (2).

R> obs <-list() R> for(i in 1:length(vars)) { + obs[[i]] <-x_det[,i] + rnorm(n,0,sigma) + } R> names(obs) <-vars
Now that we have setup the system of ODEs in a symbolic form and generated observations from the statistical model we can explore cases (a)-(b)-(c) described above: we estimate model parameters, plot model fits, and provide profile-likelihood confidence intervals.
The package uses integral-matching as a first stage and then (by default) executes nonlinear least squares optimization, namely minimizing equation (7), starting from the integral-matching estimates. The first estimation stage, i.e., the integral-matching, is based on smoothing the observations. We use by default the smooth.spline method of the stats package, with generalized cross validation (we also support kernel smoothing and performing integralmatching without smoothing the observations). Our implementation of the estimators (10)-(11) uses lsqlincon function of the pracma package, which also allows us to introduce constraints on the parameters.

Case (a): ODEs linear in the parameters
For simplicity we begin with assuming that the initial conditions x 1 (0), x 2 (0) are known to the user. Here we also assume that the kinetic parameters are all known, so that our goal is to estimate the vector θ L = (α 1 , β 1 , α 2 , β 2 ) . The code for doing so is:    (1).
The call to simode returns a simode object containing the parameters estimates obtained using integral-matching as well as those obtained using nonlinear least squares optimization starting from the integral-matching estimates. An implementation of the generic plot function for simode objects can be used to plot the fit obtained using these estimates, either the fit after integral-matching and nonlinear least squares optimization (the default), just the integral-matching based fit, or both; see Figure 3. In this case, we can also plot the fit against the true curves since we know the true values of the parameters that were used to generate the observations. The same plot function can also be used to show the estimates obtained as demonstrated in the sequel.
In the code above we have defined which parameters are linear by lin pars and which are not by nlin pars. Then we fixed the set of parameters and initial conditions that we do not want to estimate. However, note that it is not mandatory for the user to know which parameters are linear and which are not. For instance, here is the result of running the estimation without this knowledge: As can be seen, the simode function generates error messages that point out exactly the nonlinear parameters. This is due to the fact that the simode function considers ODEs linear in the parameters as its default, namely im method = "separable", unless im method = "non-separable" is defined. This added simode feature makes it very useful for handling ODEs with linear features in case the mathematical knowledge for characterizing them is lacking. Now we generate and plot confidence intervals for the parameters using profile likelihood, see Figure 4. In case nonlinear optimization for the point estimates was used, then the profiling is done using a Gaussian based likelihood with fixed sigma which we estimate in the background.

Case (b): ODEs semi-linear in the parameters
Now consider ODEs semi-linear in the parameters where separable nonlinear least squares might be used as in equation (12). Thus, the nonlinear parameters θ N L = (g 12 , h 11 , g 21 , h 22 ) = (1, 0.5, 0.1, 1) are not assumed to be known and their estimation is needed. Estimating nonlinear parameters R> plot(profile_lin, mfrow=c(2,2))  requires nonlinear optimization. The function simode uses the optim function, thus we need to provide initial guess for optimization. In our example, the true parameter values are known, so we generate random initial guess in the vicinity of the true nonlinear parameters. The code and estimation results are given below. We can plot the resulting integral-maching and nonlinear least squares parameter estimates and compare them visualy to their true values, see Figure  5. So far we have used the default optimization of simode function which executes separable least squares. However, we could ignore the fact that the ODE we consider has linear features in its parameters. In that case we execute classical nonlinear optimization of the integral-matching criterion function for all the parameters, θ L , and θ N L .

Case (c): ODEs nonlinear in the parameters
There are cases where the system of ODEs is not separable, meaning that there are only nonlinear parameters. For instance, consider our toy example where the linear parameters, namely the coefficients are known. Still we may want to use integral-matching since this way we bypass the need for solving numerically the ODEs, at least in the first stage of optimization. Doing so leads to minimization of (9) as it is. The resulting code is:

Initial conditions x(0) are unknown
In the examples above, for simplicity of presentation, we considered the initial conditions to be known. Here is an example of estimating the initial conditions using the separability property of the ODEs. We extend case (b) above by adding the names of the unknown x 0 variables to the list of parameters to estimate.

Partially observed systems of ODEs
In general, inference using integral-matching requires a fully observed system. However, in some cases, integral-matching can be applied to a partially observed system. For example, if it is possible to reconstruct the unobserved variables using estimates of the system parameters. We demonstrate this using an example of an ODE system describing the spread of seasonal influenza in multiple age-groups across mutiple seasons. A discrete-time version of the model and a two-stage estimation procedure similar to the one used in simode, was described in details in [16]. Here we present the model and show how to employ the simode package in order to estimate its parameters. The model is an SIR-type (Susceptible-Infected-Recovered) model. The epidemic in each age-group 1 ≤ a ≤ M and each season 1 ≤ y ≤ L can be described using two equations for the proportion of susceptible (S) and infected (I) in the population (the proportion of recovered is given by 1 − S − I): S a,y (t) = −S a,y (t)κ y M j=1 β a,j I j,y (t), I a,y (t) = S a,y (t)κ y M j=1 (β a,j I j,y (t)) − γI a,y (t).
The parameters of the model include the M × M transmission matrix β, the recovery rate γ and κ 2,...,L which signify the relative infectivity of the influenza virus strains circulating in seasons 2, ..., L compared to season 1 (κ 1 is fixed as 1). As shown in [16], taking into account separability characteristics of this model is advantageous.
The simode package includes an example dataset called sir example, containing pre-made structures for testing this example. In this dataset there are two age-groups and five influenza seasons, so in total there are 10 equations for S and 10 equations for I. R> sir_example$S0 S1_1 S2_1 S1_2 S2_2 S1_3 S2_3 S1_4 S2_4 S1_5 S2_5 0.56 0.57 0.49 0.45 0.56 0.32 0.56 0.47 0.47 0.41

R> data(sir_example) R> summary(sir_example)
The dataset contains noisy observations for the I variables created using the parameter and initial condition values given in the example according to model (14), where the measurement errors have a Gaussian distribution with σ = 0.001. There are no observations of the S variables in the example as the proportion of susceptible in the population is typically unknown. Time is given in weeks and include 18 weeks of observations. The values of the parameters β and γ are also given assuming a time unit of weeks.

Additional functionalities of the package
In this section we demonstrate some additional functionalities of the package. We do that using other systems of ODEs in order to further explore the package usability.

User defined likelihood function
Consider the case where the user has her own likelihood function to be used in the second stage of optimization, meaning after the integral-matching stage. The default optimization of the package implements in the second stage the nonlinear least squares loss function (7). This default estimation procedure suffices for the method to result in consistent estimators. However, one may prefer to implement a specific likelihood function, for example a Gaussian distribution with known or unknown variance. We demonstrate this option using as an example the FitzHugh-Nagumo spike potential equations where the ODE model is given by See [5] for further explanation of this model. The above system is linear in a, b but nonlinear in c. The following code sets the equations, parameters and the 'true' values for initial conditions and parameters: R> pars <-c('a','b','c') R> vars <-c('V','R') R> eq_V <-'c*(V-V^3/3+R)' R> eq_R <-'-(V-a+b*R)/c' R> equations <-c(eq_V,eq_R) R> names(equations) <-vars R> x0 <-c(-1,1) R> names(x0) <-vars R> theta <-c(0.2,0.2,3) R> names(theta) <-pars The following code generates observations from a Gaussian measurement error model: Here we implement a Gaussian distribution (negative-log) likelihood function, which will be passed in the call to simode using the calc nll argument. The function should receive as arguments the parameter values, time points, observations and output of the solutions of the ODEs (calculated using the estimated parameter values in the current iteration of the optimization). Additional arguments can be passed to calc nll by passing them in the call to simode using the ellipsis construct (in this case the parameter sigma is passed as well). Note that the user-defined likelihood will also be used in the call to profile and the calculation of confidence intervals based on the likelihood profiles. Here is the user defined (negative-log) Gaussian likelihood function, }))) + } Now we demonstrate the usage of the likelihood function defined above, where the variance is assumed to be known and therefore is given as a fixed parameter. In this example the nonlinear parameters are assumed to be known. The resulting model fits are presented in Figure 9. The parameter estimates are  Now the above example is explored but with unknown variance which requires estimation as well. The user likelihood function has to be slightly modified as follows.    [20] demonsrated that system decoupling combined with data smoothing may lead to better reconstruction of the underlying dynamic system, and to better estimation of parameters. We implemented this functionality within the simode package. We use the ODE system of the previous example (15) for demonstrating decoupling functionality. The only difference is that we now set decouple equations=T in the simode function. The resulting integral-matching estimated values are stored in the returned simode object (as usual). If there is a parameter shared across system equations, the simode function will use the mean value (across system equations) of these parameter estimates (parameter 'c' in this example

Monte Carlo simulations
Conducting Monte Carlo experiments for ODEs can be an intensive computational task. The simode function can be given multiple sets of observations of a system from Monte Carlo simulations and fit each of these sets separately, returning a list of simode objects with the parameter estimates obtained from each fit. We demonstrate this using as an example the predator-prey Lotka-Volterra model ( [2]) given by the equations (X=prey, Y =predator): The Lotka-Volterra system is linear in all of its parameters. The following code sets the equations, parameters and the 'true' values for intial conditions and parameters: R> pars <-c('alpha','beta','gamma','delta') R> vars <-c('X','Y') R> eq_X <-'alpha*X-beta*X*Y' R> eq_Y <-'delta*X*Y-gamma*Y' R> equations <-c(eq_X,eq_Y) R> names(equations) <-vars R> x0 <-c(0.9,0.9) R> names(x0) <-vars R> theta <-c(2/3,4/3,1,1) R> names(theta) <-pars Next, we generate ten Monte Carlo sets of observations from a Gaussian measurement error model. The returned list.simode object has its own implementation of plot(). Setting the parameter plot mean sd=T the function plots the mean and standard deviation of the fitted curves (not shown here), or parameter estimates as in Figure 10, obtained from the multiple fits.

Multiple Subjects
There are cases where it is reasonable to consider a model where some parameters are assumed to be the same for all experimental subjects while other parameters are specific to an individual subject; see [31] for an example of mixed-effects modeling. While this is possible to do by defining a large ODE model where each individual has his own equations within this model, here we present an easier way to handle such a scenario using the simode package. We consider the Lotka-Volterra system (two equations)  Figure 10: Mean and standard deviation of the parameter estimates from ten Monte-Carlo simulations of the Lotka-Volterra model (16).
with N individuals. We assume that all individuals share the same system parameter values, whereas each individual has its own initial values. Hence, we would like to use the information from all individuals for estimating the system parameters. We call simode with the parameter obs containing a list of length N (where each member of this list is another list containing the observations for this specific subject), and set the parameter obs sets to N . We set the control parameter obs sets fit in simode.control to 'separate x0' to indicate we want to fit the same parameter values to all subjects but allow different initial conditions. Running simode returns an object of class list.simode, where each simode object contains the estimates for one individual. The parameter estimates in this case will be the same in each of the objects while the initial conditions estimates may be different. In the example below we set N = 5. Note that one can use the decoupling option here as well by setting decouple equations=T. In this case, the decoupling is designed such that the parameters of each equation will be estimated using the relevant data from all individuals together.

Models with an external input function
The ODE equations given to simode can include any function of time (e.g., forcing functions) using the reserved symbol 't'. To demonstrate this we expand the predator-prey Lotka-Volterra model from the previous section to include seasonal forcing of the predation rate, using two additional parameters that control the amplitude ( ) and phase (ω) of the forcing: X (t) = αX(t) − β(1 + sin(2π(t/T + ω)))X(t)Y (t), Y (t) = δ(1 + sin(2π(t/T + ω)))X(t)Y (t) − γY (t).
The parameter T sets the periodic time scale and is assumed to be known. The following code sets the equations, parameters and the 'true' values for the initial conditions and parameters assuming T = 50: Setting as nonlinear we can now proceed with the fit. The package enables the user to modify the optimization methods as is possible in the usual R implementation of optim function. Here we use the simplex ('Nelder-Mead') method instead of the gradient-based 'BFGS' method (the default), which yields the model fits presented in Figure 11.
R> nlin_pars <-c('epsilon','omega') R> nlin_init <-c(0.3,0.3) R> names(nlin_init) <-nlin_pars R> pars_min <-c(0,0) R> names(pars_min) <-nlin_pars R> pars_max <-c(1,1) R> names(pars_max) <-nlin_pars R> lv_force1 <-simode( + equations=equations, pars=pars, fixed=c(x0), time=time, obs=obs, + nlin_pars=nlin_pars, start=nlin_init, lower=pars_min, upper=pars_max, + simode_ctrl=simode.control(im_optim_method='Nelder-Mead', + nls_optim_method='Nelder-Mead')) R> summary(lv_force1)$est  The above example is based on an explicit input function sin(t). However, one can think of cases where we would like to use a general input, not necessarily for which there is a closed form. In such a case one would incorporate the external input by adding it to the list of observations (the obs parameter) and referencing it in the equations using the same name: several real-life scenarios such as partially observed systems, multiple subjects, models with an external input function, and user defined likelihood function. The package enables both splines and kernel smoothing methods. In addition, we implemented automatic system decoupling which in many cases leads to more stable numerical optimization. Furthermore, it is not mandatory for the user to know which parameters are linear and which are not. This added simode feature makes it very useful for handling ODEs with linear features in case the mathematical knowledge for characterizing them is lacking. Confidence intervals for the ODEs parameters can be calculated as well using profile likelihood, hence a full estimation pipeline is implemented. Finally, simode supports parallel Monte Carlo simulations. All above simode package functionalities are demonstrated using a variety of ODEs models: a biochemical system (S-system), SIR-type (Susceptible-Infected-Recovered), FitzHugh-Nagumo spike potential equations, and Lotka-Volterra. All the numerical examples presented in the paper are implemented as demos in the simode package. Future abilities are now developed, among them the use of kernel smoothing with automatic bandwidth selection as in [17], and a computational efficient parameter estimation for high dimensional systems of ordinary differential equations.