**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:

$$\begin{cases} dY^{(1)i}(t) = -\left(\beta\_{11}b\_{11}^{i}\left(Y^{(1)i}(t) - a\_{1}\right) + \beta\_{12}b\_{12}^{i}\left(Y^{(2)i}(t) - a\_{2}\right)\right)dt + \Sigma\_{1i}d\mathcal{W}^{(1)i}(t), \quad Y\_{0}^{i} = y\_{0}^{(1)i}, i = 1, \ldots, M\\ dY^{(2)i}(t) = -\left(\beta\_{21}b\_{21}^{i}\left(Y^{(1)i}(t) - a\_{1}\right) + \beta\_{22}b\_{22}^{i}\left(Y^{(2)i}(t) - a\_{2}\right)\right)dt + \Sigma\_{22}d\mathcal{W}^{(2)i}(t), \quad Y\_{0}^{i} = y\_{0}^{(2)i}, i = 1, \ldots, M \end{cases} \tag{14}$$

$$\begin{split} & \text{With } Y^i(t) = \begin{pmatrix} Y^{(1)i}(t) \\ Y^{(2)i}(t) \end{pmatrix}; \beta = \begin{pmatrix} \beta\_{11} & \beta\_{12} \\ \beta\_{21} & \beta\_{22} \end{pmatrix}; a = \begin{pmatrix} a\_1 \\ a\_2 \end{pmatrix}; \Sigma = \begin{pmatrix} \Sigma\_{11} & 0 \\ 0 & \Sigma\_{22} \end{pmatrix}; \\ & W^i(t) = \begin{pmatrix} W^{(1)i}(t) \\ W^{(2)i}(t) \end{pmatrix}; Y^i(0) = \begin{pmatrix} Y^{(1)i}(0) \\ Y^{(1)i}(0) \end{pmatrix} \text{ and } b^i = \begin{pmatrix} b^i\_{11} & b^i\_{12} \\ b^i\_{21} & b^i\_{22} \end{pmatrix} \\ & \Gamma \begin{pmatrix} r\_{ll'} \cdot r\_{ll'}^{-1} \end{pmatrix}, \quad l, l' = \mathbf{1}, 2; \ i = \mathbf{1}, \ldots, M. \end{split}$$

We rewrite the system in matrix notation under the Itô formula; we denote by ð Þ� the elementwise multiplication:

$$d\boldsymbol{Y}^{i}(t) = \boldsymbol{\beta} \cdot \boldsymbol{b}^{i} (\boldsymbol{a} - \boldsymbol{Y}^{i}(t))dt + \Sigma d\boldsymbol{W}^{i}(t), \ \boldsymbol{Y}^{i}\_{0} = \boldsymbol{y}^{i}\_{0}, \ \boldsymbol{i} = \mathbf{1}, \ldots, \boldsymbol{M}. \tag{15}$$

Here, the random effects *bi* are a matrix and not a vector in order to have a uniform dimension writing of the Eq. (15) and are assumed mutually independent and independent of *Y<sup>i</sup>* <sup>0</sup> and *W<sup>i</sup>* . The fixed parameter vector is *θ* ¼ *β*11, *β*<sup>12</sup> ð , *β*21, *β*22, *α*1, *α*2, Σ11, Σ22Þ, and the population parameter vector is Ψ ¼ ð Þ *r*11, *r*12, *r*21, *r*<sup>22</sup> . The exact transition density of model Eq. (15) for a given realization of the random effects is a bivariate normal:

$$P\_Y\left(\mathbf{Y}\_{t\_j}^i, \Delta\_j^i | \mathbf{Y}\_{t\_{j-1}}^i, b^i, \theta\right) = \left(2\pi\right)^{-1} |\boldsymbol{\varepsilon}|^{\frac{-1}{2}} \exp\left(\frac{-\left(\mathbf{Y}\_{t\_j}^i - \mu\right)^{\prime} \boldsymbol{\varepsilon}^{-1} \left(\mathbf{Y}\_{t\_j}^i - \mu\right)}{2}\right),\tag{16}$$

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

$$\begin{aligned} \text{with mean vector } \mu &= a + \left( Y\_{t\_{j-1}}^{i} - a \right) \exp \left( - \left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right) \Delta\_{j}^{i} \right) \text{ and covariance matrix }\\ \boldsymbol{\xi} &= \boldsymbol{\tau} - \left( \exp \left( - \left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right) \Delta\_{j}^{i} \right) \quad \boldsymbol{\tau} \cdot \exp \left( - \left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right) \Delta\_{j}^{i} \right) \right), \text{ where} \\\\ \boldsymbol{\tau} &= \frac{1}{2 \text{tr}\left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right) | \boldsymbol{\beta} \boldsymbol{\beta}^{i} |} \left( \left| \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right| \boldsymbol{\Sigma} \boldsymbol{\Sigma}^{\prime} + \left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} - \text{tr}\left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right) \boldsymbol{I} \right) \boldsymbol{\Sigma} \boldsymbol{\Sigma}^{\prime} \left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} - \text{tr}\left( \boldsymbol{\beta} \boldsymbol{\beta}^{i} \right) \boldsymbol{I} \right)' \right). \end{aligned} \tag{17}$$

We assume that the matrices *β:b<sup>i</sup>* and Σ have full rank and the real parts of the eigenvalues of *β:bi* are positive definite in order that a stationary solution to Eq. (15) exists. Under these assumptions, we derive from Eqs. (14) and (11) the following approximated transition density of Y:

$$P\_{Y}^{(d)}\left(Y\_{t\_{l}}^{i},\Delta\_{j}^{i}|Y\_{t\_{l-1}}^{i},b^{i},\theta\right) = \left(2\sqrt{\Pi\Delta\_{j}}\right)^{-2}(\Sigma\_{11}\Sigma\_{22})^{-\frac{1}{2}} \ast \exp\left(-\frac{1}{4\Delta\_{j}}\left(\Sigma\_{11}^{-1}\right)\left[Y\_{t\_{l}}^{(1)\flat}-Y\_{t\_{l-1}}^{(1)\flat}+\left(\beta\_{11}b\_{11}^{i}\left(Y\_{t\_{l-1}}^{(1)\flat}-a\_{1}\right)\right)^{2}\right]\right) \tag{10.10}$$

$$\left(1+\beta\_{12}b\_{12}^{i}\left(Y\_{t\_{l-1}}^{(2)\flat}-a\_{2}\right)\right)\Delta\_{j}^{i}\right)^{2} + \left(\Sigma\_{22}^{-1}\right)\left[Y\_{t\_{l}}^{(2)\flat}-Y\_{t\_{l-1}}^{(2)\flat}+\left(\beta\_{21}b\_{21}^{i}\left(Y\_{t\_{l-1}}^{(1)\flat}-a\_{1}\right)+\beta\_{22}b\_{22}^{i}\left(Y\_{t\_{l-1}}^{(2)\flat}-a\_{2}\right)\right)\Delta\_{j}^{i}\right]^{2}\right). \tag{10.10}$$

In **Figure 1**, we illustrate in (a) the simulation of the OU process using the Euler scheme [45] with the following set of parameters:

#### **Figure 1.**

*A sample path of the OU process in the third graph of (a) for the given parameters set with the initial condition:* Y0 ¼ *(3,3) and time interval [3,10]; and the transition density for a transition from* Yj *to* Yjþ<sup>1</sup> *in (b).*

In the simulation studies paragraphs, the GA is implemented using the MATLAB

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.

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], *EN* ¼ 4, *MP* ¼ 0*:*2, *CP* ¼ 0*:*8 *and SR* ¼ 1*=*3, and the search spaces are around the confidence interval of the minimal model

parameters (see [42] and references therein).

*Numerical Modeling and Computer Simulation*

**3. Simulation studies for SDME model**

The model is defined as follows:

<sup>11</sup> *Y*ð Þ<sup>1</sup> *<sup>i</sup>*

<sup>21</sup> *Y*ð Þ<sup>1</sup> *<sup>i</sup>*

*Y*ð Þ<sup>1</sup> *<sup>i</sup>* ð Þ*t Y*ð Þ<sup>2</sup> *<sup>i</sup>* ð Þ*t*

!

; *Y<sup>i</sup>*

ðÞ� *t α*<sup>1</sup> � �

ðÞ�*t α*<sup>1</sup> � � <sup>þ</sup> *<sup>β</sup>*12*b<sup>i</sup>*

<sup>þ</sup> *<sup>β</sup>*22*b<sup>i</sup>*

� � � �

� � � �

ð Þ¼ <sup>0</sup> *<sup>Y</sup>*ð Þ<sup>1</sup> *<sup>i</sup>*

ð Þ*<sup>t</sup>* � �*dt* <sup>þ</sup> <sup>Σ</sup>*dW<sup>i</sup>*

<sup>0</sup> ¼ 1, 2; *i* ¼ 1, … , *M*.

<sup>0</sup> and *W<sup>i</sup>*

realization of the random effects is a bivariate normal:

<sup>¼</sup> ð Þ <sup>2</sup>*<sup>π</sup>* �<sup>1</sup>

<sup>12</sup> *Y*ð Þ<sup>2</sup> *<sup>i</sup>*

<sup>22</sup> *Y*ð Þ<sup>2</sup> *<sup>i</sup>*

; *<sup>β</sup>* <sup>¼</sup> *<sup>β</sup>*<sup>11</sup> *<sup>β</sup>*<sup>12</sup> *β*<sup>21</sup> *β*<sup>22</sup> � �

!

*Y*ð Þ<sup>1</sup> *<sup>i</sup>* ð Þ 0

ð Þ 0

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

Here, the random effects *bi* are a matrix and not a vector in order to have a uniform dimension writing of the Eq. (15) and are assumed mutually independent

*β*21, *β*22, *α*1, *α*2, Σ11, Σ22Þ, and the population parameter vector is Ψ ¼ ð Þ *r*11, *r*12, *r*21, *r*<sup>22</sup> . The exact transition density of model Eq. (15) for a given

> j j *ς* �1 <sup>2</sup> *exp*

ð Þ*<sup>t</sup>* , *<sup>Y</sup><sup>i</sup>*

� *<sup>Y</sup><sup>i</sup>*

0

B@

*tj* � *μ* � �<sup>0</sup>

*ς*�<sup>1</sup> *Y<sup>i</sup>*

2

*tj* � *μ* � �

1

CA, (16)

<sup>0</sup> <sup>¼</sup> *<sup>y</sup><sup>i</sup>*

. The fixed parameter vector is *θ* ¼ *β*11, *β*<sup>12</sup> ð ,

ðÞ�*t α*<sup>2</sup>

ð Þ�*t α*<sup>2</sup>

*dt* <sup>þ</sup> <sup>Σ</sup>11*dW*ð Þ<sup>1</sup> *<sup>i</sup>*

*dt* <sup>þ</sup> <sup>Σ</sup>22*dW*ð Þ<sup>2</sup> *<sup>i</sup>*

<sup>11</sup> *bi* 12

!

*bi* <sup>21</sup> *<sup>b</sup><sup>i</sup>* 22

; *<sup>α</sup>* <sup>¼</sup> *<sup>α</sup>*<sup>1</sup> *α*2 � �

and *<sup>b</sup><sup>i</sup>* <sup>¼</sup> *<sup>b</sup><sup>i</sup>*

ð Þ*<sup>t</sup>* , *<sup>Y</sup><sup>i</sup>* <sup>0</sup> ¼ *y* ð Þ1 *i*

ð Þ*<sup>t</sup>* , *<sup>Y</sup><sup>i</sup>* <sup>0</sup> ¼ *y* ð Þ2 *i*

; <sup>Σ</sup> <sup>¼</sup> <sup>Σ</sup><sup>11</sup> <sup>0</sup>

0 Σ<sup>22</sup> � �

where *b<sup>i</sup>*

<sup>0</sup>, *i* ¼ 1, … , *M:* (15)

<sup>0</sup> , *i* ¼ 1, … , *M*

<sup>0</sup> , *i* ¼ 1, … , *M* (14)

;

*ll*<sup>0</sup> *i:i:d* �

ðÞ¼� *<sup>t</sup> <sup>β</sup>*11*b<sup>i</sup>*

ðÞ¼� *<sup>t</sup> <sup>β</sup>*21*b<sup>i</sup>*

ðÞ¼ *t*

*W*ð Þ<sup>1</sup> *<sup>i</sup>* ð Þ*t*

!

*W*ð Þ<sup>2</sup> *<sup>i</sup>* ð Þ*t*

, *l*, *l*

*dY<sup>i</sup>*

and independent of *Y<sup>i</sup>*

*PY Y<sup>i</sup> tj* , Δ*<sup>i</sup> j* j*Yi tj*�<sup>1</sup> , *bi* , *θ*

**98**

ð Þ� the elementwise multiplication:

� �

ðÞ¼ *<sup>t</sup> <sup>β</sup>* � *bi <sup>α</sup>* � *<sup>Y</sup><sup>i</sup>*

With *Y<sup>i</sup>*

*dY*ð Þ<sup>1</sup> *<sup>i</sup>*

*dY*ð Þ<sup>2</sup> *<sup>i</sup>*

*W<sup>i</sup>* ðÞ¼ *t*

Γ *rll*0,*r*�<sup>1</sup> *ll*0 � �

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

ð ¼ *β*<sup>11</sup> ¼ 2*:*8, *β*<sup>12</sup> ¼ 2*:*5, *β*<sup>21</sup> ¼ 1*:*8, *β*<sup>22</sup> ¼ 2, *α*<sup>1</sup> ¼ 0*:*8, *α*<sup>2</sup> 1*:*5, Σ<sup>11</sup> ¼ 0*:*3, Σ<sup>22</sup> ¼ 0*:*5,*r*<sup>11</sup> ¼ 45,*r*<sup>12</sup> ¼ 100,*r*<sup>21</sup> ¼ 100, *r*<sup>22</sup> ¼ 125Þ, and a time step of Δ*<sup>j</sup>* ¼ 0*:*001 and represent in (b) the graph of the two transition densities given by Eqs. (16) and (18) for *Yj* to *Yj*þ<sup>1</sup> using the same set of parameters and time step.

### *3.1.1 Simulation results*

In this simulation study, we generate 1000 artificial datasets of dimension 2 � ð Þ� *n* þ 1 *M* from Eq. (15), where M is the number of subjects and n presents the number of repetitions of the experiment for each subject; then we estimate the parameters using the proposed approximated method, and we obtain 1000 sets of parameter estimates. The observations are obtained by linear interpolation from simulated trajectories using the Euler-Maruyama scheme with step size equal to 10�<sup>3</sup> .

distribution of the most approximated estimators seems to be reasonably close to a

*Empirical distribution of estimates obtained using the exact and approximated transition density.*

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

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–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

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 *σ*1*dw*1ð Þ*t* , *σ*2*dw*2ð Þ*t* , and *σ*3*dw*3ð Þ*t* :

*dt* <sup>þ</sup> *<sup>σ</sup>*1*dw*1ð Þ*<sup>t</sup>* , *<sup>G</sup>*ð Þ¼ <sup>0</sup> *<sup>G</sup>*<sup>0</sup>

where *G t*ð Þ and *I t*ð Þ are, respectively, the concentration of glucose and insulin at

concentration before the glucose injection, this injection will cause a disturbance of the concentrations according to the mechanism described in these equations, and these values are assumed known for each individual. And, *G*<sup>0</sup> and *I*<sup>0</sup> are the theoretical measure of the concentrations at glucose injection moment at the begin-

(19)

 *dt* <sup>þ</sup> *<sup>σ</sup>*2*dw*2ð Þ*<sup>t</sup>* , *<sup>X</sup>*ð Þ¼ <sup>0</sup> <sup>0</sup> *dI t*ðÞ¼ �½ � *nIt* ð Þþ ðÞ� *Ib γ*ð Þ *G t*ð Þ� *h t dt* þ *σ*3*dw*3ð Þ*t* , *I*ð Þ¼ 0 *I*0,

time t in the blood. *Gb* and *Ib* indicate the basal level of glucose and insulin

normal distribution:

**Figure 2.**

**3.2 Stochastic minimal model**

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

of the glucose-insulin kinetics in **Figure 3**.

*dG t*ðÞ¼ � *<sup>p</sup>*<sup>1</sup> <sup>þ</sup> *X t*ð Þ *G t*ðÞþ *<sup>p</sup>*1*Gb*

*dX t*ðÞ¼ �*p*2*X t*ðÞþ *p*3ð Þ *I t*ðÞ� *Ib*

ning of the experiment.

**101**

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 *θ* and Ψ are unavailable. Therefore, in both cases, either using the exact or the approximated transition density, the Hessian matrix in Laplace approximation can be obtained analytically after a tedious calculation; then we apply the GA to obtain parameter estimates. But we cannot ignore the time consumed by this algorithm to get the results because of the long and complex expressions.

**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 Ψ. Moreover, we conduct this simulation study using **Figure 2,** which shows that the empirical


#### **Table 1.**

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

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

**Figure 2.**

ð ¼ *β*<sup>11</sup> ¼ 2*:*8, *β*<sup>12</sup> ¼ 2*:*5, *β*<sup>21</sup> ¼ 1*:*8, *β*<sup>22</sup> ¼ 2, *α*<sup>1</sup> ¼ 0*:*8, *α*<sup>2</sup> 1*:*5, Σ<sup>11</sup> ¼ 0*:*3, Σ<sup>22</sup> ¼ 0*:*5,*r*<sup>11</sup> ¼ 45,*r*<sup>12</sup> ¼ 100,*r*<sup>21</sup> ¼ 100, *r*<sup>22</sup> ¼ 125Þ, and a time step of Δ*<sup>j</sup>* ¼ 0*:*001 and represent in (b) the graph of the two transition densities given by Eqs. (16) and (18) for *Yj* to

In this simulation study, we generate 1000 artificial datasets of dimension 2 � ð Þ� *n* þ 1 *M* from Eq. (15), where M is the number of subjects and n presents the number of repetitions of the experiment for each subject; then we estimate the parameters using the proposed approximated method, and we obtain 1000 sets of parameter estimates. The observations are obtained by linear interpolation from simulated trajectories using the Euler-Maruyama scheme with step size

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 *θ* and Ψ are unavailable. Therefore, in both cases, either using the exact or the approximated transition density, the Hessian matrix in Laplace approximation can be obtained analytically after a tedious calculation; then we apply the GA to obtain parameter estimates. But we cannot ignore the time consumed by this algorithm to get the results because of the long and complex

**Table 1** shows that, for the given sample size, the results can be correctly

*β***<sup>11</sup>** *β***<sup>12</sup>** *β***<sup>21</sup>** *β***<sup>22</sup>** ^*β***<sup>11</sup>** ^*β***<sup>12</sup>** ^*β***<sup>21</sup>** ^*β***<sup>22</sup>** 2.8 2.5 1.8 2 3.10–3.25 2.48–2.56 1.72–1.65 2.11–2.37

*α*<sup>1</sup> *α*<sup>2</sup> Σ<sup>11</sup> Σ<sup>22</sup> *α*^<sup>1</sup> *α*^<sup>2</sup> Σ^<sup>11</sup> Σ^<sup>22</sup> 0.8 1.5 0.3 0.5 1.06–1.13 1.57–1.64 0.28–0.37 0.56–0.45

*r*<sup>11</sup> *r*<sup>12</sup> *r*<sup>21</sup> *r*<sup>22</sup> ^*r*<sup>11</sup> ^*r*<sup>12</sup> ^*r*<sup>21</sup> ^*r*<sup>22</sup> 45 100 100 25 44.75–52.43 100–112.75 102.35–89.64 24.72–31.02

*Ornstein-Uhlenbeck model: maximum likelihood estimates from 1000 simulations of model Eq. (14), using the*

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 Ψ. Moreover, we conduct this simulation study using **Figure 2,** which shows that the empirical

**Mean and (Std) (M** ¼ **40, n** ¼ **10)**

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

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

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

identified using this estimated approach; even if some of parameters are

*Yj*þ<sup>1</sup> using the same set of parameters and time step.

*Numerical Modeling and Computer Simulation*

*3.1.1 Simulation results*

equal to 10�<sup>3</sup>

expressions.

**True values**

**Table 1.**

**100**

*exact and the approximated transition density.*

.

*Empirical distribution of estimates obtained using the exact and approximated transition density.*

distribution of the most approximated estimators seems to be reasonably close to a normal distribution:

#### **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–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 *σ*1*dw*1ð Þ*t* , *σ*2*dw*2ð Þ*t* , and *σ*3*dw*3ð Þ*t* :

$$\begin{aligned} dG(t) &= \left[ -\left( p\_1 + X(t) \right) G(t) + p\_1 G\_b \right] dt + \sigma\_1 dw\_1(t), & G(0) &= G\_0 \\ dX(t) &= \left[ -p\_2 X(t) + p\_3 (I(t) - I\_b) \right] dt + \sigma\_2 dw\_2(t), & X(0) &= 0 \\ dI(t) &= \left[ -n(I(t) - I\_b) + \gamma (G(t) - h)t \right] dt + \sigma\_3 dw\_3(t), & I(0) &= I\_0, \end{aligned} \tag{19}$$

where *G t*ð Þ and *I t*ð Þ are, respectively, the concentration of glucose and insulin at time t in the blood. *Gb* and *Ib* indicate the basal level of glucose and insulin concentration before the glucose injection, this injection will cause a disturbance of the concentrations according to the mechanism described in these equations, and these values are assumed known for each individual. And, *G*<sup>0</sup> and *I*<sup>0</sup> are the theoretical measure of the concentrations at glucose injection moment at the beginning of the experiment.

process for the entire population (all individuals). So, we have the following ran-

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

*i*

<sup>0</sup> � *N μI*<sup>0</sup> , *σI*<sup>0</sup>

Random effects are assumed to be independent with a multinormal joint density

So, here we deal with a time-inhomogeneous NLME model with SDE describing

model Eq. (20). By using the approximated transition density (Eq. (11)), we get the

� *<sup>σ</sup>*�<sup>1</sup> *<sup>I</sup>*<sup>0</sup> *I i* <sup>0</sup> � *μ<sup>I</sup>*<sup>0</sup> � �<sup>2</sup>

*l*

2 to obtain a closed form approximation to the log-likelihood function *log L*ð Þ *<sup>a</sup>* ð Þ *<sup>θ</sup>*, <sup>Ψ</sup> � � for the model Eq. (20); then by applying the GA, we get the

� *<sup>μ</sup><sup>l</sup> <sup>Y</sup><sup>i</sup>*

We have no closed form solution to this integral, so exact estimators of *θ* and Ψ are unavailable. So, we use the Laplace approximation method described in Section

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 ð Þ� *ni* þ 1 *M* from Eq. (20) using the Euler-Maruyama scheme [45] with a step size of 10�<sup>3</sup> and a set of true parameters that are chosen according to [53, 54] representing the normal range of parameters values to simulate healthy subjects (without diabetes), with M being the number of units and *ni* being the number of observations or repeated measurements collected on each unit i.

For each data from 5000 generated data sets, we estimate ð Þ *θ*, Ψ by applying the proposed method. We first assume that the number of repeated measurements

� �<sup>2</sup> � �Þ*dSG*

2 ffiffiffiffiffiffiffiffi ΠΔ*<sup>j</sup>* � � <sup>p</sup> �<sup>3</sup> � �½ � *<sup>σ</sup>*1*σ*2*σ*<sup>3</sup> �<sup>1</sup>

, *μI*<sup>0</sup> , *μG*<sup>0</sup>

the glucose-insulin kinetics (see [52] for the implementation of SDE timeinhomogeneous 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 ð Þ *<sup>θ</sup>*, <sup>Ψ</sup> given the observations *<sup>y</sup>* <sup>¼</sup> *<sup>y</sup>*1, … , *<sup>y</sup><sup>M</sup>*

following approximated likelihood function for model Eq. (20):

�1 2 Y*ni j*¼1

� *<sup>σ</sup>*�<sup>1</sup> *SI <sup>S</sup><sup>i</sup> <sup>I</sup>* � *μSI* � �<sup>2</sup>

ð Þ *σSGσSIσ<sup>I</sup>*0*σ<sup>G</sup>*<sup>0</sup>

� �<sup>1</sup> 4Δ*<sup>j</sup>* � 1 *σ*1 *Ai* 1*j* � �<sup>2</sup> þ 1 *σ*2 *Ai* 2*j* � �<sup>2</sup> þ 1 *σ*3 *Ai* 3*j* � �<sup>2</sup>�i

� �, *G<sup>i</sup>*

� � and the covariance matrix *<sup>ϕ</sup>* <sup>¼</sup>

2Þ

� *<sup>σ</sup>*�<sup>1</sup> *<sup>G</sup>*<sup>0</sup> *<sup>G</sup><sup>i</sup>* <sup>0</sup> � *μ<sup>G</sup>*<sup>0</sup>

*<sup>j</sup>*�1, *tj*�1, *<sup>θ</sup>*, *<sup>b</sup><sup>i</sup>* � �Δ*<sup>i</sup>*

h i*:* (23)

*j*

, *μI*<sup>0</sup> , *μG*<sup>0</sup> , *σSG* , *σSI* , *σI*<sup>0</sup> , *σG*<sup>0</sup> � � and *<sup>θ</sup>* <sup>¼</sup>

<sup>0</sup> � *N μG*<sup>0</sup> , *σG*<sup>0</sup>

� �*:* (21)

� � from

*i dSI i dI*<sup>0</sup> *i dG*<sup>0</sup> *i*

(22)

*<sup>G</sup>*, *Si <sup>I</sup>*, *I i* 0, *G<sup>i</sup>* 0 � �, and we assume that:

*<sup>I</sup>* � *N μSI*

, *σSI* � �, *I*

dom effect vector *bi* <sup>¼</sup> *Si*

� �, *S<sup>i</sup>*

function, with the mean *ϑ* ¼ *μSG* , *μSI*

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

*diag σSG* , *σSI* , *σI*<sup>0</sup> , *σG*<sup>0</sup> ð Þ, so we have Ψ ¼ *μSG* , *μSI*

*<sup>G</sup>* � *N μSG* , *σSG*

ð Þ *p*2, *n*, *γ*, *h*, *σ*1, *σ*2, *σ*<sup>3</sup> .

*<sup>L</sup>*ð Þ *<sup>a</sup>* ð Þ¼ *<sup>θ</sup>*, <sup>Ψ</sup> <sup>Y</sup>*<sup>M</sup>*

with

*i*¼1 ð Þ <sup>2</sup>*<sup>π</sup>* �<sup>2</sup>

> *exp* �X*ni j*¼1

> > *Ai*

approximate estimates ^*θ* and Ψ^.

*3.2.1 Simulation results*

**103**

*lj* <sup>¼</sup> *<sup>Y</sup><sup>i</sup> j* � � *l* � *<sup>Y</sup><sup>i</sup> j*�1 � �

ð *R*4

� 1 <sup>2</sup> *<sup>σ</sup>*�<sup>1</sup> *SG <sup>S</sup><sup>i</sup> <sup>G</sup>* � *μSG* � �<sup>2</sup>

*Si*

#### **Figure 3.**

*At first, glucose and insulin concentrations in the blood are described by two sets of differential equations (see [50]); at a rate* p1*, glucose leaves and enters the glucose space in proportion to the difference between the plasma glucose concentration G(t) and the basal plasma concentration* Gb *that represent the known pre-injection glucose level for each individual. Therefore, the parameter* p1 *represents the glucose's own ability to be eliminated in muscles, liver, and tissues independently of insulin which is called glucose efficiency and denoted by* SG*. Then, the glucose disappears from the glucose space at a rate proportional to insulin concentration in the insulin compartment* X tð Þ *which represents the dynamic of insulin response according to the two rates* p2 *and* p3*. These two parameters represent, respectively, the decreased glucose absorption capacity in tissues and its increased insulin dependency. For insulin secretion I(t), it is secreted by the pancreas independently of the glucose concentration, and proportional to a rate n to its own level already in the body and to the glucose level deferred from a threshold h at a rate γ when* G tð Þ *is above h, the insulin secretion does not only depend on the hyperglycemia level but also to the time spent since glucose injection.*

The insulin sensitivity is defined by combining the two rates *p*<sup>3</sup> and *p*<sup>2</sup> as *SI* ¼ *p*3*=p*2, representing insulin's ability to increase the net glucose utilization [51].

Finally, the stochastic minimal model under the Itô sense re-parameterized by *SG* and *SI* can be rewritten as:

$$dY^i(t) = \begin{pmatrix} -\left(\mathbf{S}\_G^i + X(t)^i\right)\mathbf{G}(t)^i + \mathbf{S}\_G^i\mathbf{G}\_b^i \\\ -p\_2\left(X(t)^i + \mathbf{S}\_I^i\left(I(t)^i - I\_b^i\right)\right) \\\ -n\left(I(t)^i - I\_b^i\right) + \gamma\left(\mathbf{G}(t)^i - h\right)\mathbf{t}\right) \\\ \end{pmatrix} \quad dt + \Sigma dW(t); \quad Y\_0^i = \begin{pmatrix} G\_0^i \\\ 0 \\ I\_0^i \end{pmatrix} \tag{20}$$

where *Y<sup>i</sup>* ðÞ¼ *t G t*ð Þ*<sup>i</sup> X t*ð Þ*<sup>i</sup> I t*ð Þ*<sup>i</sup>* 0 B@ 1 CA and <sup>Σ</sup> is the diagonal diffusion matrix which contains

the unknown constants *σ*1, *σ*2, and *σ*3. The parameters *Si <sup>G</sup>*, *p*2, *S<sup>i</sup> <sup>I</sup>*, *n*, *γ*, *h*, *Gi* <sup>0</sup> and *I i* <sup>0</sup> are unknown in the model and will be estimated. The parameters *Si <sup>G</sup>*, *Si <sup>I</sup>*, *I i* <sup>0</sup> and *G<sup>i</sup>* <sup>0</sup> are assumed random, because they represent individual parameters that change from an individual to another. Each person has its own insulin sensitivity *S<sup>i</sup> <sup>I</sup>* which allows to know if the cells of his body react correctly or not to the insulin and if the insulin produced by the pancreas is sufficient or not, which can make some people with T2D and others without diabetes. Also, for glucose effectiveness *Si G*, which represents the glucose's own ability to be eliminated independently of insulin, it is unique to each individual and changes from a person to another, as well as for the measurement of glucose and insulin concentration. For the rest of the parameters, we consider them fixed since they describe the common side of the glucose-insulin

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

process for the entire population (all individuals). So, we have the following random effect vector *bi* <sup>¼</sup> *Si <sup>G</sup>*, *Si <sup>I</sup>*, *I i* 0, *G<sup>i</sup>* 0 � �, and we assume that:

$$\mathcal{S}\_{\mathcal{G}}^{i} \sim \mathcal{N}(\mu\_{\mathcal{S}\_{\mathcal{G}}}, \sigma\_{\mathcal{S}\_{\mathcal{G}}}), \ \mathcal{S}\_{\mathcal{I}}^{i} \sim \mathcal{N}(\mu\_{\mathcal{S}\_{\mathcal{I}}}, \sigma\_{\mathcal{S}\_{\mathcal{I}}}), \ I\_{0}^{i} \sim \mathcal{N}(\mu\_{I\_{0}}, \sigma\_{I\_{0}}), \ G\_{0}^{i} \sim \mathcal{N}(\mu\_{G\_{0}}, \sigma\_{G\_{0}}).\tag{21}$$

Random effects are assumed to be independent with a multinormal joint density function, with the mean *ϑ* ¼ *μSG* , *μSI* , *μI*<sup>0</sup> , *μG*<sup>0</sup> � � and the covariance matrix *<sup>ϕ</sup>* <sup>¼</sup> *diag σSG* , *σSI* , *σI*<sup>0</sup> , *σG*<sup>0</sup> ð Þ, so we have Ψ ¼ *μSG* , *μSI* , *μI*<sup>0</sup> , *μG*<sup>0</sup> , *σSG* , *σSI* , *σI*<sup>0</sup> , *σG*<sup>0</sup> � � and *<sup>θ</sup>* <sup>¼</sup> ð Þ *p*2, *n*, *γ*, *h*, *σ*1, *σ*2, *σ*<sup>3</sup> .

So, here we deal with a time-inhomogeneous NLME model with SDE describing the glucose-insulin kinetics (see [52] for the implementation of SDE timeinhomogeneous 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 ð Þ *<sup>θ</sup>*, <sup>Ψ</sup> given the observations *<sup>y</sup>* <sup>¼</sup> *<sup>y</sup>*1, … , *<sup>y</sup><sup>M</sup>* � � from model Eq. (20). By using the approximated transition density (Eq. (11)), we get the following approximated likelihood function for model Eq. (20):

*<sup>L</sup>*ð Þ *<sup>a</sup>* ð Þ¼ *<sup>θ</sup>*, <sup>Ψ</sup> <sup>Y</sup>*<sup>M</sup> i*¼1 ð Þ <sup>2</sup>*<sup>π</sup>* �<sup>2</sup> ð Þ *σSGσSIσ<sup>I</sup>*0*σ<sup>G</sup>*<sup>0</sup> �1 2 Y*ni j*¼1 2 ffiffiffiffiffiffiffiffi ΠΔ*<sup>j</sup>* � � <sup>p</sup> �<sup>3</sup> � �½ � *<sup>σ</sup>*1*σ*2*σ*<sup>3</sup> �<sup>1</sup> 2Þ ð *R*4 *exp* �X*ni j*¼1 � �<sup>1</sup> 4Δ*<sup>j</sup>* � 1 *σ*1 *Ai* 1*j* � �<sup>2</sup> þ 1 *σ*2 *Ai* 2*j* � �<sup>2</sup> þ 1 *σ*3 *Ai* 3*j* � �<sup>2</sup>�i � 1 <sup>2</sup> *<sup>σ</sup>*�<sup>1</sup> *SG <sup>S</sup><sup>i</sup> <sup>G</sup>* � *μSG* � �<sup>2</sup> � *<sup>σ</sup>*�<sup>1</sup> *SI <sup>S</sup><sup>i</sup> <sup>I</sup>* � *μSI* � �<sup>2</sup> � *<sup>σ</sup>*�<sup>1</sup> *<sup>I</sup>*<sup>0</sup> *I i* <sup>0</sup> � *μ<sup>I</sup>*<sup>0</sup> � �<sup>2</sup> � *<sup>σ</sup>*�<sup>1</sup> *<sup>G</sup>*<sup>0</sup> *<sup>G</sup><sup>i</sup>* <sup>0</sup> � *μ<sup>G</sup>*<sup>0</sup> � �<sup>2</sup> � �Þ*dSG i dSI i dI*<sup>0</sup> *i dG*<sup>0</sup> *i* (22)

with

The insulin sensitivity is defined by combining the two rates *p*<sup>3</sup> and *p*<sup>2</sup> as *SI* ¼ *p*3*=p*2, representing insulin's ability to increase the net glucose utilization [51]. Finally, the stochastic minimal model under the Itô sense re-parameterized by

*At first, glucose and insulin concentrations in the blood are described by two sets of differential equations (see [50]); at a rate* p1*, glucose leaves and enters the glucose space in proportion to the difference between the plasma glucose concentration G(t) and the basal plasma concentration* Gb *that represent the known pre-injection glucose level for each individual. Therefore, the parameter* p1 *represents the glucose's own ability to be eliminated in muscles, liver, and tissues independently of insulin which is called glucose efficiency and denoted by* SG*. Then, the glucose disappears from the glucose space at a rate proportional to insulin concentration in the insulin compartment* X tð Þ *which represents the dynamic of insulin response according to the two rates* p2 *and* p3*. These two parameters represent, respectively, the decreased glucose absorption capacity in tissues and its increased insulin dependency. For insulin secretion I(t), it is secreted by the pancreas independently of the glucose concentration, and proportional to a rate n to its own level already in the body and to the glucose level deferred from a threshold h at a rate γ when* G tð Þ *is above h, the insulin secretion does not only depend on the*

*SG* and *SI* can be rewritten as:

0

BBBB@

� *Si*

�*n It*ð Þ*<sup>i</sup>*

0

B@

ðÞ¼ *t*

�*p*<sup>2</sup> *X t*ð Þ*<sup>i</sup>*

*<sup>G</sup>* <sup>þ</sup> *X t*ð Þ*<sup>i</sup>* � �

*hyperglycemia level but also to the time spent since glucose injection.*

*Numerical Modeling and Computer Simulation*

� *I i b* � �

> *G t*ð Þ*<sup>i</sup> X t*ð Þ*<sup>i</sup> I t*ð Þ*<sup>i</sup>*

<sup>þ</sup> *<sup>S</sup><sup>i</sup> <sup>I</sup> I t*ð Þ*<sup>i</sup>*

1

the unknown constants *σ*1, *σ*2, and *σ*3. The parameters *Si*

unknown in the model and will be estimated. The parameters *Si*

an individual to another. Each person has its own insulin sensitivity *S<sup>i</sup>*

T2D and others without diabetes. Also, for glucose effectiveness *Si*

*G t*ð Þ*<sup>i</sup>*

<sup>þ</sup> *<sup>γ</sup> G t*ð Þ*<sup>i</sup>*

� � � �

<sup>þ</sup> *<sup>S</sup><sup>i</sup> GG<sup>i</sup> b*

� *I i b*

� *h* � �

assumed random, because they represent individual parameters that change from

to know if the cells of his body react correctly or not to the insulin and if the insulin produced by the pancreas is sufficient or not, which can make some people with

sents the glucose's own ability to be eliminated independently of insulin, it is unique to each individual and changes from a person to another, as well as for the measurement of glucose and insulin concentration. For the rest of the parameters, we consider them fixed since they describe the common side of the glucose-insulin

*t*

1

CCCCA

*dt* <sup>þ</sup> <sup>Σ</sup>*dW t*ð Þ; *<sup>Y</sup><sup>i</sup>*

*<sup>G</sup>*, *p*2, *S<sup>i</sup>*

CA and <sup>Σ</sup> is the diagonal diffusion matrix which contains

<sup>0</sup> ¼

*<sup>I</sup>*, *n*, *γ*, *h*, *Gi*

*<sup>G</sup>*, *Si <sup>I</sup>*, *I i* <sup>0</sup> and *G<sup>i</sup>*

*Gi* 0 0 *I i* 0

1

CA ,

(20)

<sup>0</sup> are

<sup>0</sup> and *I i* <sup>0</sup> are

*<sup>I</sup>* which allows

*G*, which repre-

0

B@

*dY<sup>i</sup>* ðÞ¼ *t*

**102**

**Figure 3.**

where *Y<sup>i</sup>*

$$\boldsymbol{A}\_{lj}^{i} = \left[ \left( \mathbf{Y}\_{j}^{i} \right)\_{l} - \left( \mathbf{Y}\_{j-1}^{i} \right)\_{l} - \mu\_{l} \left( \mathbf{Y}\_{j-1}^{i}, \mathbf{t}\_{j-1}, \theta, \mathbf{b}^{i} \right) \Delta\_{j}^{i} \right]. \tag{23}$$

We have no closed form solution to this integral, so exact estimators of *θ* and Ψ are unavailable. So, we use the Laplace approximation method described in Section 2 to obtain a closed form approximation to the log-likelihood function *log L*ð Þ *<sup>a</sup>* ð Þ *<sup>θ</sup>*, <sup>Ψ</sup> � � for the model Eq. (20); then by applying the GA, we get the approximate estimates ^*θ* and Ψ^.

### *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 ð Þ� *ni* þ 1 *M* from Eq. (20) using the Euler-Maruyama scheme [45] with a step size of 10�<sup>3</sup> and a set of true parameters that are chosen according to [53, 54] representing the normal range of parameters values to simulate healthy subjects (without diabetes), with M being the number of units and *ni* being the number of observations or repeated measurements collected on each unit i.

For each data from 5000 generated data sets, we estimate ð Þ *θ*, Ψ by applying the proposed method. We first assume that the number of repeated measurements


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 *σSI* , *σ<sup>I</sup>*<sup>0</sup> , and *σ<sup>G</sup>*<sup>0</sup> ; nevertheless, the results are more satisfactory when the sample size M is large for all parameters. As mentioned before, we cannot ignore the limitations of this method, and the results could be improved by eliminating these limits and improving this approach, but for the given assumptions and tools, we can say that the results are still satisfactory even for a sample of *M* ¼ 10 with 20 measures taken on each subject. So, we can conclude that we can rely on the results obtained from small samples with a small number of measurement repetitions of at least 20 observations. In addition, although we can use a relatively small number of measurements, we need to know how to choose the time points to perform the measurements in the blood, as this could, physiologically, affect selected observations and results. Thus, it is specified here that the essential task is to know how to choose the good moments of measurement after the injection, even in small numbers, chosen according to medical knowledge. **Figure 4** shows that, in the case of ð Þ¼ *M*; *m* ð Þ 40; 60 , the empirical distribution of the most approximated

*SMM: empirical distribution of population parameter estimates obtained using (18) for (M, m) = (40,60).*

*Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential…*

*DOI: http://dx.doi.org/10.5772/intechopen.90751*

estimates seems to be reasonably close to a normal distribution.

tion to approximate the transition density of a stochastic process.

**4. Conclusions**

**105**

**Figure 4.**

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 approxima-

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

#### **Table 2.**

*Approximated maximum likelihood estimates and standard deviation from simulations of model Eq. (20), using the approximated transition density Eq. (11) with large and small DATA.*

collected on each unit is constant *ni* <sup>¼</sup> *<sup>m</sup>* and <sup>Δ</sup>*<sup>i</sup> <sup>j</sup>* ¼ Δ for all 1≤*i* ≤ *M* and 1≤*j*≤ *ni*; then we apply the GA after choosing the good algorithm parameters (N, EN, SR, CP, MP), and we get 5000 sets of parameters estimates. We repeat this for large and small data with different possibilities of repetition of the experiment ð Þ¼ *m*, *M* ð Þ 40, 60 ; 40, 10 ð Þ and 10, 20 ð Þ; for each parameter the sample mean and standard deviation are reported in **Table 2**. The simulation study on small data is treated in order to see if the size of the sample influences on the results and if the number of measures taken over time has a negligible effect or not, in other words, to see if it is possible to select only the essential measuring moments without repeating the measurements several times to well simulate a subject. We treat this issue in relation to our model and its study context, since in epidemiology the availability of data (measurements) at any point of time is an interesting constraint. We note that quantities of *Gb* and *Ib* are randomly simulated from the normal range of healthy subjects.

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 *Simulation and Parametric Inference of a Mixed Effects Model with Stochastic Differential… DOI: http://dx.doi.org/10.5772/intechopen.90751*

**Figure 4.** *SMM: empirical distribution of population parameter estimates obtained using (18) for (M, m) = (40,60).*

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 *σSI* , *σ<sup>I</sup>*<sup>0</sup> , and *σ<sup>G</sup>*<sup>0</sup> ; nevertheless, the results are more satisfactory when the sample size M is large for all parameters. As mentioned before, we cannot ignore the limitations of this method, and the results could be improved by eliminating these limits and improving this approach, but for the given assumptions and tools, we can say that the results are still satisfactory even for a sample of *M* ¼ 10 with 20 measures taken on each subject. So, we can conclude that we can rely on the results obtained from small samples with a small number of measurement repetitions of at least 20 observations. In addition, although we can use a relatively small number of measurements, we need to know how to choose the time points to perform the measurements in the blood, as this could, physiologically, affect selected observations and results. Thus, it is specified here that the essential task is to know how to choose the good moments of measurement after the injection, even in small numbers, chosen according to medical knowledge. **Figure 4** shows that, in the case of ð Þ¼ *M*; *m* ð Þ 40; 60 , the empirical distribution of the most approximated estimates seems to be reasonably close to a normal distribution.

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.
