# Variational autoencoder for inference of nonlinear mixed effect models based on ordinary differential equations

Zhe Li, Mélanie Prague, Rodolphe Thiébaut, Quentin Clairon

► To cite this version:

Zhe Li, Mélanie Prague, Rodolphe Thiébaut, Quentin Clairon. Variational autoencoder for inference of nonlinear mixed effect models based on ordinary differential equations. 2026. hal-05472212v3

HAL Id: hal-05472212

https://hal.science/hal-05472212v3

Preprint submitted on 10 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Variational autoencoder for inference of nonlinear mixed effect models based on ordinary differential equations

Zhe Li, Mélanie Prague, Rodolphe Thiébaut, Quentin Clairon

Univ. Bordeaux, INSERM BPH, U1219, Inria SISTM team, VRI, France

email: zhe.li@u-bordeaux.fr, melanie.prague@inria.fr,

rodolphe.thiebaut@u-bordeaux.fr, quentin.clairon@u-bordeaux.fr

**SUMMARY:** We propose a variationl autoencoder approach for identifiability-aware parameter estimation in nonlinear mixed-effects models based on ordinary differential equations (NLME-ODEs) using longitudinal data from multiple subjects. In moderate dimensions, likelihood-based inference via the stochastic approximation EM algorithm (SAEM) is widely used, but it relies on Markov Chain Monte-Carlo (MCMC) to approximate subject-specific posteriors. As model complexity increases or observations per subject are sparse and irregular, performance often deteriorates due to a complex, multimodal likelihood surface which may lead to MCMC convergence difficulties. We instead estimate parameters by maximizing the evidence lower bound (ELBO), a regularized surrogate for the marginal likelihood. A shared encoder amortizes inference of subject-specific random effects by avoiding per-subject optimization and the use of MCMC. Beyond pointwise estimation, we quantify parameter uncertainty using observed-information-based variance estimator and verify that practical identifiability of the model parameters is not compromised by nuisance parameters introduced in the encoder. We evaluated the method in three simulation case studies (pharmacokinetics, humoral response to vaccination, and TGF-β activation dynamics in asthmatic airways) and on a real-world antibody kinetics dataset, comparing against SAEM baselines.

**KEY WORDS:** Longitudinal population data analysis; Neural Networks; Variational bayesien inference;
# 1. Introduction

Ordinary differential equation (ODE) models are widely used in systems biology, mathematical medicine, virology, vaccinology, and pharmacometrics due to their interpretability and predictive power. However, classical ODEs typically describe an average trajectory and do not capture inter-individual variability in longitudinal population data, which are often sparse, irregularly sampled, and noisy. To address this, ODEs are commonly embedded in nonlinear mixed-effects (NLME) models, yielding NLME-ODEs [14], where subject-specific random effects capture heterogeneity while preserving a shared mechanistic structure.

More formally, NLME-ODEs consider a population of $n$ subjects. For subject $i \in \{1, \dots, n\}$, observations are collected at (possibly irregular) times $t_{ij}$, $j = 1, \dots, n_i$, where $n_i$ denotes the number of measurements for subject $i$; we write $Y_{ij}$ for the measurement taken at time $t_{ij}$. The dynamics of the subject $i$ is described by a $D_X$-dimensional ODEs:

$$ \left\{ \begin{array}{l} \frac{dX_i(t)}{dt} = f_{\theta_i}(t, X_i(t)), \\ X_i(0) = x_{i,0}. \end{array} \right. $$

where $f_{\theta_i}$ is a $D_X$-dimensional vector field defining the system dynamics, $\theta_i$ is a vector of subject-specific parameters, and $x_{i,0}$ denotes the initial condition. Inter-individual variability is modeled through a mixed-effects formulation $\theta_i = g(\bar{\theta}, b_i)$, where $\bar{\theta}$ represents the population-level mean and $b_i$ denotes the subject-specific deviation around this mean, also known as the random effects. Here, $g(\cdot)$ is a (possibly nonlinear) transformation function, for example through additive or log-normal parameterizations to enforce constraints such as positivity. We assume $b_i \sim p_\phi(b_i) = \mathcal{N}(0, \Omega)$, where $b_i = (b_{i1}, \dots, b_{id_b})^\top \in \mathbb{R}^{D_b}$ is the $D_b$-dimensional vector of random effects and $\Omega \in \mathbb{R}^{D_b \times D_b}$ is a symmetric positive-definite covariance matrix. To guarantee positive definiteness, we parameterize $\Omega$ via its Cholesky factor $\Omega = L_\Omega L_\Omega^\top$.

Our goal is to estimate the population parameters $\phi = (\bar{\theta}, L_\Omega, L_\Sigma)$ of dimension $D_\phi$ from
the noisy and discrete observations $Y_{ij}$ generated according to the following observation model:

$$Y_{ij} = h_{\theta_i}(X_i(t_{ij})) + \epsilon_{ij}, \quad \epsilon_{ij} \sim \mathcal{N}(0, \Sigma). \qquad (1)$$

Here, $h_{\theta_i}$ is the observation function mapping the $i$-th subject state represented by the ODEs solution $X_i$ to the measure $Y_{ij}$ and $\Sigma \in \mathbb{R}^{D_{n_i} \times D_{n_i}}$ is the covariance matrix of the measurement noise which is parameterized by its Cholesky factor $\Sigma = L_{\Sigma}L_{\Sigma}^{\top}$.

As in classical inference problems, NLME-ODEs parameters can be estimated by maximum likelihood, marginalizing over the random effects $b_i$ viewed as latent variables. However, maximizing the marginal likelihood requires repeated integration over $b_i$, which is generally intractable in nonlinear ODE models. Early approaches such as first-order (FO) and first-order conditional estimation (FOCE) rely on local Taylor/Laplace approximations to obtain a tractable likelihood [20]. These approximations can be inaccurate under strong nonlinearities, non-Gaussian random effects, or sparse individual data, settings common in NLME-ODEs; adaptive Gaussian quadrature can improve accuracy [22] but becomes computationally prohibitive as the number of random effects grows. A widely used method is the stochastic approximation EM (SAEM) algorithm [9; 16], which avoids explicit numerical integration by sampling random effects from their conditional distribution (typically via MCMC) and updating parameters through stochastic approximation. SAEM is often more accurate and stable than FO/FOCE or quadrature for nonlinear models and is implemented in Monolix and Pumas [30; 26]. Nonetheless, for sparse data and/or high-dimensional random effects, SAEM may become computationally expensive and less reliable, motivating scalable variational alternatives.

SAEM inherits several limitations from its reliance on MCMC in the E-step. In sparse settings, the conditional posterior of the random effects is weakly informed by the data, which can lead to slow mixing and poor exploration. As the dimension of the random effects
increases, MCMC becomes increasingly expensive and more sensitive to multimodality, potentially converging to local modes and yielding non-identifiable solutions. Bayesian MCMC-based inference or MAP estimation faces similar computational and convergence issues in the same regimes [17; 6].

More fundamentally, NLME-ODEs likelihoods are often highly non-convex because ODE solutions depend nonlinearly on the parameters, producing rugged objective landscapes and compounding identifiability and optimization difficulties. These challenges motivate optimization-based alternatives that can regularize inference and scale to complex models. Variational inference (VI) replaces stochastic sampling with a tractable, parametric approximation of the posterior, optimized by maximizing a lower bound. In particular, variational autoencoders (VAEs) use neural networks to parameterize this approximation, offering a scalable and flexible framework for NLME-ODEs inference.

This coupling between NLME-ODEs and neural networks have already been investigated for other purposes. [25] proposed to integrate expert pharmacological ODEs into Neural ODEs to provide additional insights to the clinicians in disease-progression modeling. [5] introduced low-dimensional neural ODEs implemented directly in Monolix [30], showing how inter-individual variability can be incorporated into neural ODEs parameters while preserving pharmacometric interpretability. [10] compared FO/FOCE and VI approaches for mixed-effect estimation in Deep Compartment Models. [1] developed an amortized inference framework in which a conditional normalizing flow is trained on simulated data generated from the parametric ODE model to approximate the parameter posterior distribution. [19] introduced a Deep Nonlinear Mixed-Effect [23] modeling framework in which a neural network term is embedded directly into the ODE to represent unknown system dynamics. [28] performed amortised Bayesian inference for NLME-ODEs through a conditional VAE with population-, group-, and individual-level latent blocks.
Despite these advances, most deep or variational extensions of NLME-ODEs emphasize predictive performance and scalability over statistical interpretability and identifiability. Amortized approaches [28; 1] provide efficient posterior approximations but do not explicitly address identifiability. Conversely, methods that embed neural networks into pharmacometric ODEs [5; 19] typically rely on SAEM or FOCE, and thus remain computationally demanding and sensitive to initialization; moreover, neural components is generally nonidentifiable due to redundancy. Recent ELBO-based inference methods for NLME-ODE models close to ours have been proposed by [29]. However, they do not address practical identifiability and uncertainty quantification, as they do not construct a variance estimator. We focus on these aspects because they are essential for interpretability. Moreover, they do not provide a rigorous comparison with SAEM-based methods on these points.

In this paper, we propose a interpretable, identifiability-aware amortized variational inference framework for NLME-ODEs tailored to sparse and irregular longitudinal settings. We maximize the evidence lower bound (ELBO), yielding a tractable objective that avoids the sampling step of SAEM. For each subject, the posterior distribution of random effects is approximated within a parametric family parameterized by a lightweight, shallow neural encoder to prevent overfitting and non-identifiability given the limited data; the decoder is a parametric ODE model which reconstructs the trajectory from the random effects, ensuring that the inferred representations retain a mechanistic interpretation as subject-specific biological parameters. Although the method is compatible with higher-dimensional random-effects structures, in this work we primarily evaluate its performance in data-scarce regimes and leave a systematic study of scalability with respect to the random-effect dimension to future work.

Beyond point estimation, we quantify population-level parameter uncertainty via a vari-
ance estimator based on the observed Fisher information matrix, computed by automatic differentiation at the variational optimum for population parameters.

Identifiability is fundamental for interpretable inference in NLME models. When neural network components are introduced, this property can be compromised by additional nuisance parameters or redundant representations. We further assess convergence behavior to verify the practical identifiability of the underlying ODEs, ensuring that each population parameter can be uniquely determined from the data. These steps guarantee that the resulting inference remains stable and interpretable within the NLME-ODEs framework.

The paper is organized as follows: **Section 2** introduces the proposed method, including the variational inference strategy, identifiability analysis, and uncertainty quantification procedure; **Section 3** reports simulation results in which we compare our method with the approach SAEM implemented in Monolix [30] in terms of estimation accuracy and identifiability issues under various scenarios, in particular for the estimation of NLME-ODEs of increasing complexities in terms of numbers of estimated parameters; **Section 4** presents an application to real clinical data on antibody concentration dynamics following immunization with a COVID-19 vaccine as in [7], demonstrating the practical relevance and robustness of the proposed framework.

## 2. Amortized Variational Inference for NLME-ODEs

### 2.1 Marginal Likelihood and Variational Reformulation

Classical inference for NLME models maximizes the population marginal log-likelihood

$$ \ell(\phi; \mathbf{Y}) = \sum_{i=1}^{n} \log p_{\phi}(Y_i), \quad p_{\phi}(Y_i) = \int p_{\phi}(Y_i | b_i) p_{\phi}(b_i) db_i, \quad (2) $$

where $\mathbf{Y} = \{Y_i\}_{i=1}^n$ and $Y_i = \{Y_{ij}\}_{j=1}^{n_i}$ are observations at times $\{t_{ij}\}_{j=1}^{n_i}$. Under the observation model 1:

$$ p_{\phi}(Y_i | b_i) = \prod_{j=1}^{n_i} \mathcal{N}(Y_{ij} - h_{\theta_i}(X_i(t_{ij})), \Sigma), \quad (3) $$
where $\mathcal{N}(\cdot, \Sigma)$ denotes the density of a centered Gaussian with covariance $\Sigma$ evaluated at the residual $Y_{ij} - h_{\theta_i}(X_i(t_{ij}))$. The integral over $b_i$ is generally intractable, so we adopt variational inference and introduce $q_{\psi}(b_i | Y_i) = \mathcal{N}(\mu_{\psi}(Y_i), L_{\psi}(Y_i)L_{\psi}(Y_i)^{\top})$, where $L_{\psi}(Y_i)$ is a lower-triangular Cholesky factor (with positive diagonal). This approximation to $p_{\phi}(b_i | Y_i)$ yields the individual ELBO following [13]:

$$ \mathcal{L}_{\text{elbo}_i}(\phi, \psi; Y_i) = \mathbb{E}_{q_{\psi}(b_i|Y_i)}[\log p_{\phi}(Y_i | b_i)] - \text{KL}(q_{\psi}(b_i | Y_i) \| p_{\phi}(b_i)), \quad (4) $$

and the population objective is

$$ \mathcal{L}_{\text{ELBO}}(\phi, \psi; \mathbf{Y}) = \sum_{i=1}^{n} \mathcal{L}_{\text{elbo}_i}(\phi, \psi; Y_i) \le \ell(\phi; \mathbf{Y}), \quad (5) $$

where the inequality follows from Jensen's inequality (see Supplementary Web Materials Section 1 for details).

In our model, we have two distinct sets of parameters:

(1) structural parameters of interest: $\phi = (\bar{\theta}, L_{\Omega}, L_{\Sigma}) \in \mathbb{R}^{D_{\phi}}$.

(2) nuisance parameters: the variational parameter $\psi \in \mathbb{R}^{D_{\psi}}$ that parameterizes the approximate posterior $q_{\psi}$.

So the estimator can be defined as:

$$ (\hat{\phi}, \hat{\psi}) = \arg \max_{(\phi, \psi)} \mathcal{L}_{\text{ELBO}}(\phi, \psi; \mathbf{Y}) \quad (6) $$

## 2.2 Variational Autoencoder

To optimize the ELBO, we adopt amortized variational inference as in [13], where a shared neural network encoder parameterized by $\psi$ maps each subject's longitudinal observations to the parameters of the approximate posterior. The expressiveness of the variational family can be adapted to data availability and model complexity. In this work we use a Gaussian family for computational efficiency, but richer distributions (e.g., Gaussian mixtures or normalizing flows) can be employed when greater flexibility is required. The encoder outputs $\mu_{\psi}(Y_i)$ and $L_{\psi}(Y_i)$ of the approximate posterior distribution. Coupled with an ODEs-based decoder
that reconstructs the data from the random effects, this architecture forms a VAE adapted
to NLME-ODEs.

Amortized inference offers several advantages in the context of longitudinal data with sparse available measurement per subject, such as in clinical trial setting. By sharing statistical strength across subjects, it improves inference stability even when individual trajectories are short, sparse, or irregularly sampled. Once trained, the encoder provides near-instant posterior estimates for new subjects, making the approach computationally efficient and scalable to large cohorts. Furthermore, by learning a smooth mapping from observed data to latent parameters, the amortized encoder regularizes the inverse problem, reducing overfitting and enhancing generalization to unseen clinical profiles.

The encoder architecture depends on the sampling design of the longitudinal data. For regularly sampled trajectories with fixed sequence length, we use a lightweight encoder combining one-dimensional convolutional and projection layers with smooth nonlinear activation functions such as the Gaussian Error Linear Unit (GELU). Given the regular grid, local temporal dependencies can be captured effectively by convolutions, and compared with recurrent architectures, this design uses fewer parameters and is therefore less prone to overfitting in our setting, while maintaining numerical stability and identifiability of the model parameters. For irregularly sampled trajectories, we instead use a recurrent encoder, which takes as input both the observed values and the elapsed time between successive measurements ($\Delta t_{ij} = t_{ij} - t_{i,j-1}$), together with a normalized absolute time covariate. These additional time inputs allow the network to adapt its dynamics to heterogeneous sampling intervals, which is crucial when the timing of observations carries information about the underlying process. Given the limited size of the data in our application, we use shallow architectures with small hidden dimensions smaller than 32 to limit overfitting and avoid introducing non-identifiability due to excessive network capacity.
In real-world applications, longitudinal sequences typically have variable lengths. To enable efficient batching, we pad each sequence to a common length *T* (the maximum sequence length in the dataset) and define a binary mask $m_{it} \in \{0, 1\}$ indicating whether time step *t* for subject *i* is observed ($m_{it} = 1$) or corresponds to padding/missingness ($m_{it} = 0$). The masked attention pooling is applied by aggregating variable-length sequences into a fixed-dimensional subject representation throughout the encoder. In the reconstruction term, the mask ensures that padded positions do not contribute to the objective [31].

Given a latent sample $b_i$, the decoder reconstructs the corresponding subject-specific dynamics by integrating the nonlinear ODEs system $X_i(t) = \text{ODESolve}(f_{\theta_i}, x_{i,0}, t)$. We solve the mechanistic ODE systems using adaptive time-stepping methods, selecting the solver according to the stiffness properties of the dynamics. For non-stiff or mildly stiff systems, we use an explicit Runge-Kutta method of order five. For models exhibiting moderate stiffness, we use a fifth-order diagonally implicit Runge-Kutta scheme, which provides improved stability in the presence of fast transient modes. Relative and absolute error tolerances are set adaptively and specified per model in the corresponding experimental section. All integrations are performed with the `Diffrax` [11] library.

During the training, the parameters are updated using Adam [12]. More details of the optimization, in particular the stopping criterion rules, are in the Supplementary Web Materials Section 4.

### 2.2.1 Pitfall of posterior collapse.
A well-known failure mode in VAE is posterior collapse [8], where the encoder learns a distribution too close to the prior which is fixed as $N(0, I_{D_b})$ and fails to capture individual-level variability. In our case, the risk of posterior collapse is inherently reduced due to two modeling choices:

* Learned variance of the prior $\Omega$: Instead of fixing the prior, we jointly learn its covariance
structure, which increases flexibility and prevents the prior from dominating the variational posterior;

* ODEs decoder conditioned on $b_i$: The subject-specific random effects $b_i$ directly influences the ODEs vector field, ensuring it has a non-negligible effect on the data likelihood.

## 2.3 Parameters Identifiability Analysis

Parameter identifiability is a fundamental requirement in mixed-effects and dynamical system modeling, as it determines whether model parameters can be uniquely inferred from observed data. In the context of inverse problems, accurate and interpretable parameter estimation depends critically on this property: a model is identifiable if distinct parameter values yield distinct distributions of the data. Formally, identifiability is defined as

$$ \forall \phi_1, \phi_2 \in \Phi, \quad p_{\phi_1}(\mathbf{Y}) = p_{\phi_2}(\mathbf{Y}) \Rightarrow \phi_1 = \phi_2. $$

This ensures that the mapping from parameters $\phi$ to the data distribution $p_{\phi}(\mathbf{Y})$ is injective.

We distinguish structural identifiability, which concerns the uniqueness of parameters given perfect, noise-free observations of the system's states [3; 2], from practical identifiability, which reflects the ability to recover parameters from finite and noisy data [27; 15].

In this work, we assume that the structural identifiability of the ODEs is established, that is, the underlying biological dynamics are parameterized in a way that allows unique recovery of the true parameters under ideal conditions. We assess practical identifiability empirically by repeatedly reinitializing the estimation procedure from different random parameter values, also called convergence assessement. Consistent convergence of the learned parameters across runs suggests practical identifiability and indicates that the inference algorithm avoids the local optima that often arise in classical likelihood-based methods and signal non-identifiability.
## 2.4 Uncertainty Quantification of Population-Level Parameters

By assuming that the approximate variational posterior $q_{\hat{\psi}}(b_i|Y_i)$ provides a sufficiently accurate approximation of the true posterior $p_{\phi}$, we estimate the variance of $\hat{\phi}$ using the inverse observed Fisher information derived from the marginal likelihood.

Following [18], the individual ELBO error decomposes as:

$$ \log p_{\phi}(Y_i) - \mathcal{L}_{\text{elbo},i}(\phi, \psi; Y_i) = \underbrace{[\log p_{\phi}(Y_i) - \mathcal{L}(\phi, q^*; Y_i)]}_{\text{variational gap}} + \underbrace{[\mathcal{L}(\phi, q^*; Y_i) - \mathcal{L}_{\text{elbo},i}(\phi, \psi; Y_i)]}_{\text{amortization gap}} $$

where $q^*(b_i | Y_i) \in \mathcal{Q}$ maximizes the ELBO in the chosen variational family. The amortization gap arises because a single encoder $q_{\psi}(b_i | Y_i)$ amortizes inference across subjects rather than optimizing a separate $q_i$ for each subject $i$. [18] shows that this gap can vanish when an ideal inference function exists and the encoder class is sufficiently expressive for exchangeable latent-variable models with $p_{\phi}(b, \mathbf{Y}) = \prod_{i=1}^{n} p_{\phi}(b_i) p_{\phi}(Y_i | b_i)$. Our NLME-ODE model satisfies this factorization since $b_i \sim \mathcal{N}(0, \Omega)$ are i.i.d. given $\phi$. Hence an ideal mapping $f^*(Y_i)$ to the optimal variational parameters exists, and when $q_{\psi}(b_i | Y_i)$ approximates $f^*$ within the Gaussian family, the amortization gap is negligible.

The variational gap equals $KL(q^*(b_i | Y_i) || p_{\phi}(b_i | Y_i))$ and generally does not vanish for a fixed $\mathcal{Q}$ due to variational misspecification. In our framework, we do not place a prior distribution on the population parameters $\phi$; instead, $(\phi, \psi)$ are learned by maximizing the population-level ELBO: $(\hat{\phi}, \hat{\psi}) = \arg\max_{\phi, \psi} \frac{1}{n} \sum_{i=1}^{n} \mathcal{L}_{\text{elbo},i}(\phi, \psi; Y_i)$, so $(\hat{\phi}, \hat{\psi})$ is an M-estimator with $\psi$ a shared nuisance parameter. Equivalently, $\hat{\phi}$ maximizes the profiled criterion $M_n(\phi) = \frac{1}{n} \sum_{i=1}^{n} \mathcal{L}_{\text{elbo},i}(\phi, \hat{\psi}; Y_i)$. If the variational gap (and its $\phi$-derivatives) is uniformly small near the true parameter, maximizing the ELBO is close to maximizing $\frac{1}{n} \sum_{i=1}^{n} \log p_{\phi}(Y_i)$, so that $\hat{\phi}$ behaves like a quasi-MLE.

When both gaps are small, the population ELBO tightly approximates the marginal log-likelihood. In this regime, we quantify uncertainty in $\hat{\phi}$ by approximating $\log p_{\phi}(Y_i)$ via Monte
Carlo over $b_i$ and computing the Hessian at $\hat{\phi}$, yielding the observed information

$$I_n(\hat{\phi}) = - \sum_{i=1}^{n} \nabla_{\hat{\phi}}^2 \log p_{\hat{\phi}}(Y_i), \qquad (7)$$

where $p_{\hat{\phi}}(Y_i)$ denotes the marginal likelihood for subject $i$ at $\hat{\phi}$.

To evaluate $p_{\hat{\phi}}(Y_i)$ and its derivatives, we apply a Monte Carlo approximation based on the reparameterization trick. We express the random effects $b_i$ as a transformation of a noise variable $p_{\epsilon} \sim N(0, I_{D_b})$ independent of $\phi$: $b_i = \mathcal{T}_{\hat{\phi}}(\epsilon) = L_{\hat{\Omega}}\epsilon$ where $L_{\hat{\Omega}}$ is the estimated lower-triangular Cholesky factor of the prior distribution $p_{\hat{\phi}}(b_i)$. So we can approximate the marginal likelihood using monte carlo as:

$$p_{\hat{\phi}}(Y_i) \approx \frac{1}{L} \sum_{\ell=1}^{L} p_{\hat{\phi}}(Y_i | b_i^{(\ell)}) = \frac{1}{L} \sum_{\ell=1}^{L} p_{\hat{\phi}}(Y_i | \mathcal{T}_{\hat{\phi}}(\epsilon^{(\ell)})), \quad p_{\epsilon^{(\ell)}} \sim N(0, I_{D_b}). \qquad (8)$$

The corresponding gradient and Hessian of the marginal likelihood are:

$$\nabla_{\phi} p_{\hat{\phi}}(Y_i) = \mathbb{E}_{p(\epsilon)}[\nabla_{\phi} p_{\hat{\phi}}(Y_i | \mathcal{T}_{\hat{\phi}}(\epsilon))] \approx \frac{1}{L} \sum_{\ell=1}^{L} \nabla_{\phi} p_{\hat{\phi}}(Y_i | \mathcal{T}_{\hat{\phi}}(\epsilon^{(\ell)})),$$

$$\nabla_{\phi}^2 p_{\hat{\phi}}(Y_i) = \mathbb{E}_{p(\epsilon)}[\nabla_{\phi}^2 p_{\hat{\phi}}(Y_i | \mathcal{T}_{\hat{\phi}}(\epsilon))] \approx \frac{1}{L} \sum_{\ell=1}^{L} \nabla_{\phi}^2 p_{\hat{\phi}}(Y_i | \mathcal{T}_{\hat{\phi}}(\epsilon^{(\ell)})).$$

Finally, The Hessian of the log-marginal likelihood is computed as:

$$\nabla_{\phi}^2 \log p_{\hat{\phi}}(Y_i) = \frac{\nabla_{\phi}^2 p_{\hat{\phi}}(Y_i) \cdot p_{\hat{\phi}}(Y_i) - \nabla_{\phi} p_{\hat{\phi}}(Y_i) \nabla_{\phi} p_{\hat{\phi}}(Y_i)^{\top}}{p_{\hat{\phi}}(Y_i)^2}. \qquad (9)$$

## 2.5 Subject-specific parameter estimators

In NLME models, after estimating the population parameters, likelihood-based methods provide subject-specific parameter estimates known as Empirical Bayes estimates, defined as the maximizers of the estimated a posteriori distributions $\hat{b}_i^{\text{EBE}} = \arg \max_{b_i} p_{\hat{\phi}}(b_i | Y_i)$. In our framework, the trained encoder provides the variational approximation $q_{\hat{\psi}}(b_i | Y_i)$ and yields EBEs in a single forward pass. Since $q_{\hat{\psi}}(b_i | Y_i)$ is Gaussian with mean $\mu_{\hat{\psi}}(Y_i)$, we use $\hat{b}_i^{\text{VAE}} := \arg \max_{b_i} q_{\hat{\psi}}(b_i | Y_i) = \mu_{\hat{\psi}}(Y_i)$, which provides scalable subject-specific estimates without per-subject optimization.
### 3. Simulations

To validate our approach, we evaluate its performance and compare it with a classic likelihood/SAEM based approach on several benchmarks constituted of synthetic data simulated from NLME-ODEs with known ground-truth parameters. In the simulation study, we consider the special case of independent random effects, i.e., $\Omega = \text{diag}(\omega_1^2, \dots, \omega_{d_h}^2)$, and i.i.d. observation noise with diagonal covariance $\Sigma = \text{diag}(\sigma_1^2, \dots, \sigma_{d_{n_i}}^2)$. We assess the accuracy of pointwise estimation via a Monte-Carlo approach in which $N_{MC} = 100$ trials of simulated data are constituted from which parameters are estimated. While larger Monte Carlo sample sizes (e.g., $N_{MC} = 1000$) are often recommended, we found $N_{MC} = 100$ to provide stable estimates given the substantial computational cost of each fit. From the estimated parameters $\hat{\phi}_{VAE_k}$ and $\hat{\phi}_{LL_k}$ where $k = 1, \dots, N_{MC}$ obtained from our approach and Monolix (based on SAEM) respectively, we compare pointwise estimation accuracy using three standard metrics:

* Relative Root Mean Squared Error: RRMSE = $\frac{1}{\phi}\sqrt{\frac{1}{N_{MC}}\sum_{k=1}^{N_{MC}}(\hat{\phi}_k - \phi)^2}$;

* Relative Bias: $\text{Rel.Bias} = \frac{1}{\phi}\left(\frac{1}{N_{MC}}\sum_{k=1}^{N_{MC}}\hat{\phi}_k - \phi\right)$;

* Empirical Variance: $\text{Emp.Var.} = \frac{1}{N_{MC}}\sum_{k=1}^{N_{MC}}\hat{\phi}_k^2 - \left(\frac{1}{N_{MC}}\sum_{k=1}^{N_{MC}}\hat{\phi}_k\right)^2$.

We compute the estimated variance (Est. Var) for our approach described in Section 2.4 for each obtained $\hat{\phi}_{VAE_k}$ and the classic Fisher based one available for $\hat{\phi}_{LL_k}$. The accuracy of the variance estimation is then assessed by comparison of the theoretical variance approximated this way with the empirical one described before which are supposed to estimate the same quantity. We also evaluate the frequentist coverage rate of the 95% confidence interval derived from the Est. Var and Emp. Var (Est. Cov and Emp. Cov). We systematically evaluate practical identifiability for each inference procedure.

All experiments were run on a single machine equipped with an NVIDIA A100 GPU (40 GB), using float64 precision. Importantly, the proposed VAE-based NLME-ODE approach is not GPU-intensive: we were able to run the full pipeline on a standard laptop (Intel Core
i7 CPU with an NVIDIA RTX 2000 GPU) including the variance estimation step, though it is slower than on the A100. For the SAEM benchmark, we used Monolix 2023R1 on a high-performance CPU server equipped with 2x 16-core Skylake Intel Xeon Gold 6142 @ 2.6 GHz and 384 GB of RAM. The implementation details are in the Supplementary Web Materials Section 5. Due to space constraints, the details and results for the pharmacokinetics simulation study are reported in Supplementary Web Materials Section 5.1.

## 3.1 Antibody kinetics

*Simulation setting.* We designed a simulation study to assess the performance of our VAE-based inference framework: a partially observed NLME-ODEs describing antibody kinetics following multiple vaccine injections developed in [7]. The mechanistic model for the given $i$-th subject is defined as:

$$ 
\left\{
\begin{aligned}
\dot{S}_i(t) &= \bar{f}_{M_{k,i}} e^{-\delta_V(t-t_k)} - \delta_S S_i(t), \\
\dot{Ab}_i(t) &= \vartheta_i S_i(t) - \delta_{Ab} Ab_i(t), \\
(S_i(0), Ab_i(0)) &= (0.01, 0.1)
\end{aligned}
\right.
\qquad (10) $$

with three injection times $t_{k=1,2,3}$. The above system represents a two-compartment mechanistic model capturing the post-vaccination antibody response dynamics. The state variable $S_i$ denotes the secreting cells population responsible for antibody production, while $Ab_i$ denotes blood circulating antibody concentration. The parameters $\bar{f}_{M_k}$ represents the fold-change of B-cells population magnitude able to differentiate into secreting cells after $k$-th injection compared to the first one (by definition $\bar{f}_{M_{1,i}} = 1$). The decay rate $\delta_S$ controls the lifespan of the secreting cell population, and $\delta_{Ab}$ governs the antibody degradation rate. $\delta_V$ is the degradation rate of the vaccine induced antigen. The parameter $\vartheta_i$ represents the subject-specific antibody production rate per secreting cells, capturing inter-individual variability in immune response intensity. The pair $(\delta_S, \delta_{Ab})$ is not structurally identifiable when relying solely on $Ab_i$ observations. To break the symmetry that causes structural non-identifiability,
we impose the constraint: $log(\delta_{Ab}) = log(\delta_S + \exp(\lambda))$, where $\lambda \in \mathbb{R}$ is a log-gap parameter enforcing $\delta_{Ab} > \delta_S$ in the ODEs Eq. 10 (more details in Supplementary Web Materials Section 3).

Inter-individual variability is introduced in the parameters $\vartheta_i$ and $\bar{f}_{M_{2,i}}$ according to

$$ \log(\vartheta_i) = \log(\vartheta) + b_{\vartheta,i}, \quad b_{\vartheta,i} \sim \mathcal{N}(0, \omega_{\vartheta}^2), $$

$$ \log(\bar{f}_{M_{2,i}}) = \log(\bar{f}_{M_2}) + b_{\bar{f}_{M_{2,i}}}, \quad b_{\bar{f}_{M_{2,i}}} \sim \mathcal{N}(0, \omega_{\bar{f}_{M_{2,i}}}^2), $$

where $\vartheta$ and $\bar{f}_{M_2}$ are shared across the population. In this system, only the antibody compartment is observed, whereas the upstream driver of the dynamics remains unmeasured. The observed antibody concentration is related to the latent state through a log-transformed observation model with $\sigma_{Ab}: Y_{Ab,ij} = \log_{10} Ab_{ij}(t_{ij}) + \epsilon_{Ab,ij}$, where $\epsilon_{Ab,ij} \sim \mathcal{N}(0, \sigma_{Ab}^2)$. We generated two sets of datasets: 100 datasets with 15 measurements per subject collected at equally spaced time points, and 100 datasets with 10 measurements per subject collected at irregular time points over the 400-day follow-up as illustrated in Figure 1(a) and (b). Each dataset consists of 50 subjects. The true parameter values are in the Supplementary Web Materials Table 1.

**Estimation scenario.** The parameters, such as $\delta_S$ and $\bar{f}_{M_k}$ affect the observations only indirectly. Moreover, compensatory effects between $\delta_{Ab}$ and $\delta_S$ can induce strong correlations in their estimates, further impairing identifiability. We therefore consider three inference scenarios of increasing dimensionality, reflecting progressively more complex settings: (S1) baseline scenario in which we estimate the population parameter: $\phi = (\log(\vartheta), \log(\bar{f}_{M_2}), \log(\bar{f}_{M_3}), \omega_{\vartheta}, \omega_{\bar{f}_{M_2}}, \sigma_{Ab})$ (S2) scenario S1 augmented with the estimation of the degradation rate of Secreting cells $\log(\delta_S)$ (S3) Scenario S2 is further extended with the antibody decay estimation, $\log(\delta_{Ab})$.

**Results and identifiability analysis.** In the sparse regular sampling setting, the VAE approach yields lower or comparable relative bias, RRMSE and variance of the estimations
for the population parameters compared to SAEM, confirming its accuracy in recovering fixed effects in the scenario 1. However, it performs slightly worse for the random-effect variances, where a modest underestimation is observed. This behavior is expected, as variational inference typically produces a tighter posterior approximation that tends to underestimate uncertainty components [4]. In Scenario 2, SAEM exhibits a clear inconsistency between its estimated variances and the empirical variances computed across the $N_{MC} = 100$ simulated datasets for the mean population parameter $\log(\vartheta)$, $\log(\tilde{f}_{M_2})$ and $\log(\tilde{f}_{M_3})$. The disagreement can be attributed to the FIM computation in SAEM, which relies on local curvature at the converged point. In this nonlinear setting with correlated parameters, the MCMC-based E-step may mix slowly and drive SAEM toward a suboptimal local basin, where the likelihood appears spuriously sharp in certain directions, resulting in variance underestimation.

In the more complex scenario 3, SAEM frequently converges to local minima and fails to recover several parameters as shown in Figure 2(A), preventing replication over $N_{MC} = 100$ datasets. The VAE method, by contrast, remains numerically stable and produces coherent parameter estimates, demonstrating its robustness in moderately ill-conditioned inference regimes.

In the sparse irregular sampling setting, the SAEM algorithm exhibits a loss of practical identifiability across all three scenarios. This behavior becomes more pronounced as the scenarios become more challenging, and is particularly visible for parameters governing the decay dynamics (e.g. $\delta_S$, $\delta_{Ab}$), for which SAEM frequently converges away from the true value. By contrast, the VAE-based method remains stable under the same sparse irregular design. We report only the results for Scenario 3 obtained with the VAE in the Table 1. For SAEM, a non-negligible fraction of simulation runs did not yield reliable standard error estimates because the observed Hessian was singular for key parameters. In particular, in Scenario 2 and Scenario 3, the Hessian blocks associated with $\delta_S$ and $\delta_{Ab}$ were not invertible
for approximately 15% and 20% of the datasets respectively, preventing computation of the estimated covariance matrix and downstream uncertainty summaries.

Overall, these results confirm that partial observation combined with sparse sampling amplifies the practical identifiability challenges of NLME-ODEs. In such settings, SAEM's reliance on finite-sample MCMC approximations in complex posterior spaces can lead to unstable or biased estimation. The VAE approach, by directly optimising a smoother ELBO, avoids the sampling bottleneck and better preserves identifiability as the parameter space grows.

## 3.2 Dynamics of TGF Activation in Asthmatic Airways

To further challenge the proposed inference methods, we next consider a more complex ODE system with richer nonlinear interactions. To reproduce the dynamics of airway remodeling, we adopt the reduced system proposed by [24], which captures the fast biochemical activation of transforming growth factor $\beta$ (TGF-$\beta$) mediated by airway smooth muscle (ASM) and extracellular matrix (ECM) interactions. The model describes the temporal evolution of the densities of proliferating ASM cells $p(t)$, contractile ASM cells $c(t)$, active TGF-$\beta$ concentration $a(t)$, and ECM $m(t)$ as follows:

$$ 
\left\{
\begin{aligned}
\frac{dp_i}{dt} &= \kappa_{p_i} p_i \left(1 - \frac{p_i}{p_{\max}}\right) \left(1 + \frac{\kappa_{ap} a_i p_i}{\bar{\eta}_{ap} + a_i p_i}\right) + \kappa_{cp} \frac{\gamma_c}{\gamma_p} c_i - \kappa_{pc} p_i, \\
\frac{dc_i}{dt} &= \kappa_{pc} \frac{\gamma_p}{\gamma_c} p_i - (\kappa_{cp} + \phi_c) c_i, \\
\frac{da_i}{dt} &= \left( \frac{\kappa_s}{\nu} e^{-\frac{(t-\tau_i)^2}{\nu^2}} + \frac{\kappa_{ac} a_i c_i}{\bar{\eta}_{ac} + a_i c_i} \right) \frac{c_i}{m_i} \gamma_a - \kappa_b \left( \frac{\gamma_p}{\gamma_c} p_i + c_i \right) a_i - a_i, \\
\frac{dm_i}{dt} &= \left( \kappa_{pm} + \frac{\kappa_{apm} a_i p_i}{\bar{\eta}_{apm} + a_i p_i} \right) \gamma_p p_i - \phi_m m_i \\
(p_i(0), c_i(0), a_i(0), m_i(0)) &= (0.1, 0.8, 0.01, 0.9).
\end{aligned}
\right.
\quad (11) $$

The design was implemented with one random effect on the kinetic parameter $k_p$ as:

$$ \log(k_{p_i}) = \log(k_p) + b_i \quad \text{and} \quad b_i \sim \mathcal{N}(0, \omega_{k_p}^2), $$
representing the individual rate of ASM proliferation. Only TGF-beta is observed according to the observation model: $Y_{ij} = \log_{10}(a(t_{ij})) + \epsilon_{a,ij}$ where $\epsilon_{a,ij} \sim N(0, \sigma_a^2)$ with 20 longitudinal measurements uniformly distributed over 400 days, example of randomly drawn data are represented in Figure 1(c). Each dataset consists of 50 subjects. The true parameters are in the Supplementary Web Materials Table 3.

*Estimation scenarios.* The inference task is made difficult by the strongly coupled, feedback-rich ODE structure (logistic proliferation, phenotype switching, Hill-type saturations, and multiplicative terms such as $a_ic_i$), which induces nonconvex likelihoods, strong parameter interactions, and practical near-nonidentifiabilities. We consider three inference scenarios of increasing dimensionality: (S1) a baseline setting estimating $\phi = \{\log(k_p), \omega_{k_p}, \sigma_a\}$; (S2) S1 augmented by estimating $\log(k_b)$; and (S3) S2 further extended by estimating $\log(k_{ac})$.

*Results.* As shown in Table 2, the VAE approach consistently outperforms SAEM, with smaller relative bias and lower RRMSE across all parameters. Regarding practical identifiability (Figure 2D), parameters appear identifiable in Scenarios 1–2, but SAEM becomes unstable in Scenario 2: FIM-based coverage for $k_b$ drops to 2% due to severe variance underestimation. This likely reflects a challenging nonconvex likelihood, causing SAEM to converge to suboptimal regions with spuriously high curvature in the $k_b$ direction. As a result, empirical variance is much larger (empirical coverage 86%), while VAE provides more stable estimates and uncertainty, especially for $k_b$.

The scenario 3 leads to more complex setting: adding $k_{ac}$ strengthens the positive feedback $a↑\Rightarrow c↑$, $m↑\Rightarrow a↑$, increasing effective stiffness and deepening nonconvexity. VAE remains stable with acceptable bias and coverage, whereas SAEM loses practical identifiability, consistent with well-known challenges of SAEM on stiff, highly nonlinear NLME-ODEs when curvature is poorly conditioned. We do not execute SAEM over the 100 simulated datasets in this scenario.
Overall, these experiments confirm that the proposed VAE approach can handle the nonlinear, partially stiff dynamics of the TGF-$\beta$ activation with reasonable accuracy and calibrated uncertainty, even when SAEM becomes numerically unstable.

# 4. Application on real data

In the real data case, we relied on data from a cohort of Bnt162b2 vaccine recipients, in which both antibody kinetics and neutralizing activity were measured longitudinally against SARS-CoV-2. In brief, SARS-CoV-2 naive patients were recruited in Orléans, France between August 27, 2020 and May 24, 2022. Individuals were followed for up to 483 days after their first vaccine injections (see more details on the data in [21]). The data pre-treatment and selection process is exactly the same as in [7]. We used only the antibody kinetics data in the original dataset and analysed $N = 25$ subjects who received three injections. The evolution of antibody titers (BAU/ml) is shown in Figure 1(d) over a 500-day period. For the analysis of this dataset, we employed exactly the same model design as the simulation study in the section 3.1. We estimate the following parameter set: $\phi = \{\log(\vartheta), \log(\bar{f}_{M_2}), \log(\bar{f}_{M_3}), \log(\delta_S), \log(\delta_{Ab}), \omega_\vartheta, \omega_{\bar{f}_{M_2}}, \sigma_\epsilon\}$,

## 4.1 Results

Table 3 reports the estimates on the real dataset obtained with both methods. Overall, the two methods provide consistent estimates for the population parameters and observation noise, with all corresponding 95% confidence intervals overlapping. The main difference appears in the random-effect variance associated with $\bar{f}_{M_2}$, where the VAE estimate is smaller than that obtained by SAEM. Nevertheless, the confidence intervals still overlap, indicating no statistically significant discrepancy. A plausible explanation is limited identifiability of the random-effect variance for $\bar{f}_{M_2}$ in the real dataset. In such settings, different inference schemes may yield similar population means while producing more variable estimates of vari-
ance components. The smaller VAE estimate is consistent with the well-known tendency of variational inference to underestimate posterior uncertainty, but the overlapping confidence intervals suggest that the data do not strongly distinguish between the two estimates. As shown in Figure 2(C), the VAE and SAEM yielded consistent parameter estimates regardless of starting points.

At the population level without accounting for interindividual random effects, Figure 3 compares the estimated population mean trajectories obtained by the VAE and SAEM approaches. For the antibody kinetic, the VAE-based trajectory provides a more accurate reconstruction: the SAEM trajectory consistently falls within the 95% prediction interval of the VAE approach. This indicates that the variational approach captures the underlying uncertainty more faithfully, while SAEM remains well-aligned but more dispersed. For the secreting cells, however, a small discrepancy appears in the peak region after 3rd injection, where SAEM deviates outside the VAE related confidence interval. Figure 3 also illustrates individual predictive trajectories under the VAE-based NLME-ODEs, generated using the mode of the variational posterior for the subject-specific random effects $b_i$. The predicted means match the observed data, and the 95% prediction intervals show good overall coverage, confirming that the variational posterior captures the main sources of individual uncertainty.

# 5. Discussion

In this work, we use a variational autoencoder for parameter inference in NLME-ODE models, replacing the sampling step in classical methods such as SAEM. The random-effects posterior is approximated by a Gaussian variational distribution parameterized by a lightweight amortized encoder, designed for sparse biological and clinical datasets. Through-out the design, particular attention was devoted to practical identifiability to ensure that the introduction of variational parameters does not compromise the interpretability of the model. For uncertainty quantification, we evaluate the observed FIM via direct hessian of Monte
Carlo approximated marginal likelihood using prior-based reparametrization and automatic differentiation. In terms of estimation accuracy, for fixed-effect parameters, our variational method achieved performance comparable to SAEM on real data and outperformed it in more challenging simulated scenarios involving stronger nonlinearities, partial observation. As the number of population parameters increases, the performance gap widened in favor of the variational approach, illustrating its robustness to the increased curvature and multimodality of the likelihood surface—conditions under which SAEM often suffers from slow convergence or attraction to suboptimal local maxima. Both the empirical coverage and the coverage predicted from the estimated variances remained generally close to the nominal 95% level, especially for the fixed effect parameters.

Future extensions could focus on increasing the flexibility and robustness of the variational approximation. One direction is to enrich the variational posterior, through for example normalizing flows, to alleviate underestimation of random-effect variances and better capture posterior correlations. Another avenue is to move toward a fully variational Bayes formulation by placing priors on the population-level parameters. Beyond this, we will apply the framework to higher-dimensional mechanistic models with covariates, censored data, or large sets of random effects. Extending the methodology to semi-parametric, for instance via neural augmentations of ODE components, also raises new identifiability challenges whose interaction with variational regularization deserves careful study. Finally, a theoretical analysis of the estimator, encompassing conditions for consistency and asymptotic normality of the ELBO estimator would strengthen the statistical foundations of the proposed approach.

## Acknowledgments

This work receives financial support from Inserm to the booster program Exposome, also has been supported by SOLVE, funded by the European Union's Horizon Europe Research and Innovation programme under grant n°101137185. We are deeply thankful to Delphine Planas,
Timothée Bruel, Laurent Hocqueloux, Thierry Prazuck and Olivier Schwartz for the production and the public availability of the data used in section 4. The experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr). We thank Simulations Plus, Lixoft division for the free academic use of the MonolixSuite.

## SUPPLEMENTARY MATERIAL

Tables and Figures for Sections 3 and 4 are available with this paper.

Code and data are available at https://github.com/AuroraFr/NN_NLME/tree/paper

## REFERENCES

J. Arruda, Y. Schälte, C. Peiter, O. Teplytska, U. Jaehde, and J. Hasenauer. An amortized approach to non-linear mixed-effects modeling based on neural posterior estimation. In *Proceedings of the 41st International Conference on Machine Learning*, volume 235, pages 1865–1901, 2024.

S. Audoly, G. Bellu, L. D'Angio, M. P. Saccomani, and C. Cobelli. Global identifiability of nonlinear models of biological systems. *IEEE Transactions on Biomedical Engineering*, 48(1):55–65, 2001.

R. Bellman and K. J. Åström. On structural identifiability. *Mathematical Biosciences*, 7: 329–339, 1970.

D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. *Journal of the American Statistical Association*, 112(518):859–877, 2017.

D. S. Bräm, B. Steiert, M. Pfister, B. Steffens, and G. Koch. Low-dimensional neural ordinary differential equations accounting for inter-individual variability implemented in Monolix and NONMEM. *CPT: Pharmacometrics & Systems Pharmacology*, 14:5–16, 2025.
B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell. Stan: A probabilistic programming language. *Journal of Statistical Software*, 76:1–32, 2017.

Q. Clairon, M. Prague, D. Planas, T. Bruel, L. Hocqueloux, T. Prazuck, O. Schwartz, R. Thiébaut, and J. Guedj. Modeling the kinetics of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2. *PLoS Computational Biology*, 19(8):e1011282, 2023.

B. Dai, Z. Wang, and D. Wipf. The usual suspects? reassessing blame for VAE posterior collapse. In *Proceedings of the 37th International Conference on Machine Learning*, pages 2313–2322, 2020.

B. Delyon, M. Lavielle, and E. Moulines. Convergence of a stochastic approximation version of the EM algorithm. *The Annals of Statistics*, 27:94–128, 1999.

A. Janssen, F. C. Bennis, M. H. Cnossen, and R. A. A. Mathôt. Mixed effect estimation in deep compartment models: Variational methods outperform first-order approximations. *Journal of Pharmacokinetics and Pharmacodynamics*, 51:797–808, 2024.

P. Kidger. *On Neural Differential Equations*. PhD thesis, University of Oxford, 2021.

D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In *International Conference on Learning Representations (ICLR)*, 2015.

D. P. Kingma and M. Welling. Auto-encoding variational bayes. In *International Conference on Learning Representations (ICLR)*, 2014.

M. Lavielle. *Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools*. Chapman & Hall/CRC Biostatistics Series. Chapman and Hall/CRC, 2014.

M. Lavielle and L. Aarons. What do we mean by identifiability in mixed effects models? *Journal of Pharmacokinetics and Pharmacodynamics*, 43:111–122, 2016.

M. Lavielle and F. Mentré. Estimation of population pharmacokinetic parameters of
saquinavir in HIV patients with the Monolix software. *Journal of Pharmacokinetics and Pharmacodynamics*, 34:229-249, 2007.

D. J. Lunn, A. Thomas, N. Best, and D. Spiegelhalter. Winbugs—a bayesian modelling framework: concepts, structure, and extensibility. *Statistics and Computing*, 10:325–337, 2000.

C. C. Margossian and D. M. Blei. Amortized variational inference: When and why? In *Proceedings of the 40th Conference on Uncertainty in Artificial Intelligence (UAI)*, volume 244, pages 2434–2449, 2024.

C. J. Martensen, N. Korsbo, V. Ivaturi, and S. Sager. Data-driven discovery of feedback mechanisms in acute myeloid leukaemia: Alternatives to classical models using deep nonlinear mixed effect modeling and symbolic regression. *bioRxiv*, 2024.

J. C. Pinheiro and D. M. Bates. Approximations to the log-likelihood function in the nonlinear mixed-effects model. *Journal of Computational and Graphical Statistics*, 4 (1):12–35, 1995.

D. Planas, D. Veyer, A. Baidaliuk, I. Staropoli, F. Guivel-Benhassine, M. M. Rajah, and et al. Reduced sensitivity of SARS-CoV-2 variant delta to antibody neutralization. *Nature*, 596:276–280, 2021.

M. Prague, D. Commenges, J. Guedj, J. Drylewicz, and R. Thiébaut. Nimrod: A program for inference via a normal approximation of the posterior in models with random effects based on ordinary differential equations. *Computer Methods and Programs in Biomedicine*, 111(2):447–458, 2013.

PumasAI. DeepPumas. Software, 2025. URL https://deeppumas-docs.pumas.ai/.

H. J. Pybus, R. D. O'Dea, and B. S. Brook. A dynamical model of TGF-β activation in asthmatic airways. *Mathematical Medicine and Biology: A Journal of the IMA*, 2023.

Z. Qian, W. R. Zame, L. M. Fleuren, P. Elbers, and M. van der Schaar. Integrating expert
ODEs into neural ODEs: Pharmacology and disease progression. In *Proceedings of the 35th Conference on Neural Information Processing Systems (NeurIPS)*, 2021.

C. Rackauckas, Y. Ma, A. Noack, V. Dixit, P. K. Mogensen, S. Byrne, S. Maddhashiya, J. B. Santiago Calderón, J. Nyberg, J. V. Gobburu, et al. Accelerated predictive healthcare analytics with pumas, a high performance pharmaceutical modeling and simulation platform. *BioRxiv*, pages 2020-11, 2020.

A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, and J. Timmer. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. *Bioinformatics*, 25(15):1923–1929, 2009.

G. Roeder, P. Grant, A. Phillips, N. Dalchau, and E. Meeds. Efficient amortised bayesian inference for hierarchical and nonlinear dynamical systems. In *International Conference on Machine Learning*, pages 4445–4455, 2019.

J. Rohleff, F. Bachmann, U. Nahum, D. Bräm, B. Steffens, M. Pfister, G. Koch, and J. Schropp. Redefining parameter estimation and covariate selection via variational autoencoders: One run is all you need. *CPT: Pharmacometrics & Systems Pharmacology*, 14(12):2232–2243, 2025.

Simulations Plus. Monolix 2024r1, 2024.

A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin. Attention is all you need. In *Advances in Neural Information Processing Systems*, volume 30, 2017.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Randomly selected subject trajectories from the simulated datasets and the real dataset. Points are noisy observations; solid curves show the corresponding noise-free underlying dynamics. Vertical light-grey and black dashed lines mark the second and third vaccine injections (the first occurs at $t=0$).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Convergence assessment of SAEM and VAE across three simulation scenarios (S1–S3) and one real-data analysis. Columns correspond to simulation scenarios (S1, S2, S3). Rows correspond to datasets/cases: (A) antibody kinetics with regular sampling, (B) antibody kinetics with irregular sampling, (C) real antibody dataset, and (D) TGF-β dynamics. For each parameter, points show estimates obtained from multiple random initializations for each method. In simulation panels (A,B,D), the horizontal segment indicates the true parameter value; a tight cluster around it suggests robust convergence, whereas multiple separated clusters indicate competing local optima and sensitivity to initialization. In the real-data panel (C), no true value is available, so convergence is assessed by the concentration of the estimates across initializations.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Real-world dataset. Top: population-level predictions for (a) S-cell and (b) anti-S IgG trajectories, with 95% uncertainty bands for VAE and SAEM. Bottom: subject-specific predictions from the VAE with 95% bands; points denote observed measurements.</figcaption></figure>
<table><caption>Table 1: Comparison of parameter estimates obtained with NN-NLME and SAEM methods on the antibody kinetics simulation datasets in the section 3.1.</caption><thead><tr><th>Scenario</th><th>Parameter</th><th>Method</th><th>Rel. Bias (%)</th><th>RRMSE (%)</th><th>Emp. Var. (10<sup>-2</sup>)</th><th>Est. Var. (10<sup>-2</sup>)</th><th>Emp. Cov.</th><th>Est. Cov.</th></tr></thead><tbody><tr><td colspan="9" style="text-align: center; font-weight: bold;">Antibody kinetics (regular sampling)</td></tr><tr><td rowspan="8">S1</td><td rowspan="2">log(ϑ)</td><td>NN-ELBO</td><td>-0.49</td><td>2.46</td><td>0.59</td><td>0.48</td><td>0.96</td><td>0.91</td></tr><tr><td>SAEM</td><td>-1.08</td><td>2.67</td><td>0.61</td><td>0.61</td><td>0.98</td><td>0.90</td></tr><tr><td rowspan="2">log(f&#773;<sub>M2</sub>)</td><td>VAE</td><td>1.38</td><td>7.12</td><td>1.87</td><td>1.29</td><td>0.94</td><td>0.90</td></tr><tr><td>SAEM</td><td>1.65</td><td>7.22</td><td>1.90</td><td>1.81</td><td>0.95</td><td>0.93</td></tr><tr><td rowspan="2">log(f&#773;<sub>M3</sub>)</td><td>VAE</td><td>0.34</td><td>1.33</td><td>0.14</td><td>0.14</td><td>0.94</td><td>0.94</td></tr><tr><td>SAEM</td><td>1.02</td><td>1.70</td><td>0.16</td><td>0.17</td><td>0.92</td><td>0.77</td></tr><tr><td rowspan="2">ω<sub>ϑ</sub></td><td>VAE</td><td>-3.67</td><td>11.72</td><td>0.31</td><td>0.21</td><td>0.94</td><td>0.82</td></tr><tr><td>SAEM</td><td>1.36</td><td>11.28</td><td>0.31</td><td>0.26</td><td>0.97</td><td>0.93</td></tr><tr><td rowspan="2">ω&#773;<sub>f&#773;<sub>M2</sub></sub></td><td>VAE</td><td>-2.41</td><td>10.58</td><td>0.86</td><td>0.68</td><td>0.94</td><td>0.86</td></tr><tr><td>SAEM</td><td>0.11</td><td>11.23</td><td>1.02</td><td>0.88</td><td>0.94</td><td>0.88</td></tr><tr><td rowspan="2">σ<sub>Ab</sub></td><td>VAE</td><td>1.56</td><td>2.90</td><td>0.06</td><td>0.06</td><td>0.94</td><td>0.94</td></tr><tr><td>SAEM</td><td>1.52</td><td>3.05</td><td>0.07</td><td>0.07</td><td>0.95</td><td>0.94</td></tr><tr><td rowspan="8">S2</td><td rowspan="2">log(ϑ)</td><td>VAE</td><td>-0.96</td><td>2.25</td><td>0.42</td><td>0.45</td><td>0.90</td><td>0.90</td></tr><tr><td>SAEM</td><td>-0.91</td><td>2.55</td><td>0.58</td><td>1.09</td><td>0.98</td><td>0.94</td></tr><tr><td rowspan="2">log(f&#773;<sub>M2</sub>)</td><td>VAE</td><td>1.55</td><td>6.35</td><td>1.46</td><td>1.46</td><td>0.94</td><td>0.90</td></tr><tr><td>SAEM</td><td>1.07</td><td>7.08</td><td>1.90</td><td>2.83</td><td>0.95</td><td>0.95</td></tr><tr><td rowspan="2">log(f&#773;<sub>M3</sub>)</td><td>VAE</td><td>0.60</td><td>1.39</td><td>0.13</td><td>0.14</td><td>0.95</td><td>0.92</td></tr><tr><td>SAEM</td><td>0.78</td><td>1.59</td><td>0.16</td><td>1.15</td><td>0.90</td><td>0.83</td></tr><tr><td rowspan="2">log(δ<sub>s</sub>)</td><td>VAE</td><td>0.12</td><td>0.39</td><td>0.03</td><td>0.03</td><td>0.95</td><td>0.94</td></tr><tr><td>SAEM</td><td>0.04</td><td>0.44</td><td>0.04</td><td>0.06</td><td>0.94</td><td>0.96</td></tr><tr><td rowspan="2">ω<sub>ϑ</sub></td><td>VAE</td><td>-3.01</td><td>11.54</td><td>0.31</td><td>0.20</td><td>0.94</td><td>0.85</td></tr><tr><td>SAEM</td><td>-1.51</td><td>10.38</td><td>0.26</td><td>0.26</td><td>0.95</td><td>0.95</td></tr><tr><td rowspan="2">ω&#773;<sub>f&#773;<sub>M2</sub></sub></td><td>VAE</td><td>-3.19</td><td>10.98</td><td>0.89</td><td>0.85</td><td>0.92</td><td>0.88</td></tr><tr><td>SAEM</td><td>~0.0</td><td>11.32</td><td>0.94</td><td>0.94</td><td>0.94</td><td>0.89</td></tr><tr><td rowspan="2">σ<sub>Ab</sub></td><td>VAE</td><td>1.65</td><td>3.12</td><td>0.07</td><td>0.07</td><td>0.94</td><td>0.94</td></tr><tr><td>SAEM</td><td>1.75</td><td>3.33</td><td>0.08</td><td>0.08</td><td>0.95</td><td>0.94</td></tr><tr><td rowspan="8">S3</td><td rowspan="2">log(ϑ)</td><td>VAE</td><td>-0.70</td><td>2.42</td><td>0.55</td><td>0.73</td><td>0.93</td><td>0.93</td></tr><tr><td>SAEM</td><td>0.62</td><td>6.86</td><td>1.79</td><td>1.46</td><td>0.94</td><td>0.90</td></tr><tr><td rowspan="2">log(f&#773;<sub>M2</sub>)</td><td>VAE</td><td>0.59</td><td>1.38</td><td>0.13</td><td>0.17</td><td>0.95</td><td>0.95</td></tr><tr><td>SAEM</td><td>0.14</td><td>0.52</td><td>0.05</td><td>0.06</td><td>0.93</td><td>0.95</td></tr><tr><td rowspan="2">log(δ<sub>s</sub>)</td><td>VAE</td><td>~0.0</td><td>1.97</td><td>0.25</td><td>0.24</td><td>0.95</td><td>0.95</td></tr><tr><td>SAEM</td><td>-2.81</td><td>11.43</td><td>0.31</td><td>0.25</td><td>0.95</td><td>0.90</td></tr><tr><td rowspan="2">ω<sub>ϑ</sub></td><td>VAE</td><td>-2.96</td><td>10.73</td><td>0.86</td><td>0.72</td><td>0.95</td><td>0.88</td></tr><tr><td>SAEM</td><td>1.2</td><td>2.91</td><td>0.07</td><td>0.07</td><td>0.94</td><td>0.94</td></tr><tr><td colspan="9" style="text-align: center; font-weight: bold;">Antibody kinetics (irregular sampling)</td></tr><tr><td rowspan="8">S3</td><td rowspan="2">log(ϑ)</td><td>VAE</td><td>-0.58</td><td>3.02</td><td>0.89</td><td>0.86</td><td>0.96</td><td>0.95</td></tr><tr><td>SAEM</td><td>-0.2</td><td>6.95</td><td>1.85</td><td>1.92</td><td>0.96</td><td>0.97</td></tr><tr><td rowspan="2">log(f&#773;<sub>M2</sub>)</td><td>VAE</td><td>0.14</td><td>1.63</td><td>0.22</td><td>0.29</td><td>0.95</td><td>0.96</td></tr><tr><td>SAEM</td><td>0.34</td><td>0.52</td><td>0.20</td><td>0.19</td><td>0.98</td><td>0.97</td></tr><tr><td rowspan="2">log(δ<sub>s</sub>)</td><td>VAE</td><td>-1.01</td><td>3.82</td><td>0.95</td><td>0.97</td><td>0.98</td><td>0.97</td></tr><tr><td>SAEM</td><td>-3.65</td><td>12.97</td><td>0.39</td><td>0.34</td><td>0.90</td><td>0.87</td></tr><tr><td rowspan="2">ω<sub>ϑ</sub></td><td>VAE</td><td>-4.75</td><td>14.28</td><td>1.46</td><td>1.04</td><td>0.98</td><td>0.88</td></tr><tr><td>SAEM</td><td>1.68</td><td>3.92</td><td>2.65</td><td>0.69</td><td>0.96</td><td>0.90</td></tr></tbody></table>
<table><caption>Table 2: Comparison of parameter estimates on the simulation datasets from dynamics of TGF-β activation in section 3.2 obtained with the VAE and SAEM methods.</caption><thead><tr><th>Scenario</th><th>Parameter</th><th>Method</th><th>Rel. Bias (1%)</th><th>RRMSE (1%)</th><th>Emp. Var. (10<sup>-2</sup>)</th><th>Est. Var. (10<sup>-2</sup>)</th><th>Emp. Cov.</th><th>Est. Cov.</th></tr></thead><tbody><tr><td rowspan="4">S1</td><td rowspan="2">log(k<sub>p</sub>)</td><td>VAE</td><td>-0.27</td><td>4.89</td><td>0.0047</td><td>0.0061</td><td>0.95</td><td>0.97</td></tr><tr><td>SAEM</td><td>-2.61</td><td>5.47</td><td>0.0045</td><td>0.0052</td><td>0.90</td><td>0.92</td></tr><tr><td rowspan="2">ω<sub>k<sub>p</sub></sub></td><td>VAE</td><td>~ 0.0</td><td>7.68</td><td>0.0015</td><td>0.0040</td><td>0.95</td><td>0.97</td></tr><tr><td>SAEM</td><td>-4.45</td><td>10.60</td><td>0.0023</td><td>0.0030</td><td>0.94</td><td>0.96</td></tr><tr><td rowspan="2">σ<sub>a</sub></td><td>VAE</td><td>1.65</td><td>3.12</td><td>0.07</td><td>0.07</td><td>0.94</td><td>0.94</td></tr><tr><td>SAEM</td><td>1.75</td><td>3.33</td><td>0.08</td><td>0.08</td><td>0.95</td><td>0.94</td></tr><tr><td rowspan="4">S2</td><td rowspan="2">log(k<sub>p</sub>)</td><td>VAE</td><td>0.37</td><td>6.21</td><td>0.010</td><td>0.009</td><td>0.96</td><td>0.95</td></tr><tr><td>SAEM</td><td>0.97</td><td>6.40</td><td>0.008</td><td>0.006</td><td>0.90</td><td>0.90</td></tr><tr><td rowspan="2">k<sub>b</sub></td><td>VAE</td><td>0.03</td><td>1.29</td><td>0.010</td><td>0.019</td><td>0.94</td><td>0.96</td></tr><tr><td>SAEM</td><td>1.40</td><td>1.90</td><td>0.020</td><td>0.0002</td><td>0.90</td><td>0.02</td></tr><tr><td rowspan="2">ω<sub>k<sub>p</sub></sub></td><td>VAE</td><td>-0.07</td><td>9.44</td><td>0.002</td><td>0.004</td><td>0.94</td><td>0.97</td></tr><tr><td>SAEM</td><td>-1.75</td><td>9.99</td><td>0.002</td><td>0.003</td><td>0.93</td><td>0.98</td></tr><tr><td rowspan="2">σ<sub>a</sub></td><td>VAE</td><td>1.55</td><td>3.07</td><td>0.07</td><td>0.07</td><td>0.94</td><td>0.94</td></tr><tr><td>SAEM</td><td>1.45</td><td>3.18</td><td>0.08</td><td>0.08</td><td>0.95</td><td>0.94</td></tr><tr><td rowspan="4">S3</td><td rowspan="2">log(k<sub>p</sub>)</td><td>VAE</td><td>1.12</td><td>6.48</td><td>0.01</td><td>0.009</td><td>0.94</td><td>0.93</td></tr><tr><td>SAEM</td><td>1.45</td><td>2.56</td><td>0.07</td><td>0.1</td><td>0.94</td><td>0.96</td></tr><tr><td rowspan="2">k<sub>b</sub></td><td>VAE</td><td>0.29</td><td>0.50</td><td>0.04</td><td>0.09</td><td>0.89</td><td>0.97</td></tr><tr><td>SAEM</td><td>-1.85</td><td>9.20</td><td>0.002</td><td>0.005</td><td>0.94</td><td>0.95</td></tr><tr><td rowspan="2">σ<sub>a</sub></td><td>VAE</td><td>1.55</td><td>3.07</td><td>0.07</td><td>0.07</td><td>0.94</td><td>0.94</td></tr><tr><td>SAEM</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>Table 3: VAE and SAEM parameter estimates and 95% CIs on real dataset</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="4">VAE</th><th colspan="4">SAEM</th></tr><tr><th>Est.</th><th>95% CI Lower</th><th>95% CI Upper</th><th>Est.</th><th>95% CI Lower</th><th>95% CI Upper</th><th>Est.</th><th>95% CI Lower</th><th>95% CI Upper</th></tr></thead><tbody><tr><td>log(ϑ)</td><td>3.68</td><td>3.56</td><td>3.81</td><td>3.62</td><td>3.05</td><td>4.19</td><td>3.68</td><td>3.56</td><td>3.81</td></tr><tr><td>log(f̃<sub>M<sub>2</sub></sub>)</td><td>1.36</td><td>1.32</td><td>1.40</td><td>1.35</td><td>1.01</td><td>1.69</td><td>1.36</td><td>1.32</td><td>1.40</td></tr><tr><td>log(f̃<sub>M<sub>3</sub></sub>)</td><td>2.51</td><td>2.48</td><td>2.54</td><td>2.62</td><td>2.31</td><td>2.93</td><td>2.51</td><td>2.48</td><td>2.54</td></tr><tr><td>log(δ<sub>S</sub>)</td><td>-4.72</td><td>-4.74</td><td>-4.70</td><td>-4.72</td><td>-4.92</td><td>-4.51</td><td>-4.72</td><td>-4.74</td><td>-4.70</td></tr><tr><td>log(δ<sub>Ab</sub>)</td><td>-1.89</td><td>-2.09</td><td>-1.69</td><td>-1.91</td><td>-2.61</td><td>-1.21</td><td>-1.89</td><td>-2.09</td><td>-1.69</td></tr><tr><td>ω<sub>ϑ</sub></td><td>0.47</td><td>0.32</td><td>0.61</td><td>0.46</td><td>0.34</td><td>0.58</td><td>0.47</td><td>0.32</td><td>0.61</td></tr><tr><td>ω<sub>f̃<sub>M<sub>2</sub></sub></sub></td><td>0.07</td><td>0.05</td><td>0.08</td><td>0.17</td><td>0.04</td><td>0.31</td><td>0.07</td><td>0.05</td><td>0.08</td></tr><tr><td>σ<sub>ε</sub></td><td>0.22</td><td>0.21</td><td>0.23</td><td>0.21</td><td>0.19</td><td>0.23</td><td>0.22</td><td>0.21</td><td>0.23</td></tr></tbody></table>