Stochastic Differential Equation Models

• Stochastic Differential Equation Models are mathematical frameworks that combine deterministic dynamics with random noise to accurately describe evolving systems.
• They extend classic models by incorporating random effects to distinguish intra-subject variability from inter-subject differences in experimental data.
• Techniques like Hermite expansion and Laplace approximation enable efficient parameter estimation even in high-dimensional, nonlinear systems.

Stochastic differential equation (SDE) models are mathematical frameworks that describe systems evolving in time under both deterministic dynamics and random fluctuations. These models are particularly suited to capturing the intrinsic randomness inherent in physical, biological, and engineered systems, as well as to modeling between-subject and within-subject variation in repeated experimental data. SDEs are commonly applied in fields including pharmacokinetics/pharmacodynamics (PK/PD), biomedical growth modeling, neuroscience, and population biology, and have become foundational tools where process noise and heterogeneity are both scientifically and statistically meaningful.

1. Stochastic Differential Mixed-Effects Models (SDMEMs): Motivation and Definition

To accommodate both the randomness of individual trajectories and variability across experimental units (such as individuals in a population or repeated experiments), stochastic differential mixed-effects models (SDMEMs) were developed. These extend classic SDEs by incorporating random effects—subject- or experiment-specific parameters—to statistically distinguish inter-individual (or between-group) variation from intra-individual (within-group) process stochasticity.

The general form of an SDMEM for the $i$th individual ($i = 1, \ldots, M$), with process dimension $d$, is:

$$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$

where:

• $\mathbf{X}_t^i$ is the state vector of subject $i$ at time $t$.
• $\theta$ is a vector of fixed effects (population parameters).
• $\mathbf{b}^i$ is a vector of random effects (unit-specific parameters), possibly non-Gaussian and multidimensional.
• $\mu$ and $i = 1, \ldots, M$0 are, respectively, the drift and diffusion (noise) functions.
• $i = 1, \ldots, M$1 is a standard multivariate Brownian motion, independent across $i = 1, \ldots, M$2.

Such models can be observed at subject-specific sampling times, potentially unevenly spaced, and are directly relevant in scenarios such as modeling patient drug response trajectories, plant or animal growth curves, and neural activity subject to both biological noise and heterogeneity.

2. Likelihood-based Parameter Estimation: Structure and Computational Challenges

The statistical inference goal is the joint estimation of $i = 1, \ldots, M$3 and $i = 1, \ldots, M$4 (the parameters of the random effect distribution). The marginal likelihood of the full observed data (across all individuals) is:

$i = 1, \ldots, M$5

where:

• $i = 1, \ldots, M$6 is the likelihood of data for subject $i = 1, \ldots, M$7, given random effects and population parameters.
• $i = 1, \ldots, M$8 is the density of the random effects.

$i = 1, \ldots, M$9 typically decomposes as a product of transition densities between observation times, but closed-form expressions for these transition densities are rarely available except for special cases (e.g., Ornstein-Uhlenbeck processes).

Major computational issues arise due to:

• Intractable transition densities for nonlinear/multidimensional SDEs.
• High-dimensional integration over random effects (especially with non-normal, multidimensional, or correlated effects).

3. Transition Density Approximations and Laplace Integration

Hermite Expansion for Transition Density

For SDEs where a transformation to unit diffusion (via the Lamperti transform $d$0) exists, the transition density over a time increment $d$1 can be approximated by a Hermite expansion:

$d$2

This yields a closed-form, polynomial approximation, facilitating efficient computation. Truncation at $d$3 is often sufficient for practical accuracy. The Hermite expansion generally outperforms lower-order approximations like Euler-Maruyama in low-frequency data or for highly nonlinear systems.

Laplace Approximation for Marginalization over Random Effects

High-dimensional integrals over random effects are performed using the Laplace approximation, which leverages a Taylor expansion of the integrand about its mode:

$d$4

with:

• $d$5,
• $d$6 the maximizing value (MAP estimate) for unit $d$7,
• $d$8 the Hessian of $d$9 with respect to $$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$0.

The resulting marginal likelihood for parameter estimation is then maximized via a two-level optimization:

1. Inner loop: For each unit $$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$1, maximize $$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$2 to obtain $$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$3 (using trust-region or Newton-type methods).
2. Outer loop: Maximize the Laplace-approximated log-likelihood over $$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$4 and $$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$5 (using derivative-free optimizers, e.g. Nelder-Mead, due to computational cost).

Automatic differentiation (AD) is used throughout to efficiently compute gradients and Hessians required for optimization and for evaluating the Laplace approximation.

4. Simulation Studies and Applied Examples

A. Stochastic Logistic Growth in Orange Trees

• Model: Stochastic logistic equation with state-dependent diffusion and random effects on both the asymptotic maximum and growth rate:

$$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$6

• Key results: Hermite expansion yields more accurate fixed and random effect parameter estimates compared with Euler-Maruyama, especially for sparse longitudinal data.

B. Two-Dimensional Ornstein-Uhlenbeck (OU) Process

• Model: Bivariate OU SDE with random effects on drift coefficients, relevant to tissue microvascularization and neurobiology.
• Results: With exact transition densities available, the full inferential pipeline (Hermite, Laplace) recovers population and individual-level parameters robustly, including gamma-distributed random effects.

C. Square-root Process (CIR/Feller)

• Model: SDE with square-root diffusion and random effects, as used in finance and neuroscience.
• Result: The methodology accommodates non-normal random effect distributions (log-normal, Beta), and the Hermite + Laplace approach recovers both fixed and random effect components.

5. Software Implementation Considerations

• Automatic differentiation enables efficient calculation of required derivatives, making the method practical for high-dimensional models.
• Parallelization across units (subjects) is possible, as random effect integrations and optimizations are independent.
• Trust-region methods for inner optimization and derivative-free methods for outer optimization lead to robust and convergent parameter estimation.
• Where possible, symbolic Hessian calculation is preferred for accuracy and speed.

Limitation: The default modeling approach does not formally accommodate measurement error. When measurement noise is substantial relative to system/process noise, additional methods (e.g., extensions involving state-space/Kalman filter models) are required.

6. Relevance for Biomedical and PK/PD Modeling

SDMEMs, with the described estimation methodology, are directly relevant to:

• Pharmacokinetics/Pharmacodynamics (PK/PD): Modeling dynamics of drug concentration/effect, capturing both inter-individual variability (due to random effects) and intra-individual stochasticity.
• Biomedical Growth: Describing stochastic growth under biological heterogeneity.
• Neuroscience: Analysis of firing rates and subcellular dynamics where channel noise and subject-level heterogeneity both play substantial roles.

Parameter uncertainty can be effectively propagated, and the method allows realistic modeling of biological variance components, crucial for optimal dosage regimen design and for understanding biological mechanisms.

7. Key Formulae and Model Structures

General SDMEM SDE

$$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$7

Marginal Likelihood with Laplace Approximation

$$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$8

Example Transition Density (2D OU)

$$d\mathbf{X}_t^i = \mu(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; dt + \sigma(\mathbf{X}_t^i, \theta, \mathbf{b}^i) \; d\mathbf{W}_t^i, \quad \mathbf{X}_0^i = \mathbf{x}_0^i$$9

where

$\mathbf{X}_t^i$0

8. Concluding Summary

The practical estimation of high-dimensional stochastic differential mixed-effects models is made feasible and statistically robust through the combination of Hermite expansion-based closed-form approximations for transition densities, Laplace approximation for integrating over random effects, and the use of automatic differentiation for computational efficiency. This framework addresses the two critical axes of biomedical and PK/PD inference: accurate representation of noise (within-individual and system), and flexible modeling of population heterogeneity. Key simulation studies confirm the method’s applicability, accuracy, and relevance across a range of applied domains, especially where small sample sizes and complex process noise structures are encountered.