Ornstein-Uhlenbeck model: maximum likelihood estimates from 1000 simulations of model Eq. (14), using the exact and the approximated transition density.

## Abstract

This chapter is concerned with estimation method for multidimensional and nonlinear dynamical models including stochastic differential equations containing random effects (random parameters). This type of model has proved useful for describing continuous random processes, for distinguishing intra- and interindividual variability as well as for accounting for uncertainty in the dynamic model itself. Pharmacokinetic/pharmacodynamic modeling often involves repeated measurements on a series of experimental units, and random effects are incorporated into the model to simulate the individual behavior in the entire population. Unfortunately, the estimation of this kind of models could involve some difficulties, because in most cases, the transition density of the diffusion process given the random effects is not available. In this work, we focus on the approximation of the transition density of a such process in a closed form in order to obtain parameter estimates in this kind of model, using the Fokker-Planck equation and the Risken approximation. In addition, the chapter discusses a simulation study using Markov Chain Monte Carlo simulation, to provide results of the proposed methodology and to illustrate an application of mixed effects models with SDEs in the epidemiology using the minimal model describing glucose-insulin kinetics.

### Keywords

- stochastic differential equations
- mixed effects model
- Fokker-Planck equation
- transition density
- maximum likelihood estimators
- genetic algorithm
- Markov Chain Monte Carlo simulation

## 1. Introduction

In pharmacokinetic/pharmacodynamic studies, repeated measurements are performed on a sample of individuals (units/subjects), and responses for all experimental subjects are assumed to be described by a common structural model. This model contains both fixed and random effects to distinguish between individuals in a population, leading to a mixed effects model in which fixed effects represent fixed parameters for all individuals in the population and random effects account for individual differences. Moreover, mixed effects modeling has been shown to be useful in pharmacokinetic/pharmacodynamic studies; particularly for modeling total variation in and between individuals, the models used in pharmacokinetic/pharmacodynamic analysis are often presented by a system of deterministic differential equations (ordinary, partial, or delay). However, the pharmacological processes in reality are always exposed to incomprehensible effects that are difficult to model, and ignoring these effects can affect the parameter estimation results and their interpretations. So, the introduction of stochastic components to deterministic models is an important tool of analysis [1] and is more appropriate to model the intra-individual variations rather than ODEs. In addition, the extension of ODEs to the SDEs makes it possible to explain the differences between the observations and the predictions by two types of noise: dynamic noise that enters through the dynamics of the system and that can result from its random fluctuations or from the shortages of model and the measurement error which is added in the case of an indirectly observed process, which may be due to a test error or to the existence of a disturbance and represent the uncorrelated part of the residual variability. In the theory, there are rich and developed resources for mixed effects models whether deterministic [2, 3, 4, 5, 6] or stochastic and linear or nonlinear models (see many applications of stochastic NLME models in biomedical [7, 8, 9] and in pharmacokinetic [10, 11, 12] fields).

Parameter estimation in mixed effects models with SDEs, known by stochastic differential mixed effects models, is not an obvious procedure except in some cases simpler [14] because it is often difficult to write the likelihood function in its closed form. In this context, we propose a review on estimation methods of SDME models in [13, 15] and an example case that treats a generalized linear mixed models in [16] and also an example to approximate the likelihood function of an NLME with the likelihood of a linear mixed effects model in [17]. Moreover, to strengthen knowledge on estimation methods of SDME, we refer to [18, 19] that propose an example of stochastic mixed effects model with random effects log-normally distributed with a constant diffusion term.

In general, it is difficult to obtain an explicit likelihood function because the transition density of the stochastic process is often unknown or that the integral in the likelihood given the random effects cannot be computed analytically, and although the size of the random effects increases, the complexity of the problem increases also rapidly. Therefore, this requires a significant need for approximate methods to compute the transition density in an approximate closed form and for effective numerical integration methods to compute or approximate the integral in the likelihood function, using, for example, the Laplacian and Gaussian quadrature approximation [3, 8, 20] or other approaches [21, 22]. In the literature, several solutions have been proposed to approximate the transition density and have shown their effectiveness despite certain limitations. For example, the transition density could be approximated by the solution of the partial differential equations of Kolmogorov [23]; by the derivation of a Hermite expansion of closed form at the transition density [24, 25, 26] (this method has been reviewed and applied for many known stochastic processes for one-dimensional [8] and multidimensional [20] time-homogeneous SDME model); or by simulating the process to Monte-Carlo-integrate the transition density [27, 28, 29]. These techniques are very useful and can solve the problem, but unfortunately, they involve intense calculations which make the problem always complicated.

In this work, we focus on two fundamental issues concerning the implementation of SDEs in NLME models. The first is how the transition density of an SDME model can be approximated when it is not known, and the second is about approximating methods of the likelihood function when the integral given the random effects has no analytic solution. Then, we propose an optimization algorithm to obtain maximum likelihood estimators in order to facilitate the estimation procedure for these models. Finally, the methodology is evaluated by simulation studies on the bidimensional Ornstein-Uhlenbeck model, and then it is implemented on the minimal model describing the glucose and insulin kinetics.

## 2. Theoretical tools

Consider an N-dimensional continuous and stochastic process

where

The joint density function of the vector

According to the model Eq. (1), the process Y is the same and follows the same form for each individual of the population, and the SDME model of type Itô expresses the dynamic of the individual i perturbed by the Brownian motion. So, the differences between the subjects are modeled on the one hand by the different realizations of the Brownian motion paths

The goal is to estimate the vector of fixed parameters

### 2.1 Maximum likelihood estimation in SDME model

The likelihood function of an SDME model is expressed as follows:

with

where

In this work, we propose an approximate method to estimate the transition density for a time-inhomogeneous NLME model with SDEs in a closed approximate form. So, we suggest to derive the transition density by solving the motion equations of the process using the Fokker-Planck equation, and we deal with the use of its solution given in [34]. Then, using the expression obtained, we get a closed form approximation of the likelihood function that we maximize to obtain the approximated MLEs

#### 2.1.1 Laplace approximation

For a multidimensional vector of random parameters, if the exact transition density or its closed form approximation is available, we can use the Laplace approximation method [16, 20, 35] to obtain an explicit expression of the likelihood function to maximize, despite the fact that the integral in Eq. (4) has no closed solution. So, for a q-dimensional random vector

where

and

### 2.2 Approximate transition density for multidimensional and nonlinear SDME model

Usually, a formal general solution of the stochastic differential equations as in Eq. (1) cannot be given, which makes the calculation of the transition density for this process more complicated. Moreover, this process has a lot of fluctuations that its exact position cannot be determined but can be known given a region by its probability density; with the FP equation, such a probability density can be determined. The FP equation is a differential equation for the distribution function describing a Brownian motion by which the probability density of the stochastic process can be calculated in a much simpler way by solving this equation. This motion equation is usually used for variables describing a macroscopic but small system, where the fluctuations are important as for some cases in physics, e.g., the position and the speed of the Brownian motion of a small particle. However, it can be also used for the larger system where, in spite of their small fluctuations, the stochastic description remains necessary when the deterministic equations may not be stable for this type of system.

So, under some assumptions, the transition density

where

In [34], H. Risken deals in this book with the derivation of the FP equation and its solution methods with some of its applications especially for problems of Brownian motion. So, here, we propose to characterize the transition density of the N-dimensional process

To test the effectiveness of this approach, we will guide our statistical methodology by simulation studies in order to examine the flexibility of its application to deduce its advantages and disadvantages. So, we substitute

#### 2.2.1 Approximated estimators

For a nonlinear SDME model with Gaussian random effects using Eqs. (6), (7), and (11), we obtain the following approximated likelihood function:

The MLEs of

### 2.3 Genetic algorithm

The genetic algorithm is a random search technique to look for an exact or approximated optimum points for optimization problems [36, 37, 38]. It is based on the concepts of natural genetic evolution which contains the following stages: the reproduction, the crossing, and the mutation of a constantly evolving population. It sets up the evolution of a random population of potential solutions of the N cardinal; then, the N simultaneous iterative trajectories interact with each other by following or imitating the biological evolution, for a convergence of certain elements of the population towards an optimal point of the fitness function.

The GA can search in multiple directions to explore all the search space by the possibility of jumping across them, so that the seeds spread uniformly over the whole search space. In this algorithm, we have a diversity of initial populations which gives the global optimum faster than other algorithms, where the initial value is very important and should be enough close to the global optimum. All of these features allows the GA to be regarded as a driving tool of evolution giving good results for optimization processes [37, 39]. In the literature, there were many works about the application of GA in optimizing problems specially for likelihood function [40, 41]. For the use of GA, we must first define some parameters of the algorithm: Population size N, EN, SR, CP, MP, fitness function, and convergence criteria. In the following we present the GA steps:

** 1**. Generate initial population

For

** 2**. Evaluate the Fitness function

** 3**. While (convergence criteria are not satisfied):

** Do**:

** 4**. Replacement step (by using SR and EN): At the SR rate, individuals with the worst results in step 2 of fitness function are replaced by new ones randomly generated, and a number EN of individuals is selected and accepted for the next step.

** 5**. Selection operator by using roulette wheel method, based on the fact that the more the individual has a good result of fitness function, the more likely he will be selected.

** 6**. Crossover operator by using CP and mutation operator by using MP: It is a mechanism of perturbation on the candidate individuals (parents) according to CP and MP to generate new groups of individuals, and we obtain a new

** Else**:

** 7**. Evolution stops; get GA output.

** 8**.

End For.

In the simulation studies paragraphs, the GA is implemented using the MATLAB software, where the function “ga,” to generate the genetic algorithm, requires inputs that are chosen according to the constraints of each example (see the help window in MATLAB). Moreover, the algorithm parameters are chosen according to what is early used in the literature [39],

## 3. Simulation studies for SDME model

### 3.1 The two-dimensional Ornstein-Uhlenbeck process

To apply the proposed methodology and evaluate its effectiveness, we consider the two-dimensional OU process that is very useful in pharmacokinetic/pharmacodynamic studies and in biology [43], physics, engineering, finance, and neuroscience applications [14, 44]. Indeed, the choice of this process is due to the fact that it is one of the few known multivariate SDME models with known transition density. For this reason, we choose the OU process to evaluate the methodology presented above, and we perform a comparison study between the results obtained using the proposed transition density in Eq. (11) and those obtained using the exact density. The model is defined as follows:

With

We rewrite the system in matrix notation under the Itô formula; we denote by

Here, the random effects

with mean vector

We assume that the matrices

In Figure 1, we illustrate in (a) the simulation of the OU process using the Euler scheme [45] with the following set of parameters: _{,} and a time step of

#### 3.1.1 Simulation results

In this simulation study, we generate 1000 artificial datasets of dimension

By plugging Eq. (16) in Eq. (4) and Eq. (18) in Eq. (6), we obtain a huge expression of the likelihood function but without a closed solution of integrals, so the exact estimators of

Table 1 shows that, for the given sample size, the results can be correctly identified using this estimated approach; even if some of parameters are overestimated or underestimated, the results remain acceptable because the results belong to the confidence interval. However, we believe that these results could be further proven by using other sample sizes and by adding alternative assumptions for the model that we did not consider in the methodology proposed in this chapter, which could further complicate the problem and be more time-consuming, as well as the present methodology suffers from some limitations. For the random parameters, the estimates can be provided using the optimization algorithm on Eq. (8) using the obtained estimate results of the parameters vector

True values | Mean and (Std) ( | ||||||
---|---|---|---|---|---|---|---|

2.8 | 2.5 | 1.8 | 2 | 3.10–3.25 | 2.48–2.56 | 1.72–1.65 | 2.11–2.37 |

(0.164)–(0.283) | (0.015)–(0.283) | (0.095)–(0.189) | (0.153)–(0.255) | ||||

0.8 | 1.5 | 0.3 | 0.5 | 1.06–1.13 | 1.57–1.64 | 0.28–0.37 | 0.56–0.45 |

(0.071)–(0.102) | (0.109)–(0.158) | (0.026)–(0.073) | (0.023)–(0.061) | ||||

45 | 100 | 100 | 25 | 44.75–52.43 | 100–112.75 | 102.35–89.64 | 24.72–31.02 |

(0.523)–(9.372) | (01.166)–(28.113) | (01.04)–(22.05) | (2.297)–(11.20) |

### 3.2 Stochastic minimal model

The minimal model describes the glucose-insulin kinetics and the dynamics of these processes illustrating the diabetes disease mechanisms. The diabetes is one of the most prevalent diseases in individuals and the nature and degree of its assignment changes from an individual to another and depends on certain individual characteristics, which implies that the concepts of stochastic modeling with random effects could be a good approach to modeling this disease. The Diabetes may be due to the insufficient insulin production (type 1 diabetes) or to the fact that the cells do not respond to the secreted insulin (type 2 diabetes) and the T2D patients tend to have substantially lower insulin sensitivity than healthy individuals, so that the T2D can be characterized by the level of insulin sensitivity for each individual. Therefore, to model the T2D, we observe how a person’s body responds to insulin in the process of transporting glucose to tissues by measuring his insulin sensitivity. So, we present in this section the estimation of the minimal model which represents a powerful model describing the glucose-insulin kinetics for an individual’s body in three differential equations, see the mathematical formulation of the model in [46, 47, 48, 49]. So, it is already clear that the model will contain both fixed and random effects, because the study of diabetes disease takes into account the response of each individual according to his own parameters and other common parameters that describe the process of glucose-insulin for the entire population. See the description of the glucose-insulin kinetics in Figure 3.

From the mentioned literature and Figure 3, the glucose-insulin disposal can be represented, with respect to time, by the following nonlinear stochastic differential equations, perturbed by the stochastic terms

where

The insulin sensitivity is defined by combining the two rates

Finally, the stochastic minimal model under the Itô sense re-parameterized by

where

Random effects are assumed to be independent with a multinormal joint density function, with the mean

So, here we deal with a time-inhomogeneous NLME model with SDE describing the glucose-insulin kinetics (see [52] for the implementation of SDE time-inhomogeneous model) and [53] where the maximum likelihood estimation for a time-inhomogeneous stochastic differential model of glucose dynamics was treated. The measurements in the model Eq. (20) are assumed directly observed without measurement errors as well as in the theoretical approach presented above.

So, we wish to estimate

with

We have no closed form solution to this integral, so exact estimators of

#### 3.2.1 Simulation results

We start this study with an application on artificial data generated on the intravenous glucose tolerance test principle (see [53] for a mathematical modeling of the test where glucose and insulin concentrations in plasma are subsequently sampled after an intravenous glucose injection). We generate 5000 sets of simulated artificial data of dimension

For each data from 5000 generated data sets, we estimate

True values | (M,m) | Mean and (Std) | ||||||
---|---|---|---|---|---|---|---|---|

0.074 | 0.10 | 0.0007 | 90 | (40,60) | 0.0737 (0.0014) | 0.100 (0.0013) | 0.00073 (2.02 × 10^{−4}) | 89.08 (0.031) |

(40,10) | 0.0784 (0.0059) | 0.209 (0.0295) | 0.00068 (2.54 × 10^{−4}) | 91.34 (0.712) | ||||

(10,20) | 0.0794 (0.0402) | 0.156 (0.046) | 0.00054 (0.0016) | 62.11 (1.13) | ||||

0.01 | 0.06 | 0.03 | 0.006 | (40,60) | 0.009 (0.0012) | 0.0616 (0.0011) | 0.0343 (0.0027) | 0.0061 (0.00010) |

(40,10) | 0.014 (0.0023) | 0.0708 (0.0060) | 0.0340 (0.0035) | 0.0073 (0.00030) | ||||

(10,20) | 0.015 (0.0031) | 0.0938 (1.0196) | 0.0480 (0.0108) | 0.0067 (0.0006) | ||||

0.000025 | 46 | 50 | 0.03 | (40,60) | 0.000021 (3 × 10^{−6}) | 44.98 (1.25) | 45.96 (2.11) | 0.0315 (0.0007) |

(40,10) | 0.000029 (4.4 × 10^{−6}) | 43.94 (1.75) | 45.16 (2.53) | 0.0349 (0.0036) | ||||

(10,20) | 0.000016 (0.8 × 10^{−5}) | 41.09 (2.82) | 44.64 (3.07) | 0.0178 (0.0051) | ||||

0.0002 | 95 | 320 | (40,60) | 0.00021 (1.2 × 10^{−6}) | 94.20 (1.12) | 321.2 (1.07) | ||

(40,10) | 0.00025 (1.6 × 10^{−6}) | 92.05 (2.25) | 318.6 (3.11) | |||||

(10,20) | 0.00037 (1.35 × 10^{−5}) | 122.51 (2.51) | 281.5 (4.88) |

In Table 2, we report the results obtained on large and small data by maximizing Eq. (22) using the GA. We notice that, for the case of the large data, the true values of the parameters are well identified with the exception of some in which, in all the simulations and samples, the true value does not belong to the estimated confidence interval, such as

Finally, from this simulation study where we have considered two SDME models, we can conclude that the parameters values of the models appear to be correctly identified using the proposed approach based on the Risken approximation to approximate the transition density of a stochastic process.

## 4. Conclusions

In this work, we have proposed a procedure to estimate the parameters of a mixed effects model containing stochastic differential equations, known by the SDME models, by proposing an approach to approximate its likelihood function to obtain the MLEs. This method has been evaluated by simulation studies on two SDME models in epidemiology: the two-dimensional Ornstein-Uhlenbeck process and stochastic minimal model. In fact, in models with SDEs instead of ODEs with random effects, the estimation of parameters is still not obvious even for one individual (one trajectory) because of the difficulties in deriving the transition densities, and difficulties become more interesting in using the population approach that treats the entire population simultaneously. The derivation of the exact density is not always possible for a stochastic and continuous process in an SDME model, so the search for an approximation of this density is an important step and requires an expensive calculation. This task is very interesting to give good results with good statistical properties of the estimators obtained by maximizing the likelihood function. In this work, the proposed estimation method has been applied to multidimensional and nonlinear SDME models with many random parameters normally distributed that can be extended to random parameters of any distribution. So, an approximation of the transition density

The classical inference of SDME models implies the problem of the numerical evaluation of the integral for the given random effects in the likelihood function, which becomes complicated especially when the model contains more than two random parameters. In the literature, several methods have been proposed and tested for the approximation of the integral (see references in the introduction) and the following examples: [8] which proposes the Gaussian quadrature method to solve the integrals for the case of an SDME models with a single random effect and [20] where the study was revised for a general case with several random parameters using the Laplace approximation to compute the integral in Eq. (4) and Eq. (6) numerically. For the mixed effects framework, see [3, 16, 55, 56]. In the case of using the Laplace method, as in this chapter, the calculation of the Hessian matrix can be done analytically when it is possible, as the examples in Section 2, or with the help of a symbolic calculus software or the automatic differentiation (AD) tools [57].

The results of simulation studies are satisfactory and can be obtained either by using moderate values for the number of experimental units M and of observations n taken for each experimental unit or by using a small sample size but with a number of measurements taken for each subject of 10 at least; this is relevant for applications where large sets of data are not available, such as biomedical applications where the mixed effect theory is widely applied.

The advantages of this approach, compared to those proposed in the literature for multidimensional SDME models with more than one random parameter [26], are that the computation of the approximate density is very easy and does not require a lot of time to calculate it or to program it in a software; the only task that can be time-consuming is in the optimization step to search for the optimum solution of the likelihood function, and also the proposed method is effective even with large data with a MATLAB program on a common PC (Intel Pentium IV 3.0 GHz with 512 MB of RAM). Nevertheless, the method suffers some limitations, for example, when the conditions to use Eq. (11) are not verified when, e.g., the inverse of the diffusion term does not exist and when, in certain cases, it is not obvious to derive the gradients and Hessians terms. Another limitation is that measurement errors are not considered in this work, and for a good stochastic version, it will be better to include noise on process increments and noise on observations that may be significant compared to system noise. These limitations may provide a perspective towards a more elaborate extension of the statistical study for SDME models, particularly in the field of epidemiology.

In conclusion, in this work, we proposed a method of parameter estimation for a mixed effects models with SDEs proposing an approximation method for the transition density in the case when it cannot be obtained in a closed form, with an approximation method of integral computation for the case of an SDME model with several random parameters. We have treated two examples to illustrate the effectiveness of this approach using computer tools. Indeed, we believe that this type of models is very interesting and provides a powerful and flexible modeling approach for repeated measurement studies such as biological and pharmacokinetic/pharmacodynamic and financial studies, as they combine the good characteristics of mixed effects and stochastic increments in intra-subject dynamics for a good modeling of a phenomenon.

## Acknowledgments

Z.T., H.E., and F.B conceived the presented idea and contributed to the design and implementation of the research. F.B., H.E., Z.T., and H.A contributed to the analysis of the results. Both F.B. and H.E authors developed the theory and performed the computations and the numerical simulations and verified the analytic methods and calculations.

## Abbreviations

stochastic differential equations

Markov Chain Monte Carlo

Fokker-Planck equation

maximum likelihood estimators

ordinary differential equations

nonlinear mixed effects

stochastic differential mixed effects

Ornstein-Uhlenbeck

genetic algorithm

elite number

selection rate

crossover probability

mutation probability

stochastic minimal model

type 2 diabetes