# A statistical approach to latent dynamic modeling with differential
equations

Maren Hackenberg
Institute of Medical Biometry and Statistics,
Faculty of Medicine and Medical Center, University of Freiburg
and
Astrid Pechmann
Department of Neuropediatrics and Muscle Disorders,
Faculty of Medicine and Medical Center, University of Freiburg
and
Clemens Kreutz
Institute of Medical Biometry and Statistics,
Faculty of Medicine and Medical Center, University of Freiburg;
Freiburg Center for Data Analysis and Modeling, University of Freiburg;
Centre for Integrative Biological Signaling Studies (CIBSS), University of
Freiburg
and
Janbernd Kirschner
Department of Neuropediatrics and Muscle Disorders,
Faculty of Medicine and Medical Center, University of Freiburg
and
Harald Binder
Institute of Medical Biometry and Statistics,
Faculty of Medicine and Medical Center, University of Freiburg;
Freiburg Center for Data Analysis and Modeling, University of Freiburg;
Centre for Integrative Biological Signaling Studies (CIBSS), University of
Freiburg
This work was supported by the DFG (German Research Foundation) – project ID
322977937/GRK 2344 (MH) and project ID 499552394/CRC 1597 (HB, CK and MH). AP
was supported by the Berta-Ottenstein clinician scientist program of the
University of Freiburg. Biogen and Novartis provide financial support for the
SMArtCARE registry.

###### Abstract

Ordinary differential equations (ODEs) can provide mechanistic models of
temporally local changes of processes, where parameters are typically informed
by external knowledge. While ODEs are popular in systems modeling, they are
less established for statistical modeling of longitudinal cohort data, e.g.,
in a clinical setting. Yet, modeling of local changes could also be attractive
for assessing the trajectory of an individual in a cohort in the immediate
future given its current status, where ODE parameters could be informed by
further characteristics of the individual. However, several hurdles so far
limit such use of ODEs, as compared to regression-based global function
fitting approaches. First, the potentially higher level of noise in cohort
data might be detrimental to ODEs, as the shape of a function obtained from
solving an ODE heavily depends on the — then rather noisy — initial value.
Second, larger numbers of variables multiply such problems and might be
difficult to handle for ODEs. To address this, we propose to use each
observation in the course of time as the initial value to obtain multiple
local ODE solutions and build a combined estimator of the underlying dynamics.
Neural networks are used for obtaining a low-dimensional latent space for
dynamic modeling from a potentially large number of variables, and for
obtaining patient-specific ODE parameters from baseline variables.
Simultaneous identification of dynamic models and of a latent space is enabled
by recently developed differentiable programming techniques. We illustrate the
proposed approach in an application with spinal muscular atrophy (SMA)
patients and a corresponding simulation study. In particular, modeling of
local changes in health status at any point in time is contrasted to the
interpretation of functions obtained from global fitting via regression. This
more generally highlights how different application settings might demand
different modeling strategies.

Keywords: latent representations, deep learning, differentiable programming,
longitudinal data

##  1 Introduction

Different modeling communities have developed distinct approaches to
describing temporal processes that underlie longitudinal data, such as the
progression of a disease over time. In statistics, analysis of longitudinal
cohort data is typically based on regression techniques (Fitzmaurice et al.,
2008; Bandyopadhyay et al., 2011), which provide a global model that is a good
fit on average over the entire observed time course, e.g., corresponding to
the average course of a typical individual from a given cohort. Yet, sometimes
the future trajectory of a specific individual, given the current status,
might be of interest. For example, a clinical practitioner might want to
predict how a patient’s health status will develop until the next follow-up
visit given the current status. This corresponds to focusing on relative
changes and a local perspective on the dynamics, as in ordinary differential
equations (ODEs). The latter are commonly used, e.g., in systems biology
(Bruggeman & Westerhoff, 2007; Banga, 2008; DiStefano, 2015), but less
established in the statistics community. ODEs describe a small set of
quantities by carefully specified relations, typically based on domain
knowledge.

In some settings, information about the dynamics are known beforehand. For
example, in a systems biology context, kinetic parameters of biochemical
reactions might be known from thermodynamics and initial conditions depend on
the experimental context. In a clinical cohort setting, a typical disease
progression might be determined by patients’ characteristics at baseline,
e.g., the age at diagnosis. In such settings, the observation at the first
time point could be used as an initial value, which will then strongly
influence the shape of the solution. This might be problematic particularly
for noisy longitudinal cohort data (Xue et al., 2010; Verbeke et al., 2014).
In addition, specifying the relations of a larger number of variables with
ODEs is challenging, as sufficient domain knowledge may not be available and
the resulting complex systems may be difficult to solve numerically.

To still enable modeling with ODEs for describing individual local changes in
such a setting, yet reduce the emphasis on the initial value, we propose a
statistically inspired approach based on ODEs. In addition, we propose to
integrate neural network techniques for dimension reduction of a larger number
of variables, and for inferring individual ODE parameters from the
characteristics of an individual.

Neural networks allow for building more flexible approaches, e.g., for
modeling with unobserved quantities in a latent representation. Such
combinations of data-driven modeling based on neural networks with knowledge-
driven modeling based on ODEs are explored, e.g., in the framework of
universal differential equations (Rackauckas et al., 2020) based on the idea
of neural ODEs (Chen et al., 2018; Rubanova et al., 2019), and have also been
investigated in biomedical applications, e.g., for intervention modeling (Gwak
et al., 2020), survival analysis (Moon et al., 2022), or modeling of
psychological resilience (Köber et al., 2022a, b). Such a combination of
modeling components is facilitated by differentiable programming, a paradigm
that allows for joint optimization of different model components via automatic
differentiation (Baydin et al., 2017). While originally proposed to bridge
scientific computing and machine learning (Innes et al., 2019), we have argued
that such integrative approaches can also be applied to statistical modeling
(Hackenberg et al., 2022a).

In settings with more classical statistical models, larger numbers of
variables have been reduced to a single summary score (Bandyopadhyay et al.,
2011; Fitzmaurice et al., 2008) and subsequently analyzed with univariate
methods. Yet, this means that dimension reduction and longitudinal modeling
are performed separately. In some approaches, the covariance structure between
outcomes has been included explicitly or implicitly into longitudinal
regression (Fan et al., 2007; Bandyopadhyay et al., 2011; Fitzmaurice et al.,
2008). For example, functional principal component analysis on the covariance
surface can be used to identify a smooth joint trajectory and individual
deviations (Nolan et al., 2023), or mixed-effects models can be adapted, e.g.,
for explicitly modeling intra-individual variability (Palma et al., 2023), or
for capturing individual trajectories via spatiotemporal transformations
(Schiratti et al., 2017). In contrast, auto-regressive approaches dynamically
adapt to incoming data by regressing the observed variable on its past (Van
Montfort et al., 2018). For modeling in a latent space, long-short term memory
(LSTM) models infer a latent state from observed data using neural networks,
which is then iteratively updated (Hochreiter & Schmidhuber, 1997; Graves &
Graves, 2012). Yet, these models operate in discrete time and typically rely
on a fixed time grid of observations. While there are some examples of
continuous-time modeling in statistical applications (Ramsay & Hooker, 2017;
de Haan-Rietdijk et al., 2017; Mews et al., 2022), often the data is too
sparse and noisy to allow for fitting ODEs (Van Montfort et al., 2018),
resulting in limited use in the statistics community so far.

Conversely, statistical function fitting techniques have been proposed to
estimate ODE parameters from an observed time series in systems modeling
(Liang & Wu, 2008; Gábor & Banga, 2015; Ramsay & Hooker, 2017), e.g., using
maximum likelihood estimation and uncertainty quantification via the profile
likelihood (Raue et al., 2009; Kreutz et al., 2012). Parameter estimation is
then based on longitudinal information, i.e., corresponds to the function
fitting approach of regression techniques.Initial values can also be
considered as additional parameters to account for observation noise and can
then be estimated jointly with the parameters specifying the dynamics. Yet,
such models are typically used in settings with a smaller number of variables
and not designed for modeling in a latent space. In addition, it could be
attractive to obtain the ODE parameters not based on function fitting, but
based on external information from the individuals, to model each time series
individually. This might be particularly important when modeling data with
heterogeneous individual dynamics.

We therefore propose an approach where ODE parameters are not obtained via
function fitting, but from baseline characteristics of individuals via a
neural network. For modeling with a larger number of variables, we fit the
dynamic model in a low-dimensional latent space learned by a neural network,
specifically a variational autoencoder (VAE) (Kingma & Welling, 2014), to
reflect the assumption of a lower-dimensional underlying dynamic process
driving the observed quantities. To reduce the dependence on the initial value
and increase robustness to noise, we propose to solve multiple ODEs, using
each observed value as the initial value and averaging the solutions using
time-dependent inverse-variance weights.

Before explaining the proposed approach in detail, we illustrate conceptual
differences between global function fitting approaches and ODE-based local
approaches, motivated by an application example from the SMArtCARE rare
disease registry (Pechmann et al., 2019), which highlights how different
research questions demand different modeling perspectives, thus more generally
motivating our proposed approach. After providing an in-depth description of
the approach, we use the SMArtCARE data and a corresponding simulation design
to contrast the resulting individual-specific local models of changes in
dynamics at any point in time with a global function fitting approach.
Finally, we discuss limitations and the more general potential of fusing
dynamic modeling approaches across communities.

##  2 An illustration of local and global perspectives in a disease registry
application

We illustrate the individual-level temporally local perspective on dynamic
modeling as compared to temporally global function fitting approaches in the
context of the SMArtCARE registry, a prospective multicenter cohort study
where longitudinal data on SMA patients’ disease developments is collected
during routine visits (Pechmann et al., 2019). As the registry has been set up
relatively recently, there are few time points per patient with irregular
timing and frequency of follow-up visits. In addition to an extensive baseline
characterization, motoric ability is assessed with physiotherapeutic tests at
follow-up visits.

SMA is caused by a homozygous deletion in the survival motor neuron (SMN) 1
gene on chromosome 5, which is crucial for normal motor neuron function
(Pechmann et al., 2019). As a result of SMN deficiency, muscles do not receive
signals from motor neurons, leading to atrophy, i.e., muscle degeneration.
Treatment is designed to increase production of the SMN protein to maintain
functional motor neurons. These underlying disease processes cannot be
observed directly, but are implicitly reflected in the motor function
assessments, thus motivating modeling of latent dynamics that drive the
observed measurements. We seek to identify these latent disease dynamics from
the observed motor function using a VAE-based dimension reduction and model
them conditional on patients’ characteristics such as age, treatment and the
extent of the SMN deficiency, which can be informative of the disease
dynamics.

For more closely analyzing the corresponding modeling challenge, we consider
the observed time series of a hypothetical individual patient. In this
setting, fitting a regression model corresponds to estimating a global
function that is the best fit through the observed points on average, e.g., by
minimizing squared distances (Figure 1a). In contrast, an ODE model describes
the rate of change relative to the current observed value as an explicit
mechanistic function (denoted as $g$ in Figure 1b), i.e., locally in time. For
predictions in subsequent time intervals, the ODE is solved using the last
observed value as initial condition (Figures 1b), while the average fit of all
previous observations on an absolute level is used in regression-based
function fitting (Figure 1a, function fitting approach).

Figure 1: Different perspectives on modeling dynamics: Regression-based
function fitting approaches provide global fits of observed values on an
absolute level (blue circles) (a), where the average fit (red solid line) is
extrapolated (top red dotted line) to make predictions for new time points. In
contrast, ODE-based approaches which use observations as initial values
provide local fits by specifying the rate of change relative to the current
status, which is used as the initial value for starting the ODE (b). Applying
a similar perspective to a regression approach corresponds to shifting the
regression fit to the current value and using that as initial value for
starting the extrapolation (Panel (a), bottom red dotted line).

In this simplified setting, the regression approach can also be turned into an
initial value approach by shifting the globally fitted function to the last
observed value, which can then be used as starting point for extrapolation
(Figure 1a, initial value approach). For example, in Figure 1, the slope
equation of the regression is an approximation of the derivative of $f$
modeled by the differential equation for $\Delta\to 0$, i.e., in the limit
both models are equivalent. Yet, regression coefficients do not translate to
parameters of a local rate of change in general, but represent a conceptually
distinct global function fitting approach.

Global function fitting considers the observed values as measurements with
noise, which is to be averaged out, within or often also across individuals,
before extrapolating. For example, measuring the distance a SMA patient can
walk within a given time frame will likely be subject to variability due to
the patients’ daily form and motivation. Using a global average for prediction
in a regression-based approach then improves robustness to such random
variations. In contrast, the last observation can be considered as most
accurate approximation of the true value of interest in ODE modeling. For
example, when a child acquired a new motoric skill, extrapolating the future
development based on having that skill might yield more meaningful predictions
than averaging over all past time points when the child has not yet had that
skill. In addition, the clinical practitioner might be more interested in
anticipated relative changes, which are more directly reflected in ODE models,
instead of the absolute levels of observations. A typical function fitting
approach cannot account for external changes that are not described by the
model, whereas these can be taken into account by restarting the ODE based on
the last observed value when a new observation is made. However, this comes at
the price of a higher vulnerability to noise, illustrated by the jumps in
Figure 1b, as opposed to a smooth global regression function.

Another difference becomes apparent when assuming that model parameters of
ODEs can be obtained from some other source instead of from direct function
fitting. As function fitting for regression models is performed by minimizing
the average distances of the absolute levels of observations over multiple
time points and maybe even across multiple individuals, model fits will become
more robust for longer observed time series, and less reliable when only a few
data points are available. For example, this is problematic in the early phase
of a clinical cohort study, such as the SMArtCARE registry, when only a small
number follow-up times is available. In contrast, when assuming that ODE
parameters can be obtained from external information such as baseline
characteristics, the ODE is completely specified when just one initial
observation is given as starting value, and can be used to predict subsequent
developments. For example, patients’ baseline characteristics could
potentially be informative about subsequent individual disease progression in
our SMA application, and specifically these connections might be of interest
to biomedical researchers. While such an approach of determining model
parameters based on baseline characteristics could, in principle, be developed
for global regression models, this would imply that also the absolute level of
the values could be predicted based on baseline characteristics, whereas in
ODEs, they only need to inform a model for individual changes.

Spelling out such conceptual differences can help to choose between the
different approaches, depending on the application scenario and question to be
answered. For example, the aim in our application is to anticipate an
individual child’s future motoric development based on present motoric skills,
which favors the proposed ODE approach.

##  3 Methods

###  3.1 Statistical modeling with ODEs and many initial values

Our proposed approach for estimating dynamics is based on solving a linear ODE
system multiple times, using each value observed in the course of time as the
initial condition. The resulting individual solutions are combined into an
inverse-variance weighted average using time-dependent weights, to obtain an
unbiased minimum variance estimator of the true underlying dynamics. We assume
a low-dimensional time series of observations ${\bf z}_{t_{0}},\dots,{\bf
z}_{t_{K}}$ with ${\bf z}_{t}\in\mathbb{R}^{d}$ for all $t=t_{0},\dots,t_{K}$,
which can be thought of as noisy measurements of a true, deterministic
underlying process
${\boldsymbol{\mu}}:\mathbb{R}_{+}\to\mathbb{R}^{d},t\mapsto{\boldsymbol{\mu}}(t)$,
i.e., ${\bf z}_{t}={\boldsymbol{\mu}}(t)+{\boldsymbol{\epsilon}}(t)$ with
${\boldsymbol{\epsilon}}(t)\sim\mathcal{N}(0,{\boldsymbol{\sigma}}^{2}(t))$.
For simplicity, we assume that measurement errors are independent of time,
i.e.,
${\boldsymbol{\epsilon}}(t)\equiv{\boldsymbol{\epsilon}}\sim\mathcal{N}(0,{\boldsymbol{\sigma}}^{2})$
for all $t\in\mathbb{R}_{+}$ for some fixed
${\boldsymbol{\sigma}}^{2}\in\mathbb{R}^{d}$, and thus ${\bf
z}_{t}\sim\mathcal{N}({\boldsymbol{\mu}}(t),{\boldsymbol{\sigma}}^{2})$ for
all $t=t_{0},\dots,t_{K}$. We assume that the underlying dynamics can be
described by a linear ODE system with constant coefficients, and consider for
each $k=0,\dots,K$ the initial value problem

$\displaystyle\begin{split}\frac{d}{dt}{\boldsymbol{\mu}}(t)&=A\cdot{\boldsymbol{\mu}}(t)+{\bf
c};\\\ {\boldsymbol{\mu}}(t_{k})&={\bf z}_{k},\end{split}$ (1)

where the $k$-th observation is used as the initial condition. For $A=0$, the
solution to (1) is given by ${\boldsymbol{\mu}}(t)={\bf c}\cdot t+{\bf
z}_{k}$. For $A\neq 0$, the solution can be computed analytically (see
Supplementary Section 6.2 for a detailed derivation):

$\displaystyle{\boldsymbol{\mu}}(t)=\exp(A(t-t_{k}))\cdot(A^{-1}{\bf c}+{\bf
z}_{k})-A^{-1}{\bf c}.$ (2)

Note that for non-constant coefficients, an analytical solution can be
obtained analogously but requires numerical evaluation of integral terms that
often do not have closed-form solutions. Slightly overloading notation, we
define the solutions of Equation (1) for $A\neq 0$ as ${\bf
z}_{k}(t,{\boldsymbol{\eta}}):=(A^{-1}{\bf c}+{\bf
z}_{k})\exp(A(t-t_{k}))-A^{-1}{\bf c}$ for each $k=0,\dots,K$, where
${\boldsymbol{\eta}}=\\{A,{\bf c}\\}$ and the initial value is given by ${\bf
z}_{k}$. The solutions are then combined to a time-dependent inverse-variance
weighted average as an estimator for true underlying dynamics:

$\displaystyle
g(t,{\boldsymbol{\eta}}):=\frac{\sum_{k=0}^{K}\mathrm{Var}\left[{\bf
z}_{k}(t,{\boldsymbol{\eta}})\right]^{-1}{\bf
z}_{k}(t,{\boldsymbol{\eta}})}{\sum_{k=0}^{K}\mathrm{Var}\left[{\bf
z}_{k}(t,{\boldsymbol{\eta}})\right]^{-1}}.$ (3)

As $\mathbb{E}\left[{\bf
z}_{k}(t,{\boldsymbol{\eta}})\right]={\boldsymbol{\mu}}(t_{k})$ for all
$k=0,\dots,K$, the estimator $g(t,{\boldsymbol{\eta}})$ is unbiased. Further,
it has minimum variance among all weighted average estimators (Shahar, 2017).
For calculating the weights, we have

$\displaystyle\begin{split}\mathrm{Var}\left[{\bf
z}_{k}(t,{\boldsymbol{\eta}})\right]&=\mathrm{Var}\left[\exp(A(t-t_{k}))\cdot(A^{-1}{\bf
c}+{\bf z}_{k})-A^{-1}{\bf c}\right]\\\
&=\mathrm{Var}\left[\exp(A(t-t_{k})\cdot{\bf z}_{k}\right]\\\
&=\exp(2A(t-t_{k}))\cdot\mathrm{Var}\left[{\bf z}_{k}\right]\\\
&=\exp(2A(t-t_{k}))\cdot{\boldsymbol{\sigma}}^{2},\end{split}$ (4)

which we can insert in (3) to obtain an explicit closed-form expression.

###  3.2 Modeling dynamics in a latent representation

When there is a large number of variables that are observed over the course of
time, often a model for an underlying lower-dimensional process is sought that
drives the observed measurements, instead of separately modeling the dynamics
of each observed variable. For example, in our SMA application, the underlying
disease dynamics are driven by degradation of motor neurons causing muscle
degeneration, which cannot be observed directly but is implicitly reflected in
various items of motor function tests. To model such latent processes, we
suggest to use neural networks for learning a low-dimensional representation
of observed measurements, and integrating our statistical approach for dynamic
modeling with ODEs into such a latent representation via simultaneous fitting
of the dynamic model and the neural networks.

####  3.2.1 Using VAEs for dimension reduction

For dimension reduction, we use a variational autoencoder (VAE), a generative
deep learning model that infers a compressed, low-dimensional representation
of the data using neural networks (Kingma & Welling, 2014, 2019).
Specifically, the latent space is defined as a low-dimensional random variable
${\bf Z}$ with prior distribution $P^{{\bf Z}}$. Two distinctly parameterized
neural networks, called the encoder and the decoder, map an observation ${\bf
x}\in\mathbb{R}^{p}$ to the latent space and back to data space. The
parameters of the encoder and decoder are jointly optimized to infer a low-
dimensional representation, based on which the original data can be well
reconstructed. Formally, the encoder and decoder parameterize the conditional
distributions $q_{{\bf Z}\mid{\bf x}}(\cdot,{\boldsymbol{\varphi}})$ and
$p_{{\bf X}\mid{\bf z}}(\cdot,{\boldsymbol{\theta}})$, where we abbreviate
conditional distributions and densities by writing, e.g., $p_{{\bf X}\mid{\bf
z}}({\bf x},{\boldsymbol{\theta}})$ for $p_{{\bf X}\mid{\bf Z}={\bf z}}({\bf
x},{\boldsymbol{\theta}})$. The model is trained, i.e., the parameters
${\boldsymbol{\theta}}$ and ${\boldsymbol{\varphi}}$ of the encoder and
decoder networks are optimized, by maximizing the evidence lower bound (ELBO),
a lower bound on the data likelihood $p_{{\bf X}}$, derived based on
variational inference (Blei et al., 2017):

$\mathrm{ELBO}({\bf
x},{\boldsymbol{\varphi}},{\boldsymbol{\theta}})=E_{q_{{\bf Z}\mid{\bf
x}}(\cdot,{\boldsymbol{\varphi}})}[\log(p_{{\bf X}\mid{\bf z}}({\bf
x},{\boldsymbol{\theta}}))]-D_{\mathrm{KL}}(q_{{\bf Z}\mid{\bf
x}}(\cdot,{\boldsymbol{\varphi}})\|p_{{\bf Z}}).$

It can be shown that maximizing the ELBO is equivalent to minimizing the
Kullback-Leibler divergence between the true but intractable posterior
$p_{{\bf Z}\mid{\bf x}}(\cdot,{\boldsymbol{\theta}})$ and its approximation by
a member of a parametric variational family $\\{q_{{\bf Z}\mid{\bf
x}}(\cdot,{\boldsymbol{\varphi}})\\}$, typically assumed as Gaussian with
diagonal covariance matrix (Blei et al., 2017; Kingma & Welling, 2014, 2019).
The first term of the ELBO can be interpreted as a reconstruction error, while
the second term acts as a regularizer that encourages densities close to the
standard normal prior.

Training the VAE then corresponds to maximizing the ELBO as a function of the
parameters ${\boldsymbol{\theta}}$ and ${\boldsymbol{\varphi}}$ of the encoder
and decoder, specifying the conditional distributions, which yields both an
approximate maximum likelihood estimate for ${\boldsymbol{\theta}}$ and an
optimal variational density $q_{{\bf Z}\mid{\bf
x}}(\cdot,{\boldsymbol{\varphi}})$. Such optimization is typically performed
by stochastic gradient descent (Kingma & Ba, 2015), where the so-called
reparameterization trick is used to obtain gradients of the ELBO w.r.t. the
variational parameters ${\boldsymbol{\varphi}}$ (Kingma & Welling, 2014,
2019).

####  3.2.2 Incorporating dynamics into the loss function

Building on an approach that we have proposed previously (Hackenberg et al.,
2022b), we integrate the dynamic model described in Section 3.1 in the VAE
latent space. We use observed time series of patients’ measurements as input
to the model, represented as a matrix ${\bf x}_{i}$ of $T_{i}+1$ measurements
of $p$ variables for individual $i$, where $t_{0}$ is the common baseline time
point, and $t_{1}^{i},\dots,t_{T_{i}}^{i}$ are individual-specific subsequent
measurement time points. In our application setting from clinical cohort
studies, typically a more extensive patient characterization at the baseline
time point is available, such as age at symptom onset or treatment start,
which we assume to be informative of intra-individual differences in
underlying disease dynamics.

We specify the variational posterior as multivariate Gaussian with diagonal
covariance matrix, parameterized by the VAE encoder. Specifically, the encoder
maps an observed time series column-wise to the posterior mean
${\boldsymbol{\mu}}_{i}=({\boldsymbol{\mu}}_{i}^{t_{0}},{\boldsymbol{\mu}}_{i}^{t_{1}^{i}},\dots,{\boldsymbol{\mu}}_{i}^{t_{T_{i}}^{i}})\in\mathbb{R}^{m\times
T_{i}+1}$ and standard deviation
${\boldsymbol{\sigma}}_{i}=({\boldsymbol{\sigma}}_{i}^{t_{0}},{\boldsymbol{\sigma}}_{i}^{t_{1}^{i}},\dots,{\boldsymbol{\sigma}}_{i}^{t_{T_{i}}^{i}})\in\mathbb{R}^{m\times
T_{i}+1}$. We assume that the dynamics of the latent posterior mean are
governed by a linear ODE system and use the approach described in Section 3.1
to calculate the estimator from Equation (3) and define

$\widetilde{{\boldsymbol{\mu}}_{i}}(t,{\boldsymbol{\eta}}_{i}):=\frac{\sum_{k=0}^{T_{i}}\mathrm{Var}[\widetilde{{\boldsymbol{\mu}}}_{i,k}(t,{\boldsymbol{\eta}}_{i})]^{-1}\widetilde{{\boldsymbol{\mu}}}_{i,k}(t,{\boldsymbol{\eta}}_{i})}{\sum_{k=0}^{T_{i}}\mathrm{Var}[\widetilde{{\boldsymbol{\mu}}}_{i,k}(t,{\boldsymbol{\eta}}_{i})]^{-1}}.$
(5)

Here $\widetilde{{\boldsymbol{\mu}}}_{i,k}(t,{\boldsymbol{\eta}}_{i})$ is the
analytical ODE solution obtained according to Equation (2), using the encoded
value ${\boldsymbol{\mu}}_{i}^{t_{k}^{i}}$ from the $k$-th time point of the
$i$-th individual as initial value.

To estimate the unknown variance
$\mathrm{Var}[\widetilde{{\boldsymbol{\mu}}}_{i,k}(t,{\boldsymbol{\eta}}_{i})]$
at a time point $t$ of the ODE solution with initial value
${\boldsymbol{\mu}}_{i}^{t_{k}^{i}}$, we use all encoded values
${\boldsymbol{\mu}}_{i}^{t_{j}^{i}}$ for $t_{k}^{i}<t_{j}^{i}<t$, i.e., from
time points between the current initial time point ${t_{k}^{i}}$ and the time
point $t$ of interest, and calculate the sample variance of the corresponding
ODE solutions
$\widetilde{{\boldsymbol{\mu}}}_{i,j}(t,{\boldsymbol{\eta}}_{i})$ for all $j$
with $t_{k}^{i}<t_{j}^{i}<t$.

To obtain personalized dynamics conditional on baseline information, we use an
additional neural network to map the baseline variables to individual ODE
parameters ${\boldsymbol{\eta}}_{i}$. We subsequently use the estimator
$\widetilde{{\boldsymbol{\mu}}_{i}}(t,{\boldsymbol{\eta}}_{i})$ of the
underlying dynamics for sampling ${\bf
z}_{i}^{t_{k}^{i}}\sim\mathcal{N}(\widetilde{{\boldsymbol{\mu}}_{i}}(t_{k}^{i},{\boldsymbol{\eta}}_{i}),{\boldsymbol{\sigma}}_{i}^{t_{k}^{i}})$
for $k=0,\dots,T_{i}$, and passing it to the VAE decoder to obtain a
reconstructed time series $\widehat{{\bf x}}_{i}$.

For jointly optimizing the VAE and the neural network for the ODE parameters,
we maximize the ELBO, using
$\widetilde{{\boldsymbol{\mu}}_{i}}(t,{\boldsymbol{\eta}}_{i})$ as the latent
posterior mean. We add the squared Euclidean distance between the posterior
mean ${\boldsymbol{\mu}}_{i}$, obtained directly from the encoded time series,
and the posterior mean constrained to smooth dynamics
$\widetilde{{\boldsymbol{\mu}}_{i}}(t,{\boldsymbol{\eta}}_{i})$ to encourage
consistency of the latent representation before and after solving the ODEs. To
regularize the inverse-variance weights of the ODE solution estimator, we
further add the log-differences between the sample variances $s^{2}$ of all
solutions
$\widetilde{{\boldsymbol{\mu}}}_{i,k}(t,{\boldsymbol{\eta}}_{i}),k=0,\dots,T_{i}$,
and of all encoded values ${\boldsymbol{\mu}}_{i}$, summed across all observed
time points, where an offset of $0.1$ is added before taking the logarithm to
avoid values near zero. The final loss function is given by

$\begin{split}\mathcal{L}({\bf
x}_{i},{\boldsymbol{\eta}}_{i},{\boldsymbol{\theta}},{\boldsymbol{\varphi}})&=D_{\mathrm{KL}}(\widetilde{q}_{{\bf
Z}\mid{\bf
x}_{i}}(\cdot,{\boldsymbol{\eta}}_{i},{\boldsymbol{\varphi}})\|p_{{\bf
Z}}(\cdot,{\boldsymbol{\theta}}))-\mathrm{E}_{\widetilde{q}}[\log(p_{{\bf
X}\mid{\bf z}_{i}}({\bf x}_{i},{\boldsymbol{\theta}}))]\\\
&+\alpha\|{\boldsymbol{\mu}}_{i}-\widetilde{{\boldsymbol{\mu}}_{i}}\|_{2}^{2}\\\
&+\beta\sum_{k=0}^{T_{i}}\log(s^{2}((\widetilde{{\boldsymbol{\mu}}}_{i,j}(t_{k},{\boldsymbol{\eta}}_{i}))_{j=0,\dots,T_{i}})+0.1)-\log(s^{2}(({\boldsymbol{\mu}}_{i}^{t_{j}})_{{j=0,\dots,T_{i}}})+0.1)\end{split}$
(6)

where $\alpha,\beta\in[0,1]$ are hyperparameters balancing the loss
components.

We optimize the joint loss function from Equation (6), which simultaneously
incorporates all model components, i.e., the dynamic model, the VAE for
dimension reduction, and the neural network for obtaining the ODE parameters
from baseline characteristics, by stochastic gradient descent. We use
automatic differentiation (Baydin et al., 2017) to simultaneously obtain
gradients with respect to all parameters, including the VAE encoder and
decoder parameters ${\boldsymbol{\theta}},{\boldsymbol{\varphi}}$ and the
individual-specific dynamic model parameters ${\boldsymbol{\eta}}_{i}$. This
can be realized efficiently using the Zygote.jl package (Innes et al., 2019)
in the Julia programming language (Bezanson et al., 2017), a flexible
framework for automatic differentiation that allows to implement the model and
loss function with only minimal code adaptation for automatically obtaining
gradients. The model is implemented as a publicly available Julia package
(https://github.com/maren-ha/LatentDynamics.jl), including Jupyter notebooks
to illustrate the approach. Further implementation details can be found in the
supplementary material, Section 6.1.

##  4 Evaluation

###  4.1 Simulation

We empirically evaluate our approach in a simulation design and in an
application with the SMArtCARE rare disease registry on SMA patients’
development of motoric ability. A general challenge for evaluation is that the
latent representation is invariant to affine linear transformations such, as
shifting, scaling and rotation, as these can be reversed in the decoder, such
that scaling in latent space is arbitrary. In this setting, our approach
should, e.g., distinguish patterns with different monotonicity behavior as a
minimum requirement. For example, in SMA there might be two underlying
processes corresponding to, e.g., motor neuron degradation and muscle
function, and a model should be able to distinguish between groups of patients
with different development trends.

To investigate this property in a simulation study, we adapt our previous
design from Hackenberg et al. (2022b). There, we defined two groups of
individuals with distinct underlying development patterns, corresponding to
two distinct sets of ODE parameters, based on a homogeneous two-dimensional
linear ODE system, i.e., with four unknown parameters. The two ODE systems are
given by

$\begin{split}\frac{d}{dt}\begin{pmatrix}u_{1}\\\
u_{2}\end{pmatrix}(t)&=\begin{pmatrix}-0.2&0.1\\\
-0.1&0.1\end{pmatrix}\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}(t);\\\
\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}(0)&=\begin{pmatrix}3\\\
1\end{pmatrix}\end{split}$

$\begin{split}\frac{d}{dt}\begin{pmatrix}u_{1}\\\
u_{2}\end{pmatrix}(t)&=\begin{pmatrix}-0.2&-0.1\\\
0.1&-0.2\end{pmatrix}\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}(t);\\\
\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}(0)&=\begin{pmatrix}3\\\
1\end{pmatrix}.\end{split}$

We simulated $n=100$ individuals split equally into the two groups. For each
individual $i=1,\dots,n$, we sampled a random number $T_{i}$ of between $1$
and $8$ follow-up observations after the common baseline time point $t_{0}$ at
random time points $t_{1}^{i},\dots,t_{T_{i}}^{i}$, sampled uniformly from
$[1.5,10]$. At each time point $t_{k}^{i}$, we simulated measurements of
$p=10$ variables by adding a variable-specific measurement error
$\delta_{j}^{t_{k}^{i}}\sim\mathcal{N}(0,\sigma_{\mathrm{var}}^{2})$ for
$j=1,\dots,p$ and an individual-specific measurement error
$\varepsilon_{i,j}^{t_{k}^{i}}\sim\mathcal{N}(0,\sigma_{\mathrm{ind}}^{2})$
for $i=1,\dots,n$, $j=1,\dots,p$ to the true value of the ODE solution. For an
individual $i$, we then obtained simulated observations ${\bf
x}_{i}\in\mathbb{R}^{p\times T_{i}+1}$ by defining

$x_{i,j}^{t_{k}}:=u_{1}(t_{k})+\delta_{j}^{t_{k}}+\varepsilon_{i,j}^{t_{k}}\quad\text{for
}k=0,\dots,T_{i}+1\text{ and }j=1,\dots,5,$

and

$x_{i,j}^{t_{k}}:=u_{2}(t_{k})+\delta_{j}^{t_{k}}+\varepsilon_{i,j}^{t_{k}}\quad\text{for
}k=0,\dots,T_{i}+1\text{ and }j=6,\dots,10.$

We additionally simulated baseline variables by sampling with variance
$\sigma_{\mathrm{info}}^{2}$ from the simulated individual’s true ODE
parameters to obtain $10$ informative baseline variables and add $40$ noise
variables sampled from a centered Gaussian distribution with variance
$\sigma_{\mathrm{noise}}^{2}$, such that we ended up with a total of $q=50$
baseline variables. We set $\sigma_{\mathrm{ind}}^{2}=0.5$ and
$\sigma_{\mathrm{var}}^{2}=\sigma_{\mathrm{info}}^{2}=\sigma_{\mathrm{noise}}^{2}=0.1$.

On this data, we trained the model described in the previous section and
extracted the learnt trajectories. The focus of our approach is to locally
predict changes in latent dynamics, given the last measurement.
Correspondingly, we visualize the learnt trajectories by starting at the
learnt latent representation of each observed time point consecutively and
solving an ODE with the learnt parameters until the next time point, thus
reflecting our aim of predicting changes in the immediate future based on the
current status until a new measurement becomes available. In Figure 2, we show
the true underlying ODE solutions of both groups in panel (a) and the fitted
latent trajectories (solid lines) based on the mean of the latent
representation for $12$ exemplary simulated individuals in panel (b). The
colored bands around the trajectory correspond to the learnt standard
deviation of the latent representation. The results show that the approach
allows for recovering group-specific underlying dynamic patterns (dashed
lines) and the local predictions at the subsequent time point mostly fall
within the range of one standard deviation (colored bands) of the next encoded
observation (filled circles). Note that the fitted trajectories sometimes
appear shifted or downscaled (in particular in the second dimension), due to
the model freedom in structuring its latent representation mentioned above and
the $\mathcal{N}(0,1)$ prior on the latent variable.

Figure 2: (a) Ground-truth trajectories for two distinct development patterns
based on which data was simulated. (b) Exemplary fits of simulated individuals
from both groups. One panel corresponds to one individual, the $x$ axis shows
the time and the $y$ axis the values of the two dimensions of the latent
representation. Circles depict the mean of the latent variables, as obtained
directly from the encoder, while solid lines indicate the ODE solutions, with
parameters obtained from the individual’s baseline variables, using the value
at the last observed time point as initial condition and solved until the
subsequent time point. Colored bands show the range of one standard deviation
of the latent variable around its mean. The two leftmost columns show $6$
individuals from the first group, corresponding to the left ground-truth
trajectory in (a), while the two rightmost columns show $6$ individuals from
the second group, corresponding to the pattern shown on the right in (a).

In addition, we wanted to empirically verify that our approach indeed
coincides with a linear regression fit when the underlying dynamics are
simple, i.e., it defaults to a least-squares fit (see Section 2). We thus
considered a simpler simulation scenario with a constant two-dimensional ODE-
system, i.e., resulting in linear fits, where estimating the rate parameter of
the ODE and the initial condition is equivalent to fitting a simple linear
regression with time as dependent variable separately for each latent
dimension. In this setting, local ODE fits from our model indeed closely match
a global linear regression fit of the encoded time series (see Supplementary
Figures 4 and 6). To fit this global regression model, the complete time
series has to be observed, while ODE solutions are obtained using baseline
information and the value at the last observed time point only. If only this
limited information is available and extrapolations have to be calculated
based on previously observed values only, the regression approach performs
worse, as reflected by a substantial drop in prediction performance compared
to the ODE approach, both in latent space and on the level of the
reconstructed items (see Supplementary Section 6.3).

###  4.2 Rare disease registry application

To illustrate our approach with real data, we use data from the SMArtCARE rare
disease registry. Specifically, we selected patients treated with Nusinersen
(Schorling et al., 2020) who completed the Revised Upper Limb Module (RULM)
test (Mazzone et al., 2017). This test comprises $20$ items with a maximum sum
score of $44$ to evaluate motor function of upper limbs and can be conducted
for all patients older than two years of age with the ability to sit in a
wheelchair (Pechmann et al., 2019). We used all test items as time-dependent
variables and all available baseline information, including, e.g., SMA
subtype, age at symptom onset and first treatment, and genetic test results.
For each patient, we removed outlier time points where the difference in RULM
sum score to the previous time point exceeded two times the interquartile
range of all sum score differences between adjacent time points. We then
filtered out $63$ patients with less than two observation time points, and
additionally filtered out $176$ patients with sum score variance smaller than
$1$, i.e., with nearly constant trajectories across all items. This resulted
in a dataset of $399$ patients with observations at between $2$ and $13$ time
points (median $7$ time points), corresponding to in total $2797$ observations
of $20$ time-dependent variables, in addition to $24$ baseline variables.
Integer-valued test items were rescaled and a logit transformation was applied
to account for the Gaussian generative distribution parameterized by the VAE
decoder.

Figure 3: Exemplary fits from SMArtCARE patients. Each panel shows the learnt
latent space trajectories of one patient, where the $x$ axis shows the time in
months and the $y$ axis the values of the two dimensions of the latent
representation. Circles depict the latent variable means obtained directly
from the VAE encoder, solid lines depict the fitted ODE solution using the
last observed value as initial condition, and solved until the subsequent
observation. Colored bands show a range of one standard deviation of the
latent variables around their mean. Least squares regression fits of all
previously observed latent values, shifted to the current observed time point
and extrapolated until the subsequent observation, are shown in pink.

In Figure 3, we show fitted latent trajectories from a two-dimensional ODE
system for the latent space for exemplary SMArtCARE patients. Analogous to
Figure 2, we display the ODE solutions for the posterior mean of the latent
representation with the learnt parameters, using the first encoded value as
the starting point and solving until the subsequent observation time point,
where the ODE is restarted with the new value as initial condition.

For comparison, we computed a least squares regression fit of all previously
observed encoded values at each observed time point , shifted the trajectory
to start at the current observation and extrapolated until the next
observation (red lines), as an adaptation of a regression-based approach to
local predictions analogous to Figure 1. Especially at the beginning of the
time series, the extrapolated values are often far off from the observed
value. For example, for the first observed time point, the current value can
only be carried forwards in time, i.e., predicting a constant trajectory. At
later time points, changes in monotony behavior lead to significant prediction
errors (e.g., in the third panel from left in the bottom row), while in
general predictions tend to get more accurate at later time points, when more
information has already been accumulated (e.g., as in the leftmost panel in
the second row).

This illustrates the conceptual difference in perspective to the ODE-based
approach, where the accuracy of the prediction does not depend on the amount
of previously observed values, as it relies on information from baseline
characteristics and acts locally as opposed to global models. Here, predicted
values mostly fall within the range of one standard deviation of the mean
trajectory (colored bands). Our proposed approach is able to capture
individual-specific patterns as reflected by individuals with various
different trends and both very dynamic and more constant patterns. Joint
optimization helps to regularize the latent representation, as the ODE model
imposes a smoothness constraint, thus encouraging a representation where
encoded values fall on a smooth trajectory. While this is effective in most
cases, there are still some outlier time points (e.g., in the second panel in
the top row) or large jumps that cannot be captured by a smooth model.

Quantitatively, our proposed approach yields lower prediction errors at
subsequent time points both in latent space and on item-level decoder
reconstructions of latent values and model predictions (mean squared error
averaged across all time points and patients of 1.386 (ODE) vs. 3.523 (least-
squares) in latent space, and 5.181 (ODE) vs. 20.89 (least-squares) on
reconstructed items).

##  5 Discussion

While regression techniques are commonly used for modeling dynamic processes
in statistics, ODE-based approaches dominate in the systems modeling
community. We have shown that different approaches correspond to conceptually
different perspectives on dynamic modeling and have highlighted several
aspects that might be considered for deciding on a modeling approach in a
given application scenario. In particular, our intention has been to
illustrate how ideas can be combined across communities, as in adapting ODE-
based modeling for longitudinal cohort data. Specifically, we were motivated
by an SMA rare disease registry as a prototypical example in a clinical
setting, where underlying disease dynamics are to be inferred from a larger
number of observed variables and modeled in a latent space, using individual-
specific models. There, a particular challenge is that observations are noisy,
irregular in time, and heterogeneous.

Regression approaches typically correspond to function fitting techniques,
where the starting point for subsequent assessments is an extrapolation of an
average fit of all previously observed data on an absolute level, within and
often across individuals, whereas in ODE-based models, relative changes to a
starting point are modeled. We have argued that this difference is not merely
technical, but reflects different analysis intentions, relating to what is
considered the closest approximation of an underlying truth. This is closely
linked to the conceptualization of variability. Function fitting approaches
implicitly or explicitly assume measurements as noisy, such that using an
average increases robustness and is considered more accurate, whereas ODE
approaches using the last observed value as starting point consider this last
observed value as most representative and are thus more sensitive to noise.
While challenging when dealing with noisy data, as frequently encountered in
biomedical applications, predicting future developments conditional on a
patient’s current status locally in time with such an ODE approach is
attractive, as it might more closely match the perspective of a clinician who
seeks to predict immediate changes in the near future based on a patient’s
current status, rather than a global average. We have argued that the local
perspective in ODE-based approaches further allows for modeling each patient’s
time series individually based on external information. This is attractive in
particular when heterogeneity of data can be explained by baseline
characteristics, and in settings where not many follow-up time points are
available yet.

To nonetheless incorporate the advantages of statistical function fitting
approaches, e.g., for dealing with noise, we have thus proposed a statistical
approach based on ODEs, where we decrease dependence on the initial condition
by combining multiple ODE solutions into an inverse-variance weighted unbiased
estimator of the underlying dynamics. To allow for personalized trajectories,
we have used patients’ baseline characteristics to infer individual-specific
ODE parameters. For modeling with a larger number of variables, we have
combined our approach with dimension reduction by a VAE, which allowed for
simultaneous fitting of all model components.

In a simulation design and an application on SMA patient’s motor function
development data, we have shown that the approach allows for inferring
individual-specific trajectories in the VAE latent space and for predicting
subsequent observations. Notably, the approach provides personalized
predictions on a patient’s development immediately upon study entry, as their
baseline measurements and first measurement of the time-dependent variables
allow for fully specifying an individual-specific ODE system. Joint
optimization of the VAE and the dynamic model, facilitated by differentiable
programming, allows for adapting the latent representation to the dynamic
model, which acts as a regularizing smoothness constraint.

Yet, the latent representation is only identifiable up to an affine linear
transformation, due to the encoder neural network structure. While the dynamic
model parameters can be interpreted in the context of the latent
representation, it might be challenging to link this interpretation to the
observed data. For facilitating clinical insight, the latent representation
could be constrained more strongly or combined with a post-hoc explainability
technique. Alternatively, latent trajectories can be decoded back to
observation space and the data-level trajectories can be queried for clinical
patterns. More generally, interpreting the latent representation and choosing
a suitable ODE system is challenging. In the present application, we used a
simple linear ODE system with two components to reflect the underlying changes
in motor neuron functionality and muscle fitness as potential driving forces,
and also preferred a model with few parameters as the number of observations
was rather small. Alternatively, specification of the ODE system could be
addressed by a model selection strategy, e.g., by comparing prediction
performance of models with different degrees of complexity in a stepwise
approach. A further current limitation of the proposed approach is that time
points of follow-up visits are assumed to be independent of the latent state.
It might be interesting to relax this assumption in future work. Then,
prediction uncertainty could be used as a criterion for determining when a
patient should ideally be seen again.

In summary, the proposed approach provides a flexible modeling alternative for
assessing individual-specific dynamics in longitudinal cohort data, as
illustrated in a prototypical application scenario, and more generally
exemplifies the benefits of integrating advantages of different modeling
strategies across communities.

##  6 SUPPLEMENTARY MATERIAL

###  6.1 Implementation details

We use differentiable programming for simultaneous optimization of the dynamic
model and the VAE encoder and decoder. Differentiable programming (Innes,
2019) is a new paradigm for flexibly coupling optimization of diverse model
components. For an overview, see our own work Hackenberg et al. (2022a).
Briefly, differentiable programming facilitates a combination of different
modeling components, such as neural network and differential equations, by
allowing for joint optimization of an overall loss function using automatic
differentiation (Baydin et al., 2017). Thus, gradients with respect to all
model components can be obtained and optimized simultaneously. While
originally proposed to integrate scientific computing and machine learning
(Innes et al., 2019; Rackauckas et al., 2020), differentiable programming can
also be useful for flexible statistical modeling (Hackenberg et al., 2022a)
and has been combined with statistical approaches, e.g., in an adaptation of
neural ODEs for a Bayesian framework (Dandekar et al., 2021), or combinations
of ODEs and non-linear mixed effects models for pharmacometric modeling
(Rackauckas et al., 2022).

In our approach, we use the flexible automatic differentiation framework
provided by Zygote.jl. (Innes et al., 2019) for jointly optimizing the dynamic
model and the VAE for dimension reduction, such that the two components can
influence each other. In this way, a latent representation can be found that
is automatically adapted to the underlying dynamics and the ODE system
structures and regularizes the representation. As the ODEs have analytical
solutions, differentiation through the latent dynamics estimator does not
require backpropagating gradients through a numerical ODE solving step.
However, differentiable programming also allows for doing that efficiently if
necessary, e.g., using the adjoint sensitivity method Chen et al. (2018).

The model is implemented in the Julia programming language (Bezanson et al.,
2017) of version 1.7.2. with the additional packages CSV.jl (v0.10.9),
DataFrames.jl (v1.2.2), Distributions.jl (v0.25.24), Flux.jl (v0.12.8) (Innes,
2018), GLM.jl (v1.5.1), LaTeXStrings.jl (v1.3.0), MultivariateStats.jl
(v0.8.0), Parameters.jl (v0.12.3), Plots.jl (v1.23.5) and StatsBase.jl
(v0.33.21). The model is available as a Julia package at
https://github.com/maren-ha/LatentDynamics.jl. The repository also provides
Jupyter notebooks for illustrating usage of the approach and reproducing the
results in the manuscript.

In the VAE, the encoder and decoder each have one hidden layer with the number
of hidden units equal to the number of input dimensions and a
$\tanh$-activation functions. The latent space is two-dimensional, and the
latent space mean and variance are obtained as affine linear transformations
of the hidden layer values without a non-linear activation. Similarly, the
decoder parameterizes mean and variance of a Gaussian distribution calculated
from the decoder hidden layer using an affine linear transformation.

The network mapping baseline variables to ODE parameters has two hidden layers
in addition to input and output layer. In the first hidden layer, the number
of hidden units corresponds to the number of input variables and a
$\tanh$-activation is used. In the second hidden layer, the number of hidden
units corresponds to the number of ODE parameters and the activation function
is a sigmoid function, shifted by $-0.5$ and scaled by $0.5$ (resp. $1$ in the
simulation design). This effectively acts as a prior for the range of the
estimated ODE parameters. Deviations from this range are possible, as an
affine linear transformation with a diagonal matrix is added as the final
output layer.

In the loss function, the KL-divergence between the prior and posterior is
scaled by a factor of $0.5$ to slightly reduce the regularizing effect of the
$\mathcal{N}(0,1)$ prior. The sum of the squared decoder parameters is added
with a weighting factor of $0.01$ to the loss function, corresponding to a
commonly used penalty term to prevent exploding decoder parameters.

Additionally, the variance penalty term described in Section 3.2.2 is added to
the loss function with a weight of $2$ in the SMArtCARE example. This is not
necessary for the simulation, as the variance is already controlled by the
data generating process. The squared Euclidean distance between the posterior
means from the encoded time series and from the individual ODE solutions for
consistency of the latent representation before and after solving the ODEs is
added with a weight of $1$ in the SMArtCARE application and a weight of $0.5$
in the simulation. In the notation of Equation (6), we thus have
$\alpha=1,\beta=2$ for the SMArtCARE example and $\alpha=0.5,\beta=0$ in the
simulation.

The models are trained based on stochastic gradient descent using the ADAM
optimizer (Kingma & Ba, 2015) with a learning rate of $0.0005$ and $30$
training epochs, chosen by monitoring convergence of the loss function and
visualizing exemplary fits.

###  6.2 Derivation of the analytical solution to the linear ODE system

In the following, we show how to arrive at the analytical solution in Equation
(2) to the initial value problem from Equation (1). For
$A\in\mathbb{R}^{d\times d},A\neq 0$ and $c\in\mathbb{R}^{d}$, we have

$\displaystyle z^{\prime}(t)$ $\displaystyle=A\cdot z(t)+c;$ $\displaystyle
z(t_{0})$ $\displaystyle=z_{0}.$

Subtracting $A\cdot z(t)$ and multiplying by $\exp(-At)$ yields

$\displaystyle(z^{\prime}(t)-Az(t))\exp(-At)=c\cdot\exp(-At).$ (7)

According to the product rule,

$\frac{d}{dt}(z(t)\exp(-At))=z^{\prime}(t)\exp(-At)-Az(t)\exp(-At)=(z^{\prime}(t)-Az(t))\exp(-At),$

and hence, substituting this identity in Equation (7) we obtain

$\frac{d}{dt}(z(t)\exp(-At))=c\exp(-At).$

Integrating both sides yields

$\displaystyle\begin{split}z(t)\exp(-At)&=\int_{t_{0}}^{t}c\exp(-As)\mathrm{d}s+K\text{~{}~{}(for
some constant }K\in\mathbb{R}^{d})\\\
&=-A^{-1}c\exp(-At)+A^{-1}c\exp(-At_{0})+K,\end{split}$

and thus

$\displaystyle\begin{split}z(t)&=\left(-A^{-1}c\exp(-At)+A^{-1}c\exp(-At_{0})+K\right)\exp(At)\\\
&=(A^{-1}c\exp(-At_{0})+K)\exp(At)-A^{-1}c.\end{split}$ (8)

Evaluating Equation (8) at $t=t_{0}$, we obtain

$\displaystyle\begin{split}z(t_{0})&=(A^{-1}c\exp(-At_{0})+K)\exp(At_{0})-A^{-1}c\\\
&=A^{-1}c+K\exp(At_{0})-A^{-1}c\\\ &=K\exp(At_{0}),\end{split}$

and thus $K=z(t_{0})\exp(-At_{0})=z_{0}\exp(-At_{0})$. Inserting this into
Equation (8), we finally arrive at the analytical solution from Equation (2):

$\displaystyle z(t)$
$\displaystyle=(A^{-1}c\exp(-At_{0})+z_{0}\exp(-At_{0}))\exp(At)-A^{-1}c$
$\displaystyle $ $\displaystyle=(A^{-1}c+z_{0})\exp(A(t-t_{0}))-A^{-1}c.$

###  6.3 Simulation design with a simpler ODE system

In addition to the simulation using a homogeneous linear ODE with four unknown
parameters from Section 4.1, we performed a simulation using a simpler two-
dimensional system with just two unknown parameters. As before we defined two
groups of individuals with distinct underlying development patterns. The
corresponding ODE systems are given by

$\begin{split}\frac{d}{dt}\begin{pmatrix}u_{1}\\\
u_{2}\end{pmatrix}(t)&=\begin{pmatrix}-0.2\\\ 0.1\end{pmatrix};\\\
\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}(0)&=\begin{pmatrix}2\\\
1\end{pmatrix}\end{split}$

$\begin{split}\frac{d}{dt}\begin{pmatrix}u_{1}\\\
u_{2}\end{pmatrix}(t)&=\begin{pmatrix}-0.1\\\ -0.2\end{pmatrix};\\\
\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}(0)&=\begin{pmatrix}2\\\
1\end{pmatrix}.\end{split}$

We simulated data for $100$ individuals and trained our model on the data as
described in Section 4.1. Exemplary fits from $12$ simulated individuals are
shown in Figure 4. Additionally, we show the fits of piecewise and global
linear regression fits in Figures 5 and 5, as described in Section 4.1.

For the global regression model in Figure 6, we have calculated a regression
fit using the complete observed time series of each patient, to empirically
verify that our approach indeed matches a simple least-squared fit when
underlying dynamics are simple.

For the piecewise regression model in Figure 5, we have calculated
extrapolations based on previously observed values only, as discussed in 4.1.
Naturally, this more limited observation leads to a drop in prediction
performance, both when comparing predictions on subsequent time points in
latent space and when decoding the respective predictions back to the level of
the original items using the VAE decoder, as summarized below in Table 1.

Model | Prediction in latent space | Prediction on reconstructed values
---|---|---
ODE | 0.406 | 1.129
Regression | 0.981 | 3.686
Table 1: Mean squared error averaged across all time points and all simulated
patients for predicting the subsequent time point both in latent space and on
the level of the reconstructed items, using the ODE approach and a linear
regression, fit on all previously observed data points. Figure 4: (a) Ground-
truth trajectories for two distinct development patterns based on which data
was simulated. (b) Exemplary fits of simulated individuals from both groups.
One panel corresponds to one individual, the $x$ axis shows the time and the
$y$ axis the value of the latent representation. Circles depict the mean of
the latent variable as obtained directly from the encoder, while solid lines
show the ODE solution with parameters obtained from the individual’s baseline
variables, using the value at the last observed time point as initial
condition and solved until the subsequent time point. Colored bands show the
range of one standard deviation of the latent variable around its mean. The
two leftmost columns show $6$ individuals from the first group corresponding
to the left ground-truth trajectory in (a), while the two rightmost columns
show $6$ individuals from the second group, corresponding to the pattern shown
on the right in (a).  Figure 5: Exemplary fits of simulated individuals from
Figure 4, with a piecewise linear regression fit shown in pink. For each
simulated individual, at each time point in each latent dimension, a linear
regression fit of all latent space values from previous time points of that
individual and dimension was calculated. The regression line was shifted to
the value of the last observation used for calculation, and extrapolated until
the subsequent time point (corresponding to the adaptation of an initial value
approach for a regression model discussed in Section 2). Figure 6: Exemplary
fits of simulated individuals from Figure 4, with global linear regression
fits shown in pink. For each simulated individual, in each latent dimension a
linear regression fit of the complete time series of that individual and
dimension was calculated (corresponding to the classical function fitting
approach for regression models discussed in Section 2).

## References

  * (1)
  * Bandyopadhyay et al. (2011) Bandyopadhyay, S., Ganguli, B. & Chatterjee, A. (2011), ‘A review of multivariate longitudinal data analysis’, Statistical Methods in Medical Research 20(4), 299–330.  PMID: 20212072.
  * Banga (2008) Banga, J. R. (2008), ‘Optimization in computational systems biology’, BMC Systems Biology 2(1), 1–7.
  * Baydin et al. (2017) Baydin, A. G., Pearlmutter, B. A., Radul, A. A. & Siskind, J. M. (2017), ‘Automatic differentiation in machine learning: A survey’, J. Mach. Learn. Res. 18(1), 5595–5637.
  * Bezanson et al. (2017) Bezanson, J., Edelman, A., Karpinski, S. & Shah, V. B. (2017), ‘Julia: A fresh approach to numerical computing’, SIAM Review 59(1), 65–98.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A. & McAuliffe, J. D. (2017), ‘Variational inference: A review for statisticians’, Journal of the American Statistical Association 112(518), 859–877.
  * Bruggeman & Westerhoff (2007) Bruggeman, F. J. & Westerhoff, H. V. (2007), ‘The nature of systems biology’, Trends in Microbiology 15(1), 45–50.
  * Chen et al. (2018) Chen, T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. (2018), Neural ordinary differential equations, in S. Bengio, H. M. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi & R. Garnett, eds, ‘Advances in Neural Information Processing Systems 31’, pp. 6572–6583.
  * Dandekar et al. (2021) Dandekar, R., Chung, K., Dixit, V., Tarek, M., Garcia-Valadez, A., Vemula, K. V. & Rackauckas, C. (2021), ‘Bayesian neural ordinary differential equations’.  arXiv preprint: https://arxiv.org/abs/2012.07244.
  * de Haan-Rietdijk et al. (2017) de Haan-Rietdijk, S., Voelkle, M. C., Keijsers, L. & Hamaker, E. L. (2017), ‘Discrete-vs. continuous-time modeling of unequally spaced experience sampling method data’, Frontiers in psychology 8, 1849.
  * DiStefano (2015) DiStefano, J. (2015), Dynamic systems biology modeling and simulation, Academic Press.
  * Fan et al. (2007) Fan, J., Huang, T. & Li, R. (2007), ‘Analysis of longitudinal data with semiparametric estimation of covariance function’, Journal of the American Statistical Association 102(478), 632–641.
  * Fitzmaurice et al. (2008) Fitzmaurice, G., Davidian, M., Verbeke, G. & Molenberghs, G. (2008), Longitudinal data analysis, CRC press.
  * Graves & Graves (2012) Graves, A. & Graves, A. (2012), ‘Long short-term memory’, Supervised sequence labelling with recurrent neural networks pp. 37–45.
  * Gwak et al. (2020) Gwak, D., Sim, G., Poli, M., Massaroli, S., Choo, J. & Choi, E. (2020), ‘Neural ordinary differential equations for intervention modeling’.  arXiv preprint: https://arxiv.org/abs/2010.08304.
  * Gábor & Banga (2015) Gábor, A. & Banga, J. R. (2015), ‘Robust and efficient parameter estimation in dynamic models of biological systems’, 9(1), 74.
  * Hackenberg et al. (2022a) Hackenberg, M., Grodd, M., Kreutz, C., Fischer, M., Esins, J., Grabenhenrich, L., Karagiannidis, C. & Binder, H. (2022a), ‘Using differentiable programming for flexible statistical modeling’, The American Statistician 76(3), 270–279.
  * Hackenberg et al. (2022b) Hackenberg, M., Harms, P., Pfaffenlehner, M., Pechmann, A., Kirschner, J., Schmidt, T. & Binder, H. (2022b), ‘Deep dynamic modeling with just two time points: Can we still allow for individual trajectories?’, Biometrical Journal .
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. & Schmidhuber, J. (1997), ‘Long short-term memory’, Neural computation 9(8), 1735–1780.
  * Innes (2018) Innes, M. (2018), ‘Flux: Elegant machine learning with Julia’, Journal of Open Source Software 3(25), 602.
  * Innes (2019) Innes, M. (2019), ‘What is differentiable programming?’.  Blogpost available at https://fluxml.ai/blogposts/2019-02-07-what-is-differentiable-programming/ (date of retrieval: December 19, 2022).
  * Innes et al. (2019) Innes, M., Edelman, A., Fischer, K., Rackauckas, C., Saba, E., Shah, V. B. & Tebbutt, W. (2019), ‘A differentiable programming system to bridge machine learning and scientific computing’.  arXiv preprint: https://arxiv.org/abs/1907.07587.
  * Kingma & Ba (2015) Kingma, D. P. & Ba, J. (2015), Adam: A method for stochastic optimization, in Y. Bengio & Y. LeCun, eds, ‘3rd International Conference on Learning Representations (ICLR), Conference Track Proceedings’.
  * Kingma & Welling (2014) Kingma, D. P. & Welling, M. (2014), Auto-encoding variational Bayes, in Y. Bengio & Y. LeCun, eds, ‘2nd International Conference on Learning Representations (ICLR), Conference Track Proceedings’.
  * Kingma & Welling (2019) Kingma, D. P. & Welling, M. (2019), ‘An introduction to variational autoencoders’, Foundations and Trends® in Machine Learning 12(4), 307–392.
  * Köber et al. (2022b) Köber, G., Kalisch, R., Puhlmann, L., Chmitorz, A., Schick, A. & Binder, H. (2022b), ‘Deep learning and differential equations for modeling changes in individual-level latent dynamics between observation periods’, arXiv preprint arXiv:2202.07403 .
  * Köber et al. (2022a) Köber, G., Pooseh, S., Engen, H., Chmitorz, A., Kampa, M., Schick, A., Sebastian, A., Tüscher, O., Wessa, M., Yuen, K. S. et al. (2022a), ‘Individualizing deep dynamic models for psychological resilience data’, Scientific Reports 12(1), 1–10.
  * Kreutz et al. (2012) Kreutz, C., Raue, A. & Timmer, J. (2012), ‘Likelihood based observability analysis and confidence intervals for predictions of dynamic models.’, BMC System Biology 6, 120.
  * Liang & Wu (2008) Liang, H. & Wu, H. (2008), ‘Parameter estimation for differential equation models using a framework of measurement error in regression models’, Journal of the American Statistical Association 103(484), 1570–1583.
  * Mazzone et al. (2017) Mazzone, E. S., Mayhew, A., Montes, J., Ramsey, D., Fanelli, L., Young, S. D., Salazar, R., Sanctis, R. D., Pasternak, A., Glanzman, A., Coratti, G., Civitello, M., Forcina, N., Gee, R., Duong, T., Pane, M., Scoto, M., Pera, M. C., Messina, S., Tennekoon, G., Day, J. W., Darras, B. T., Vivo, D. C. D., Finkel, R., Muntoni, F. & Mercuri, E. (2017), ‘Revised upper limb module for spinal muscular atrophy: Development of a new module’, Muscle Nerve 55(6), 869–874.
  * Mews et al. (2022) Mews, S., Langrock, R., Ötting, M., Yaqine, H. & Reinecke, J. (2022), ‘Maximum approximate likelihood estimation of general continuous-time state-space models’, Statistical Modelling .
  * Moon et al. (2022) Moon, I., Groha, S. & Gusev, A. (2022), Survlatent ode : A neural ode based time-to-event model with competing risks for longitudinal data improves cancer-associated venous thromboembolism (vte) prediction, in ‘Machine Learning for Healthcare’, Vol. 182 of Proceedings of Machine Learning Research, p. 1–28.
  * Nolan et al. (2023) Nolan, T. H., Goldsmith, J. & Ruppert, D. (2023), ‘Bayesian Functional Principal Components Analysis via Variational Message Passing with Multilevel Extensions’, Bayesian Analysis -1(-1), 1–27.
  * Palma et al. (2023) Palma, M., Keogh, R. H., Carr, S. B., Szczesniak, R., Taylor-Robinson, D., Wood, A. M., Muniz-Terrera, G. & Barrett, J. K. (2023), ‘Factors associated with within-individual variability of lung function for people with cystic fibrosis: a longitudinal registry study’.
https://www.medrxiv.org/content/10.1101/2023.05.12.23289768v1

  * Pechmann et al. (2019) Pechmann, A., König, K., Bernert, G., Schachtrup, K., Schara, U., Schorling, D., Schwersenz, I., Stein, S., Tassoni, A., Vogt, S., Walter, M. C., Lochmüller, H. & Kirschner, J. (2019), ‘Smartcare - a platform to collect real-life outcome data of patients with spinal muscular atrophy’, Orphanet Journal of Rare Diseases 14(18).
  * Rackauckas et al. (2020) Rackauckas, C., Ma, Y., Martensen, J., Warner, C., Zubov, K., Supekar, R., Skinner, D. & Ramadhan, A. (2020), ‘Universal differential equations for scientific machine learning’.  arXiv preprint: https://arxiv.org/abs/2001.04385.
  * Rackauckas et al. (2022) Rackauckas, C., Ma, Y., Noack, A., Dixit, V., Mogensen, P. K., Elrod, C., Tarek, M., Byrne, S., Maddhashiya, S., Calderón, J. S., Hatherly, M., Nyberg, J., Gobburu, J. V. S. & Ivaturi, V. (2022), ‘Accelerated predictive healthcare analytics with pumas, a high performance pharmaceutical modeling and simulation platform’, bioRxiv p. 2020.11.28.402297.
  * Ramsay & Hooker (2017) Ramsay, J. & Hooker, G. (2017), Dynamic data analysis, Springer.
  * Raue et al. (2009) Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmüller, U. & Timmer, J. (2009), ‘Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood’, Bioinformatics 25(15), 1923–1929.
  * Rubanova et al. (2019) Rubanova, Y., Chen, T. Q. & Duvenaud, D. (2019), ‘Latent ODEs for irregularly-sampled time series’.  arXiv preprint: https://arxiv.org/abs/1907.03907.
  * Schiratti et al. (2017) Schiratti, J.-B., Allassonnière, S., Colliot, O. & Durrleman, S. (2017), ‘A Bayesian Mixed-Effects Model to Learn Trajectories of Changes from Repeated Manifold-Valued Observations’, Journal of Machine Learning Research 18(133), 1–33.
  * Schorling et al. (2020) Schorling, D. C., Pechmann, A. & Kirschner, J. (2020), ‘Advances in treatment of spinal muscular atrophy - new phenotypes, new challenges, new implications for care’, Journal of neuromuscular diseases 7(1), 1–13.
  * Shahar (2017) Shahar, D. J. (2017), ‘Minimizing the variance of a weighted average’, Open Journal of Statistics 7(2), 216–224.
  * Van Montfort et al. (2018) Van Montfort, K., Oud, J. H. & Voelkle, M. C. (2018), Continuous time modeling in the behavioral and related sciences, Springer.
  * Verbeke et al. (2014) Verbeke, G., Fieuws, S., Molenberghs, G. & Davidian, M. (2014), ‘The analysis of multivariate longitudinal data: a review.’, Stat Methods Med Res 23(1), 42–59.
  * Xue et al. (2010) Xue, H., Miao, H. & Wu, H. (2010), ‘Sieve estimation of constant and time-varying coefficients in nonlinear ordinary differential equation models by considering both numerical error and measurement error.’, Ann Stat 38(4), 2351–2387.