# Stochastic differential equations based models

(diff) ← Version précédente | Voir la version courante (diff) | Version suivante → (diff)

## Introduction

Diffusion models are known to be a relevant tool for modeling stochastic dynamic phenomena, and are widely used in various fields including finance, physics, biology, physiology and control. In a population approach, a mixed-effects diffusion model describes each individual series of observations using a system of stochastic differential equations (SDE) while also taking into account variability between individuals.

For the sake of simplicity we will consider first a diffusion model for a single individual, and illustrate it with a very general dynamical system with linear transfers and PK examples. We will then show that the extension to mixed diffusion models is fairly straightforward.

Note that the conditional distribution $\qcypsi$ of the observations usually does not have a closed-form expression. When the underlying system is a Gaussian linear dynamical one, the conditional pdf of the observations, $\pcypsi(y_i|\psi_i)$ can be computed using the Kalman filter (KF). When the system is not linear, the extended Kalman filter (EKF) provides an approximation of the conditional pdf.

## Diffusion model

We assume that one diffusion trajectory is observed with noise at discrete time points $t_1<\ldots<t_j<\ldots<t_n$. Let us note $(X(t),t>0) \in \Rset^d$ the underlying dynamical process and $y_j \in \Rset$ a noisy function of $X(t_j)$, $j=1,\ldots,n$. The general form of the diffusion model is given by:

 $$\left\{ \begin{array}{lll} dX(t) &=& b(X(t),\psi)dt + \gamma(X(t),\psi)dW(t)\\[0.2cm] y_{j} &=& c(X(t_{j}),\psi) + \varepsilon_{j} \\[0.2cm] \varepsilon_{j} &\underset{i.i.d.}{\sim}& \mathcal{N}(0,a^2(\psi)), \quad j=1,\ldots,n , \end{array} \right.$$ (1)

with the initial condition $X(t_1) = x \in \Rset^d$. Here, $(W(t),t>0)$ is a standard Wiener process in $\Rset^d$ and $\varepsilon_j \in \Rset$ represents the measurement error occurring at the $j^{\mathrm{th}}$ observation, independent of $W(t)$. The measurement function $c: \ \Rset^d \times \Rset^p \rightarrow \Rset$, the drift function $b: \ \Rset^d \times \Rset^p \rightarrow \Rset^d$ and the diffusion function $\gamma: \ \Rset^d \times \Rset^p \rightarrow \mathcal{M}_d(\Rset)$, where $\mathcal{M}_d(\Rset)$ is the set of $d \times d$ matrices with real elements, are known functions that depend on an unknown parameter $\psi \in \Rset^p$.

We can in fact consider an SDE-based model as a ODE-based one with a stochastic component.

Example: IV bolus with linear elimination

The ordinary differential equation
 $$dA_c(t) = -k A_c(t) dt$$ (2)

is usually used to describe the kinetics of a drug administered by rapid injection (IV bolus) into plasma. In bolus-specific compartmental models, plasma is treated as the single compartment of the human body. $A_c(t)$ represents the amount of a drug ingredient in plasma at time $t$ after injection, and $k$ is the elimination rate constant. The figure below displays the typical evolution of the amount found in the central compartment when $k=4$.

 Drug concentration evolution for ODE diffusion example

Imagine now that we aim to describe the evolution of the drug amount over time by means of stochastic differential equations rather than ordinary differential equations, in order to better describe the intra-individual variability of the observed process. We can assume for example that the system (2) is randomly perturbed by an additive Wiener process:

 $$dA_c(t) = -k A_c(t) dt + \gamma dW(t).$$ (3)

The figure below displays four kinetics for the amount in the central compartment, simulated from this model with $k=4$ and $\gamma=2$.

 Drug concentration evolution for SDE diffusion example

These kinetics are clearly stochastic. Nevertheless, they are not realistic because:

• they give an overly erratic description of the evolution of the drug concentration within the compartments of the human body.
• they do not comply with certain constraints on biological dynamics (sign, monotony).

A more relevant model might consider that some parameters of the model randomly fluctuate over time, rather than the observed variable itself, modeling for example the elimination rate "constant" $k$ as a stochastic process $k(t)$ that randomly varies around a typical value $k^\star$.

More generally, we can describe the fluctuations within a linear dynamical systems by considering the transfer rates, described below, as diffusion processes rather than the observed processes themselves.

## Diffusion models for dynamical systems with linear transfers

Dynamical systems have applications in many fields. They can be used to model viral dynamics, population flows, interactions between cells, and drug pharmacokinetics. Dynamical systems involving linear transfers between different entities are usually modeled by means of a system of ODEs with the following general form:

 $$dA(t) = K\, A(t)dt,$$ (4)

where $A(t)$ is a vector whose $l^{\textrm{th}}$ component represents the condition of the $l^{\textrm{th}}$ entity at time $t$ and $K=(K_{l,l^\prime} \, 1\leq l , l^\prime \leq d)$ a deterministic matrix defined as:

 $$\left\{ \begin{array}{ll} K_{l,l^\prime} = k_{l,l^\prime} & \textrm{if} \quad l \neq l^\prime\\ K_{l,l} = - k_{l,0} - \sum_{l^\prime} k_{l,l^\prime} , \end{array} \right.$$ (5)

where $k_{l,l^\prime}$ represents the transfer rate from entity $l$ to entity $l^\prime$, and $k_{l,0}$ the elimination rate from entity $l$. An example of such a dynamical system with $3$ components is schematized below.

 A dynamical system with $3$ components (circles) and linear transfers between components (arrows)

In this particular example, matrix $K$ would be defined as

$$K = \begin{pmatrix} -k_{10} -k_{12} -k_{13} & k_{21} & k_{31}\\ k_{12} & -k_{20} -k_{21} -k_{23} & k_{32}\\ k_{13} & k_{23} & -k_{30} -k_{31} -k_{32} \end{pmatrix}.$$

The model defined by equations (4) and (5) is a deterministic model which assumes that transfers take place at the same rate at all times. This is often a restrictive assumption since in reality, dynamical systems usually exhibit some random behavior. It is therefore reasonable to consider that transfers are not constant but randomly fluctuate over time. This new assumption leads to the following dynamical system:

 $$dA(t) = K(t)A(t)dt,$$ (6)

where $K$ has the same structure as in (5) but now some components $k_{l,l^\prime}$ are stochastic processes which take non-negative values and randomly fluctuate around a typical value $k_{l,l^\prime}^\star$.

Let us now illustrate the construction of such diffusion models using some specific examples in pharmacokinetics.

Example 1: IV bolus administration with stochastic linear elimination

We will first extend the ODE based model defined in (2) by assuming that $k$ is a diffusion process which takes non-negative values and fluctuates around a typical value $k^\star$.

In this example, non-negativity of $k(t)$ is ensured by defining the logarithm of the transfer rate as an Ornstein-Uhlenbeck diffusion process:

$$d\log k(t) = - \alpha \left( \log k(t) - \log k^\star \right) dt + \gamma d W(t),$$

where $W$ is a standard one-dimensional Wiener process. This results in the following diffusion system:

$$dX(t) = b(X(t))dt + \gamma(X(t))dW(t),$$

where

$$X(t) = \begin{pmatrix} A_c(t) \\ \log k(t) \end{pmatrix}, \ \ \ \ b(x) = \begin{pmatrix} -x_1 \exp(x_2) \\ -\alpha (x_2-\log k^{\star}) \end{pmatrix}, \ \ \ \ \gamma(x) = \begin{pmatrix} 0 & 0 \\ 0 & \gamma \end{pmatrix}.$$

Note that in this specific example, the Jacobian matrix of the drift function $b$ has a simple form:

$$B(x)=\begin{pmatrix} - \exp(x_2) & -x_1 \exp(x_2)\\ 0 & -\alpha \end{pmatrix}.$$

The two figures below display four simulated processes $k(t)$ and the associated amount processes $A_c(t)$.

We measure the concentration at times $(t_{j}, 1\leq j \leq n)$:

$$y_j = \displaystyle{\frac{A_c(t_{j})}{V} } + a \, \teps_j .$$

The parameter vector of the model is therefore $\psi = (V, k^\star, \alpha, \gamma, a)$. We see in this example that the simulated kinetics are much more realistic than those obtained with the previous model, because:

• the elimination rate process $k(t)$ is a stochastic process that takes non-negative values,
• even though the amount process is stochastic, it is smooth and decreases monotonically with time.

Example 2: Oral administration with first-order absorption and stochastic linear elimination

Oral PK models with first-order absorption and linear elimination are widely used to describe the time-course of a drug orally administered to a unique compartment of the human body. The drug is administrated in a depot compartment, absorbed by the central compartment with absorption rate $k_a$ and eliminated with elimination rate $k_e$. Such a model is described by the following system of ODEs:
 $$\displaystyle{ \frac{d}{dt} } \begin{pmatrix} A_d(t) \\ A_c(t) \end{pmatrix} \ \ = \ \ \begin{pmatrix} -k_a & 0\\ k_a & -k_e\end{pmatrix} \begin{pmatrix} A_d(t) \\ A_c(t) \end{pmatrix},$$ (7)

where $A_d(t)$ and $A_c(t)$ respectively represent the amounts of drug at time $t$ in the depot and central compartments. Assume now that the elimination constant is driven by a stochastic process, solution to the stochastic differential equation

$$d k_e(t) = - \alpha (k_e - k_e^\star ) dt + \gamma \sqrt{k_e(t)} dW(t),$$

where $W$ is a standard one-dimensional Wiener process. Then (7) becomes:

$$dX(t) = b(X(t))dt + \gamma(X(t))dW(t).$$

Here,

$$X(t)= \begin{pmatrix} A_d(t) \\ A_c(t) \\ k_e(t) \end{pmatrix}, \ \ \ \ b(x) = \begin{pmatrix} -k_a x_1 \\ k_a x_1 -x_3 x_2 \\ -\alpha(x_3-k_e^\star ) \end{pmatrix}, \ \ \ \ \gamma(x) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0\\ 0 & 0 & \gamma \sqrt{x_3}\end{pmatrix} ,$$

and the parameter vector of the model is $\psi = (V, k_a, k^\star, \alpha, \gamma, a) .$

In both examples, the diffusion model can be easily extended to a population approach by defining the system's parameters $\psi$ as an individual random vector.

## Mixed-effects diffusion models

Let us now consider model (1) with observations coming from several subjects. An adequate adaptation of model (1) in such a context consists of considering as many dynamical systems as individuals, and defining the parameters of the individual dynamical systems as independent random variables, in such a way to correctly reflect variability between the different trajectories. To standardize notation, we consider $N$ different subjects randomly chosen from a population and note $n_i$ the number of observations for individual $i$, so that $t_{i1}<\ldots<t_{i,n_i}$ are subject $i$'s observation time points. $(X_i(t),t>0) \in \Rset^d$ and $y_{ij} \in \Rset$ will respectively denote individual $i$'s diffusion and the observation $X_i(t_{ij})$. The $y_{ij}$, $i=1,\ldots,N$, $j=1,\ldots,n_i$ are governed by a mixed-effects model based on a $d$-dimensional real-valued system of stochastic differential equations with the general form:

 $$\left\{ \begin{array}{l} dX_i(t) = b(X_i(t),\psi_i)dt + \gamma(X_i(t),\psi_i)dW_i(t),\\[0.2cm] y_{ij} = c(X_i(t_{ij}),\psi_i) + \teps_{ij},\\[0.2cm] \teps_{ij} \underset{i.i.d.}{\sim} \mathcal{N}(0,a^2(\psi_i)) \; , \; j=1,\ldots, n_i \; , \; i=1,\ldots,N,\\ \end{array} \right.$$ (8)

with initial condition $X_i(t_1) = x_{i1} \in \Rset^d$ for $i=1,\ldots,N$. The $\psi_i$'s are unobserved independent $d$-dimensional random subject-specific parameters, drawn from a distribution $\qpsi$ which depends on a set of population parameters $\theta$, $(W_1(t),t>0), \ldots, (W_N(t),t>0)$ are standard independent Wiener processes, and the $\teps_{ij}$ are independent Gaussian random variables representing residual errors such that the $\psi_i$, $W_i$ and $\teps_{ij}$ are mutually independent. The measurement function $c$, the drift function $b$ and the diffusion function $\gamma$ are known functions that are common to the $N$ subjects and depend on the unknown parameters $\psi_i$.

Assuming that the $N$ individuals are independent, the joint pdf is given by:

 $$\pcypsi(y_1,\ldots,y_N | \psi_1,\ldots,\psi_N) = \prod_{i=1}^{N}\pcyipsii(y_i | \psi_i).$$ (9)

Computing the conditional distribution $\pcyipsii$ of the observations for any individual $i$ requires here to compute the conditional distribution of each observation given the past:

$$\begin{eqnarray} \pcyipsii(y_i | \psi_i) &=& \pyipsiONE(y_{i1} | \psi_i)\prod_{j=2}^{n_i} p(y_{i,j} | y_{i,1},\ldots,y_{i,j-1} | \psi_i) . \end{eqnarray}$$

Except in some very specific classes of mixed-effects diffusion models, the transition density $\pmacro(y_{i,j}|y_{i,1},\ldots,y_{i,j-1} | \psi_i)$ does not have a closed-form expression since it involves the transition densities of the underlying diffusion processes $X_i$. When the underlying system is a Gaussian linear dynamical system, this density is a Gaussian density whose mean and variance can be computed using the Kalman filter. When the system is not linear, a first solution consists in approximating this density by a Gaussian density and using the extended Kalman filter for quickly computing the mean and the variance of this density. On the other hand, particle filters do not make any approximations of the transition density, but are very demanding in terms of simulation volume and computation time.

## Bibliography

Delattre, M., Lavielle, M. - Coupling the SAEM algorithm and the extended Kalman filter for maximum likelihood estimation in mixed-effects diffusion models

Statistics and Its Interface ,2013
Ditlevsen, S., De Gaetano, A. - Mixed Effects in Stochastic Differential Equation Models
REVSTAT Statistical Journal 3:137-153,2005
Donnet, S., Samson, A. - Parametric Inference for Mixed Models Defined by Stochastic Differential Equations
ESAIM: Probability and Statistics 12:196-218,2008
Doucet, A., Johansen, A. M. - A tutorial on particle filtering and smoothing: Fifteen years later
Oxford Handbook of Nonlinear Filtering ,2011
Klim, S., Mortensen, S. B., Kristensen, N. R., Overgaard, R. V., Madsen, H. - Population stochastic modelling (PSM)-an R package for mixed-effects models based on stochastic differential equations
Computer methods and programs in biomedicine 94:279-289,2009
Kristensen, N. R., Madsen, H., Ingwersen, S. H. - Using Stochastic Differential Equations for PK/PD Model Development
Journal of Pharmacokinetics and Pharmacodynamics 32:109-141,2005
Mazzoni, T. - Computational aspects of continuous-discrete extended Kalman-filtering
Computational Statistics 23:519-39,2008
Mortensen, S., Klim, S. - Population Stochastic Modelling (PSM): Model definition, description and examples
,2008
http://www2.imm.dtu.dk/projects/psm/
Mortensen, S. B., Klim, S., Dammann, B., Kristensen, N. R., Madsen, H., Overgaard, R. V. - A Matlab framework for estimation of NLME models using stochastic differential equations - Applications for estimation of insulin secretion rates
Journal of Pharmacokinetics and Pharmacodynamics 34:623-642,2007
Overgaard, R. V., Jonsson, N., Tornøe, C. W., Madsen, H. - Non-Linear Mixed-Effects Models with Stochastic Differential Equations: Implementation of an Estimation Algorithm
Journal of Pharmacokinetics and Pharmacodynamics 32:85-107,2005
Picchini, U., De Gaetano, A., Ditlevsen, S. - Stochastic Differential Mixed-Effects Models
Scandinavian Journal of Statistics 37:67-90,2010
Picchini, U., Ditlevsen, S. - Practical Estimation of High Dimensional Stochastic Differential Mixed-Effects Models
Computational Statistics and Data Analysis 55(3):1426-1444,2011
Tornøe, C. W., Overgaard, R. V., Agersø, H., Nielsen, H. A., Madsen, H., Jonsson, E. N. - Stochastic Differential Equations in NONMEM: Implementation, Application, and Comparison with Ordinary Differential Equations
Pharmaceutical Research 22:1247-1258,2005