This article illustrates how Mathematica can be employed to model stochastic processes via stochastic differential equations to compute trajectories and their statistical features. In addition, we discuss parameter estimation of the model via the maximum likelihood method with global optimization. We consider handling modeling error, system noise and measurement error, compare the stochastic and deterministic models and use the likelihood-ratio test to verify the possible improvement provided by stochastic modeling. The Mathematica code is simple and can be easily adapted to similar problems.

Introduction

Stochastic differential equations (SDEs) have received great attention in a large number of fields, including finance, physics, system biology, biochemical processes and pharmacokinetics. SDEs serve as a natural way of introducing uncertainty into a deterministic model represented by ordinary differential equations (ODEs). In contrast to the classical approach where uncertainty only exists in the measurements, SDEs can provide a more flexible framework to account for deviation in states and parameters that describe the underlying system. For an introduction to the theory and numerical solution of SDEs, see [1].

In this article, we illustrate how Mathematica deals with the problem of simulation and parameter estimation of such systems, represented here by the FitzHugh–Nagumo model describing excitable media. The first section sets up the model equations. The second section deals with the simulation of the stochastic system and computes the trajectories, their distribution at a specified time point and the change of the statistical data of these distributions in time, namely the trajectory of the mean value and standard deviation. The third section describes the parameter estimation of the model via the maximum likelihood method (ML). The last section analyzes the goodness of fit between deterministic and stochastic models using the likelihood ratio test (LRT).

Modeling

The FitzHugh–Nagumo model for excitable media is a nonlinear model describing the reciprocal dependencies of the voltage across an exon membrane and a recovery variable summarizing outward currents. The model was developed in [2]. The model is general and can also model excitable media, for example, heart tissue. The deterministic ODE model (white-box model) is described by the following system of ordinary differential equations:

(1)

(2)

with parameters , and and initial condition and .

White-box models are mainly constructed on the basis of knowledge of physics about the system. Solutions to ODEs are deterministic functions of time, and hence these models are built on the assumption that the future value of the state variables can be predicted exactly. An essential part of model validation is the analysis of the residual errors (the deviation between the true observations and the one-step predictions provided by the model). This validation method is based on the fact that a correct model leads to uncorrelated residuals. This is rarely obtainable for white-box models. Hence, in these situations, it is not possible to validate ODE models using standard statistical tools. However, by using a slightly more advanced type of equation, this problem can be solved by replacing ODEs with SDEs that incorporate the stochastic behavior of the system: modeling error, unknown disturbances, system noise and so on.

The stochastic SDE gray-box model can be considered as an extension of the ODE model by introducing system noise:

(3)

(4)

where is a Wiener process (also known as Brownian motion), a continuous-time random walk. The next section carries out the numerical simulation of the SDE model using the parameter settings , and .

Simulation

Equations (3) and (4) represent an Ito-stochastic process that can be simulated in Mathematica employing a stochastic Runge–Kutta method.

Single Realization

First, a single realization is simulated in the time interval .

Figure 1. The trajectories of the state variables (blue) and (brown) in the case of a single realization of the Ito process.

The values of and can be found at any time point. Here is an example.

Many Realizations

This computes the trajectories for 100 realizations.

Figure 2. The trajectories of the state variables (blue) and (brown) in the case of 100 realizations of the Ito process.

Slice distributions of the state variables can be computed at any time point. First, let us simulate the trajectories in a slightly different and faster way, now using 1000 realizations.

At the time , we can compute the mean and standard deviation for both state variables and , as well as their histograms.

Figure 3. The distribution of (left) and (right) in the case of 100 realizations of the Ito process, at time .

The mean value and the standard deviation of the trajectories along the simulation time can also be computed and visualized. (This may take a few minutes.)

Figure 4. , the mean value of (red) with its standard deviation, (blue).

Parameter Estimation

Parameter estimation is critical since it determines how well the model compares to the measurement data. The measurement process itself may also have serially uncorrelated errors due to the imperfect accuracy and precision of the measurement equipment.

Measurement Values

Write the measurement equation as

(5)

where . The voltage is assumed to be sampled between and at discrete time points , where , , with an additive measurement noise . To get , we consider a single realization and sample it at time points . Then we add white noise with .

Figure 5. A single realization of (brown) and the simulated measured values (black points).

Likelihood Function

As we have seen, the solutions to SDEs are stochastic processes that are described by probability distributions. This property allows for maximum likelihood estimation. Let the deviation of the measurements from the model be

(6)

where . Assuming that the density function of can be approximated reasonably well by Gaussian density, the likelihood function to be maximized is

(7)

For computation, we use its logarithm. Now we accept the model parameter and estimate and from the SDE model, employing the maximum likelihood method.

Here are the values of the .

This defines the logarithmic likelihood function.

Maximization of the Likelihood Function

The optimization of the likelihood function is not an easy task, since the objective function is often flat, with nondifferentiable terms and many local extrema. In addition, the model takes a long time to evaluate. Instead of using direct global optimization, first we compute the values of the objective function on a grid to use parallel computation to speed up the evaluation. We are looking for the parameters in the range . Let us create the grid points.

This computes the function values at the grid points in parallel.

Now apply interpolation for the grid point data and visualize the likelihood function.

Figure 6. The likelihood function of the parameter estimation problem.

Employing different global optimization methods, we compute the parameters. The third one generates a warning message that we suppress with Quiet.

There are many local maxima. Let us choose the global one.

Simulation with the Estimated Parameters

Here is a new parameter list.

This computes 500 realizations.

This visualizes the result and the measurement data.

Figure 7. The simulated process with the estimated parameters as in Figure 4.

Deterministic versus Stochastic Modeling

In the previous section, Parameter Estimation, we illustrated the technique of stochastic modeling. The measurement data was simulated with the correct model parameters , no assumed modeling error, existing system noise and measurement error . The model to be fitted had two free parameters, and , since the system noise was chosen for the model.

Now consider a different situation. Suppose that we have modeling error, since the measured values are simulated with ; however, in the model to be fitted, we use . In addition, let us increase the measurement error to . To compare the efficiencies of the deterministic and stochastic models, we should handle as a free parameter in the stochastic model to be fitted. So then we have three free parameters to be estimated: , and . Let us carry out the parameter estimation for different values of .

The results can be seen in Table 1; corresponds to the deterministic model. In order to demonstrate that the stochastic model can provide significant improvement compared with the deterministic model, the likelihood-ratio test can be applied [3]. The likelihood-ratio test has been used to compare two nested models. In our case, one of the models is the ODE model, with fewer parameters than the SDE model. The null hypothesis is that the two models are basically the same. The test statistic is

(8)

where the subscripts and stand for the deterministic and stochastic models, respectively.

The distribution of is , where is the difference in the number of parameters between the two models; in our case, . Here is the critical value for at a confidence level of 95%.

If the value of is less than the critical value, then the null hypothesis can be accepted. In our case, . Therefore, we reject the null hypothesis, which means the SDE model can be considered as a different model providing a significant improvement compared with the ODE model.

For changing parameters to simulate the measured values, you should modify the content of the list in the Simulation section and the value of in the Measurement Values subsection. To change the parameter of the model to be fitted, change the value of in the list in the Likelihood Function subsection.

In reality, the three parameters should be handled simultaneously during the optimization process.

Assuming there is no modeling error, using the SDE model, we can separate the measurement error from the system noise represented by the estimated of the fitted model; see [4].

Conclusion

The advantage of using stochastic modeling is that the actual states in the model are predicted from data. This allows the prediction to stay close to the data even when the parameters in the model are imprecise. It has been demonstrated that Mathematica can significantly support a user carrying out stochastic modeling. There are many built-in functions that help in the computation of the trajectories of stochastic differential equations and their statistics, global optimization for parameter estimation and likelihood-ratio test for model verification.

Acknowledgments

The author is grateful for the EU FP7 International Research Staff Exchange Scheme, which provided the financial support that partly facilitated his stay at the Department of Mechanical Engineering of the University of Canterbury, NZ, where this article was completed.

References

[1]	E. Allen, Modeling with Itô Stochastic Differential Equations, Dordrecht, Netherlands: Springer, 2007.
[2]	R. FitzHugh, “Impulses and Physiological States in Theoretical Models of Nerve Membrane,” Biophysics Journal, 1(6), 1961 pp. 445–466. doi:10.1016/S0006-3495(61)86902-6.
[3]	H. Madsen and P. Thyregod, Introduction to General and Generalized Linear Models, Boca Raton: CRC Press, 2010.
[4]	C. W. Tornøe, J. L. Jacobsen, O. Pedersen, T. Hansen and H. Madsen, “Grey-Box Modelling of Pharmacokinetic/Pharmacodynamic Systems,” Journal of Pharmacokinetics and Pharmacodynamics, 31(5), 2004 pp. 401–417. doi.10.1007/s10928-004-8323-8.
B. Paláncz, “Stochastic Simulation and Parameter Estimation of the FitzHugh–Nagumo Model,” The Mathematica Journal, 2016. dx.doi.org/doi:10.3888/tmj.18-6.

About the Author

Dr. Béla Paláncz is professor emeritus at the Department of Photogrammetry and Geoinformatics at the Faculty of Civil Engineering of the Budapest University of Economics and Technology, in Budapest, Hungary. His main interests are mathematical modeling and numeric-symbolic computations. He is the coauthor of several books and more than a hundred papers.

Dr. Béla Paláncz:
Budapest University of Economics and Technology
1111 Műegytem rkp. 3
Budapest, Hungary
palancz@epito.bme.hu