# Curve registration for non-linear mixed effect ordinary differential equation models

Quentin Clairon, John Fricks, Mélanie Prague

► To cite this version:

Quentin Clairon, John Fricks, Mélanie Prague. Curve registration for non-linear mixed effect ordinary differential equation models. 2025. hal-05412063

HAL Id: hal-05412063

https://inria.hal.science/hal-05412063v1

Preprint submitted on 11 Dec 2025

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 Attribution 4.0 International License
# Curve registration for non-linear mixed effect ordinary differential equation models

Quentin Clairon

Université de Bordeaux, Inserm, Inria Bordeaux Sud-Ouest, SISTM team, Bordeaux, France, Vaccine research Institute, University of Paris-Est, Créteil, France

John Fricks

School of Mathematical and Statistical Sciences Arizona State University, Tempe, USA

Mélanie Prague

Université de Bordeaux, Inserm, Inria Bordeaux Sud-Ouest, SISTM team, Bordeaux, France, Vaccine research Institute, University of Paris-Est, Créteil, France

## Abstract

Curve registration refers to the alignment of temporal trajectories when individuals experience similar dynamics at different speeds. We study this problem in the context of nonlinear mixed-effects models based on ordinary differential equations (NLME-ODEs), which are commonly used to describe biological mechanisms in longitudinal data. We extend NLME-ODEs by introducing subject and state specific time-warping functions that represent local accelerations or decelerations of the dynamics. The goal is to propose mechanistic longitudinal population data analysis which incorporates and locates model error presence to mitigate its effect on model parameter estimation and acknowledge the action of mechanisms ignored by the ODE structure. We then focus on the treatment effect inference problem where these warping functions help to locate as accurately and exhaustively as possible the effects of a given therapy from a followed treated population. From the original model, we quantify treatment effects on well-identified mechanisms specified as ODE parameter covariates, now robustly estimated in model error presence. From the time-warping functions, we infer its action on mechanisms not included in the assumed ODE but appearing as distorted dynamics for treated subjects compared to the untreated ones. We formulate the joint estimation of ODE parameters and warping functions in a mixed-
effects likelihood framework and implement it in standard softwares such as Monolix implementing SAEM algorithm. Simulation studies show that including time warping improves estimation accuracy. Application to simian immunodeficiency virus (SIV) pre-clinical trial data demonstrates that the method can detect treatment effects on both modeled mechanisms and unmodeled processes. This framework combines mechanistic interpretability with descriptive flexibility, enhancing robustness in the analysis of longitudinal trial data.

*Keywords:* Mechanistic model inference, Functional data analysis, Nonlinear mixed-effects models, Longitudinal data analysis, Population approach

# 1. Introduction

We focus on the problem of longitudinal population data analysis from a mechanistic perspective, that is based on the specification of an ordinary differential equation (ODE) model. More particularly, we propose an approach which better handles model misspecification compared to classic inference procedures by providing parameter estimation more robust to its presence and by locating where (i.e. on which state variable dynamic) this model error significantly occurs. We aim to do so via the incorporation of techniques coming from functional data analysis (FDA), notably curve registration procedures, into the classic mechanistic setting.

ODEs are widely used in pharmacometric/pharmacodynamic, viral dynamic analysis and also vaccinology [1, 2, 3, 4, 5]. Their success can be explained by three important points. First, their interpretability, such models describe causal links between variables of interest in which the parameters are understood as the mathematical translation of well-identified biological mechanisms [6]. Second, these models provide predictive ability for temporal forecasting far outside the restricted available observation interval. Third, their ability to reconstruct unobserved variables from the causal links they prescribed between them and the available measurements. This last point makes the mechanistic framework particularly suited for data analysis in a partially observed setting.

More formally, we consider a population of *n* subjects in which the time-varying state of the i-th one is represented by the *D* interacting variables
24 $X_i = (X_i^{(1)}, ..., X_i^{(D)})$, solution of the ODE:

$$ \left\{ \begin{aligned} \dot{X}_i &= f_{\xi_i}(X_i, t) \\ X_i(0) &= x_{i,0}(\xi_i), \end{aligned} \right. \qquad (1.1) $$

25 depending on the $D_\xi$-dimensional parameter $\xi_i$. Then, the observation model:

$$ Y_{ij} = C(X_i(t_{i,j})) + \epsilon_{i,j} \quad \text{with} \quad \epsilon_{i,j} \sim N(0, \Sigma), \qquad (1.2) $$

26 is associated to ODE (1.1) to link, through the observation function $C$:
27 $\mathbb{R}^D \to \mathbb{R}^{D_{obs}}$, the i-th subject state to the partial, noisy and discrete observations $\{Y_{ij}\}_{j=1,...,n_i}$ available for it and made at the $n_i$ timepoints $\{t_{i,j}\}_{j=1,...,n_i}$.
28 Additionally, because of the observational frameworks imposed by most population settings, $\xi_i$ has to be specified in such a way that 1/it can be estimated
29 from sparse data at the subject level while acknowledging 2/the existence of
30 subgroups within the population stratified according to specific features such
31 as age, sex, exposure or treatment and 3/account for the presence of intrinsic
32 inter-subject variability. This is done by adopting a nonlinear mixed-effect
33 (NLME) setting, leading to NLME-ODEs, where $\xi_i$ is decomposed as:
34
35

$$ \xi_{i,l} = \bar{\xi}_l + \sum_{m=1}^{M} \pi_{G_m}^{\xi_l} \mathbb{I}_{\{i \in G_m\}} + b_i^{\xi_l} \quad \text{for } l = 1, \dots, D_\xi, \qquad (1.3) $$

36 with:

37 $\bar{\xi}_l$: the mean value for $\xi_{i,l}$ common to the whole population,

38 $\pi_{G_m}^{\xi_l}$: the covariate term applied to $\bar{\xi}_l$ for subjects belonging to the group $G_m$ (with possibly $\pi_{G_m}^{\xi_l} = 0$),
39

40 $b_i^{\xi_l} \sim N(0, \sigma_{\xi_l}^2)$: the subject-specific variations (with possibly $\sigma_{\xi_l}^2 = 0$).

41 Indeed, the parametrization (1.3) makes it possible to aggregate information from the whole population to estimate $\bar{\xi} := (\bar{\xi}_1, ..., \bar{\xi}_{D_\xi})$ and $\pi^\xi := (\pi_{G_m}^{\xi_1}, ..., \pi_{G_m}^{\xi_{D_\xi}})_{m=1,...,M}$ to counteract the sparsity of available measurement
42 per subject while allowing for heterogeneity through the inclusion of $b_i^\xi := (b_i^{\xi_1}, ..., b_i^{\xi_{D_\xi}}) \sim N(0, \Sigma_{b^\xi} = \text{diag}(\sigma_{\xi_l}^2))$.

43 The point of such an analysis is generally to enlighten and quantify the
44 existing differences regarding the temporal evolution of some variables of interest between groups $\{G_m\}_{m=1,...,M}$ presenting altered features compared to
45
49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

a base one, taken by convention as $G_0$. In this setting, this question is turned into the problem of estimating $(\bar{\xi}, \pi^{\xi}, \Sigma, \Sigma_{b^{\xi}})$ with frequentist or Bayesian approaches which are classical in their core principles but still adapted to the computational specificity of NLME-ODE inference [7, 8, 9]. This NLME-ODE formalism indeed automatically provides interpretable analyses of the existing group differences in terms of altered biological mechanisms through $\pi^{\xi}$ inference. Moreover, the effects of these altered features are allowed to be indirect. They can primarily act on unmeasured variables which then interact with the observed ones through the causal links between variables prescribed by the structure of $f_{\xi_i}$.

This approach suffers from two critical issues, both coming from the limited precision of the ODE (1.1) constructed from significant theoretical simplification of the real biological processes at play, both in terms of represented variables and interactions between them. First, it limits the number of mechanisms for which we can observe group induced alterations, potentially leading to an underestimation of the existing differences. Second, model error causes a loss of estimation accuracy for inference methods which does not account for its presence [10, 11, 12], such as in the classically employed likelihood or Bayesian inference methods [7, 8, 9]. This is critical for $\pi^{\xi}$ estimation which may lead to false conclusion about group effect on a given mechanism.

To the contrary, many descriptive model-free procedures are not concerned with model error related pitfalls. Among them, some devices coming from FDA, which consider functions as the basic unit data [13], may be particularly suited for our purpose, since we analyze group effects on time-varying variables $Y_i : t \mapsto C(X_i(t))$. Among the different areas tackled by FDA such as functional regression [14], functional principal components analysis [15], or clustering [16], we frame our problem within the curve registration setting [17, 18]. That is, we aim to infer time-warping functions $\{h_i\}_{i=1,...,n}$ distorting the temporal evolution of the regression functions $\{t \mapsto Y_i(h_i(t))\}_{i=1,...,n}$ from noisy and discrete observations $Y_i(h_i(t_{i,j}))+\epsilon_{i,j}$. The point being to make relevant comparisons between the properly re-aligned curves $\{Y_i\}_{i=1,...,n}$ but also inspect the kinetic distortions they have been affected by through the direct analysis of the functions $\{h_i\}_{i=1,...,n}$. In our case, we can use curve registration techniques to precisely quantify the differences in terms of magnitude between two groups $G_m$ and $G_{m'}$ by comparing the re-aligned outcomes $\{Y_i\}_{i \in G_m}$ and $\{Y_i\}_{i \in G_{m'}}$, learned from the observations of their distorted versions, but also in terms of what could be thought of as kinetic speed by focusing on $\{h_i\}_{i \in G_m}$ and $\{h_i\}_{i \in G_{m'}}$ comparisons. Still, avoiding model error
via direct use of FDA devices is done at the expense of explainability and the
possibility to explore the presence of group specific alterations on unobserved
variables dynamics given by the mechanistic perspective.

That is why in this work we integrate curve registration techniques into
NLME-ODEs (1.1)-(1.2) to propose an inference procedure accounting for
model error to 1/mitigate its effect on parameter estimation accuracy and
2/locates on which part of the model this error is present in particular and
how its shape differs from one group to another. By doing so, we can quantify
group specific differences that address limitations encountered by descriptive
approaches and mechanistic ones when used separately. For this, we con-
sider that the state of each subject, represented by $X_i \circ h_i$, has its dynamics
locally altered by temporal distortion, represented by $h_i$, which perturbs its
evolution compared to the one described by the ODE (1.1). Moreover, we
assume that $h_i$ is state-variable specific, i.e. $h_i(t) := (h_i^{(1)}(t), \dots, h_i^{(D)}(t))$, to
acknowledge that the different state variables are not submitted to the same
sources of perturbations, potentially leading to an asynchronous evolution of
the components of $X_i \circ h_i$ comparing to ODE (1.1). So, our perspective is a
bit different from the usual FDA setting, here $h_i$ represents local perturba-
tions acting on the state-variable dynamics and not a temporal shift directly
applied on the resulting motion, as it is classically done with time-warping.

To stay consistent with our population context, we assume the func-
tions {$h_i$}$_{i=1,...,n}$ can be decomposed in a way mirroring $\xi_{i,l}$, that is with
a mean value $\bar{h} := (\bar{h}^{(1)}, \dots, \bar{h}^{(D)})$ common to everyone, group specific dis-
tortions $\pi^h := (\pi_{G_1}^{h^{(1)}}, \dots, \pi_{G_1}^{h^{(D)}}, \dots, \pi_{G_M}^{h^{(1)}}, \dots, \pi_{G_M}^{h^{(D)}})$ and subject-specific random
variations $b_i^h = (b_i^{h^{(1)}}, \dots, b_i^{h^{(D)}})$. In this case, a subsequent analysis leading to
$\pi_{G_M}^{h^{(d)}}(t) \neq 0$ would indicate a distorted dynamic for the d-th state variable
of the subjects belonging to $G_M$ compared to the control group evolution,
represented by $\bar{h}^{(d)}$. This gives us a way to infer alterations external to what
be can accounted for the initial ODE (1.1) in an observation setting known
as sparse functional data analysis in FDA [19, 20]. So, from this merged
descriptive/mechanistic approach, we intend to describe as exhaustively as
possible the effect of group specific features on both 1/known mechanisms
ruling some precise part of the state variables dynamic and 2/the global evo-
lution of some state variables but through mechanisms ignored by the original
model. This also opens the possibility to quantify with descriptive methods
global group specific effects on the dynamics of unobserved state variables
linked to data through the interaction specified by the ODE (1.1).

To define a unified statistical model, we rely on a modified version of the differential equation representation of $h_i$ proposed by Ramsay & Li [17] to construct a global NLME-ODE embedding both $\{\xi_i\}_{i=1,...,n}$ and $\{h_i\}_{i=1,...,n}$ as parameters. We then use finite dimensional approximations of $(\bar{h}, \pi^h, b_i^h)$ to treat $(\xi_i, h_i)$ inference as a classic estimation problem in a mixed-effect framework. We take advantage of this classic parametric setting to jointly estimate them via likelihood based inference procedure by using the SAEM algorithm implemented in Monolix software [7] and quantify the related estimation uncertainty using the Fisher Information matrix.

We then use our procedure for the specific problem of treatment effect analysis where the $M$ groups $\{G_m\}_{m=1,...,M}$ represent different subpopulations submitted to different treatments ($G_0$ being the placebo). We simulate data from a time-warped version of the ODE proposed by Perelson et al. [21] describing the joint evolution of CD4+T cells and HIV viral load. We also investigate the effect of model error by simulating data from a more complex ODE model than the one used for parameter estimation. In both cases, adding time-warping allows us to account for 1/external treatment effects ignored by the original ODE and 2/internal effects on ODE parameters despite model misspecification. We finally proceed to the real data analysis of viral load concentration measurement made on SIV infected non-human primates (NHPs) [22, 23]. These measurements are made after the injections of different treatments, the point being to infer their effects in terms of viral load control and reduction. Our procedure manages to quantify treatment actions on an identified immune response mechanism, but also on the general dynamics of unobserved CD4+T cells through mechanisms ignored by the original model.

In section 2, we present the unified model accounting for both the original ODE (1.1) and time-distortion phenomenon as well as its parametrization which is well-suited for estimation purposes. In section 3, we proceed to numerical experiments, both in the presence and absence of model error. Then, in section 4, we consider the real data analysis problem. We finally conclude in section 5 by discussing some possible extensions of this work.

## 2. Time-warped-NLME-ODEs

In this section, we start with the presentation of a unified model for data whose conditional mean function is determined by a parametrized system of
ODE but with stochastic variations affecting its dynamics due to the inclusion of randomly drawn subject-specific parameters and time-warping functions. We then introduce a specific form for the time-warping functions that will facilitate the application of software typically used to fit NLME-ODEs.

## 2.1. Unified model

We assume the ODE (1.1) gives a simplified description of the dynamics governing the evolution of the i-th subject. To account for the deviations from these dynamics, we consider that at each time $t \in [0, T]$, the current subject state is more accurately described by $X_{h_i, i}(t) := X_i(h_i(t))$ than by the unwarped ODE solution value $X_i(t)$ where $h_i: [0, T] \to [0, +\infty)^{(D)}$ is the time-warping function quantifying the discrepancy between the parameterized ODE model alone and the actual evolution. That is, we consider that the actions of perturbations and/or mechanisms ignored by the ODE (1.1) can be represented by local temporal perturbations making each state variable accelerate or decelerate compared to the assumed dynamic. Additionally, since each state-variable present in the ODE (1.1) might not be submitted to the same amount of model misspecification, we consider state-specific distortion: $X_{h_i, i}(t) = (X_i^{(1)}(h_i^{(1)}(t)), \dots, X_i^{(D)}(h_i^{(D)}(t)))$.

But so far, we have just described at each time the current state of the i-th subject. In order to be useful for inference and predictive purposes, we need to derive the evolution of $t \to X_{h_i, i}(t)$. For this, we focus on the behavior of the d-th perturbed state variable $X_i^{(d)} \circ h_i^{(d)}$ and we make the additional assumption that, at a given time and subject's state value, the ODE is still a relevant representation for the local transition from one state to another, once the effect of temporal distortion has been considered. From this, we derive by differentiation of a composed function:

$$ \frac{d}{dt} (X_i^{(d)} \circ h_i^{(d)}(t)) = \dot{h}_i^{(d)}(t) \frac{dX_i^{(d)}(h_i^{(d)}(t))}{dh_i^{(d)}(t)} \\ \approx \dot{h}_i^{(d)}(t) f_{\xi_i}(X_{h_i, i}(t), h_i^{(d)}(t)) $$

as it is assumed that the current subject state is better approximated by $X_{h_i, i}(t)$ than $X_i(t)$. By gathering all the differential equations for $d = 1, \dots, D$, we derive the initial value problem:

$$ \begin{cases} \dot{X}_{h_i, i}(t) = \dot{h}_i(t) \odot f_{\xi_i}(X_{h_i, i}(t), h_i(t)) \\ X_{h_i, i}(0) = x_i(h_i(0)) \end{cases} $$
followed by the warped function $X_{h_i, i}$ in which $\odot$ represents the component-wise product.

As required in a curve registration setting, we need to ensure $h_i^{(d)}$ is strictly increasing on $[0, T]$ for each $d = 1, \dots, D$. For this, we use the representation proposed by Ramsay & Li [17], expressing the time-warping function as the solution of the second order linear ODE $\ddot{h}_i^{(d)} = w_i^{(d)} \dot{h}_i^{(d)}$ which can be equivalently reformulated as a standard first order one:

$$ \left\{ \begin{aligned} \dot{g}_i^{(d)} &= w_i^{(d)} g_i^{(d)} \\ \dot{h}_i^{(d)} &= g_i^{(d)} \end{aligned} \right. \qquad (2.1) $$

via the introduction of $w_i^{(d)}$. It is enough to have $h_i^{(d)}$ strictly increasing if $\dot{h}_i^{(d)}(0) > 0$ without any conditions on $w_i^{(d)}$. This allows us to turn to a constrained inference problem of $h_i$ into the unconstrained one of $w_i = (w_i^{(1)}, \dots, w_i^{(D)})$. We then gather all equations for $d = 1, \dots, D$ into one to describe the evolution of the time-warped NLME-ODE:

$$ \left\{ \begin{aligned} \dot{X}_{h_i, i} &= g_i \odot f_{\xi_i}(X_{h_i, i}, h_i) \\ \dot{g}_i &= w_i \odot g_i \\ \dot{h}_i &= g_i \\ (X_{h_i, i}(0), g_i(0), h_i(0)) &= (x_i(0), \dot{h}_i(0), 0) \end{aligned} \right. \qquad (2.2) $$

in which we have added the constraint $h_i(0) = 0$ to enforce the absence of instantaneous temporal shift for all subjects.

In equation (2.2), setting $w_i^{(d)} = 0$ and $\dot{h}_i^{(d)}(0) = 1$ leads to $h_i^{(d)}(t) = \dot{h}_i^{(d)}(0)t = t$ and corresponds to the absence of temporal distortion. So, the original NLME-ODE (1.1) is embedded into the set of time-distorted models (2.2). Moreover, the particular case $\tilde{w}_i = w_i^{(d)} = \dots = w_i^{(D)}$ and $\tilde{h}_i(0) = \dot{h}_i^{(d)}(0) = \dots = \dot{h}_i^{(D)}(0)$ is contained as well. In this setting, the ODE (2.2) describes the behavior of $X_{h_i, i}(t) = X_i(\tilde{h}_i(t))$, the registered version of the original ODE solution, in which the state variables evolve in a synchronous way, closer to the classic functional data analysis setting where warping is applied to the observed motion. This is in contrast to our general case where the temporal distortion is locally applied to the law of evolution. This new model is thus a generalization accounting for processes with dynamics ruled by the laws described by the original model (1.1) but in locally accelerated or decelerated ways. This allows us to reach a broader class of possible system
behaviors compared to ODE (1.1), in particular by introducing asynchronous evolution among the variables.

Remark 1. We could have also adopted a differential geometry perspective here. By identifying the set of possible time warping functions as a Lie group $\mathcal{L}$ perturbing/acting on the differential manifold $\mathcal{M}$ induced by ODE (1.1), our problem can be seen as the identification within the quotient space $\mathcal{M}/\mathcal{L}$ of some orbits consistent with the data. Still, this is beyond the scope of this paper. For the sake of interpretability and in view of our sparse data setting, we rather identify the functions $h_i$ as model error and/or covariates in the statistical paradigm classically employed for longitudinal population data analysis. For a more geometrical point of view for FDA see for example [24, 25, 26]. Interestingly in [25], the authors also combine curve registration with estimation of ODE parameters by estimating time varying parameters of a Frenet-Serret ODE to study variations among multiple curves to identify an appropriate mean one. Still, their approach was to fit such curves by minimizing a metric to find the time warping functions and the overall mean contrary to our likelihood/mixed-effect based method.

## 2.2. Specification of $(w_i, \dot{h}_i(0))$

We are then left with the choice of $w_i^{(d)}$ and $\dot{h}_i^{(d)}(0)$ for each $d = 1, \dots, D$. Since there are no particular constraints in terms of interpretability on $w_i^{(d)}$ in our case, we can select any finite basis decomposition: $w_i^{(d)}(t) = \sum_{k=1}^{K-1} \beta_{ik}^{(d)} \phi_k(t)$. For simplicity, we choose:

$$ \phi_k(t) = \begin{cases} \frac{1}{t_{k+1}-t_k} \max(t_{k+1}-t_k, t-t_k) \mathbb{I}_{[t_k, +\infty]}(t) & \text{for } t \le t_K, \\ 0 & \text{for } t > t_K. \end{cases} \quad (2.3) $$

Thus $w_i^{(d)}$ is a piecewise continuous linear function on $[0, t_K[$, with slope $\frac{\beta_{ik}^{(d)}}{t_{k+1}-t_k}$ between $t_k$ and $t_{k+1}$, and equal to 0 after $t_K$ where $\{t_1 = 0, \dots, t_K \le T\}$ are $K$ anchor points increasingly reparted on $[0, T]$.

To properly separate the sources of time-distortion across groups and subjects, we use the following decomposition for $\beta_{ik}^{(d)}$, similar to (1.3), to account for whole population, subgroup, and subject specific time-warping:

$$ \left\{ \begin{array}{l} \beta_{i1}^{(d)} = \overline{\beta_1^{(d)}} + \sum_{m=1}^{M} \pi_{G_m}^{\beta_1^{(d)}} \mathbb{I}_{i \in G_m} + b_i^{\beta_1^{(d)}}, \\ \beta_{ik}^{(d)} = \beta_{i1}^{(d)} + \delta \beta_{ik}^{(d)} \text{ for } k \ge 2, \\ \delta \beta_{ik}^{(d)} = \overline{\delta \beta_k^{(d)}} + \sum_{m=1}^{M} \pi_{G_m}^{\delta \beta_k^{(d)}} \mathbb{I}_{i \in G_m} + b_i^{\delta \beta_k^{(d)}}, \end{array} \right. \quad (2.4) $$
with the parameters:

1. $(\overline{\beta_1^{(d)}}, \overline{\delta\beta_k^{(d)}})$: common to the whole population,

2. $(\pi_{G_m}^{\beta_1^{(d)}}, \pi_{G_m}^{\delta\beta_k^{(d)}})$: applied to subjects belonging to $G_m$,

3. $b_i^{\beta_1^{(d)}} \sim N(0, \sigma_{\beta_1^{(d)}}^2)$ and $b_i^{\delta\beta_k^{(d)}} \sim N(0, \sigma_{\delta\beta_k^{(d)}}^2)$: the subject-specific ones.

In equation (2.4), we choose the slope variation $\delta\beta_{ik}^{(d)}$ as the parameter to estimate rather than $\beta_{ik}^{(d)}$ for $k \ge 2$ for the sake of parsimony in the sparse data setting. By doing so, we can set $\pi_{G_m}^{\delta\beta_k^{(d)}} = 0$ and $b_i^{\delta\beta_k^{(d)}} = 0$ and still propagate group and subject-specific variations coming from $(\pi_{G_m}^{\beta_1^{(d)}}, b_i^{\beta_1^{(d)}})$ to $\beta_{ik}^{(d)} = \beta_{i1}^{(d)} + \delta\beta_{ik}^{(d)}$ and thus to $w_i^{(d)}(t)$ for all $t \ge 0$. Moreover, the parameter decomposition (2.4) can be transposed to $w_i^{(d)}$:

$$w_i^{(d)}(t) = \overline{w^{(d)}}(t) + \sum_{m=1}^{M} \pi_{G_m}^{w^{(d)}}(t) \mathbb{I}_{i \in G_m} + b_i^{w^{(d)}}(t) \text{ for } t \le t_K$$

emphasizing the different cause of time-distortions for the d-th state variable, with

1. $\overline{w^{(d)}}(t) = \sum_{k=1}^{K-1} (\overline{\beta_1^{(d)}} + \overline{\delta\beta_k^{(d)}}) \phi_k(t)$: the ones due to model misspecification at population scale,

2. $\pi_{G_m}^{w^{(d)}}(t) = \sum_{k=1}^{K-1} (\pi_{G_m}^{\beta_1^{(d)}} + \pi_{G_m}^{\delta\beta_k^{(d)}}) \phi_k(t)$: the additional ones appearing for the m-th group comparing to the control one,

3. $b_i^{w^{(d)}}(t) = \sum_{k=1}^{K-1} (b_i^{\beta_1^{(d)}} + b_i^{\delta\beta_k^{(d)}}) \phi_k(t)$: the subject-specific ones.

From the decomposition of $w_i^{(d)}$, we can account for the presence of model error differently depending on the scale of interest in the model. Adopting the terminology of Kennedy and O'Hagan [27], $\overline{w^{(d)}}$ quantifies model inadequacy, the discrepancy at the population level between the mean model response and real world process, whereas $b_i^{w^{(d)}}$ accounts for residual variability, that is subject specific perturbations missed by the ODE model (1.1) but vanishing when averaging over the global population.
We are left with the parametrization of $\dot{h}_i^{(d)}(0)$, since the only constraint is to ensure its positiveness, we choose the log-transformation $\log \dot{h}_i^{(d)}(0) = \overline{\beta_0^{(d)}} + \sum_{m=1}^{M} \pi_{G_m}^{\beta_0^{(d)}} \mathbb{I}_{i \in G_m} + b_i^{\beta_0^{(d)}}$ with $b_i^{\beta_0^{(d)}} \sim N(0, \sigma_{\beta_0^{(d)}}^2)$.

From the proposed parametrization, we can turn the problem of assessing the presence of model misspecification and/or group effects ignored by the original model as statistical tests, in a way similar to the work of Hooker & Ellner [28]. There is presence of model misspecification (resp. treatment effect for the $m$-th group) for the $d$-th state variable evolution if $(\overline{\beta_1^{(d)}}, \overline{\delta\beta_k^{(d)}})$ (resp. $(\pi_{G_m}^{\beta_1^{(d)}}, \pi_{G_m}^{\delta\beta_k^{(d)}})$) is significantly different from 0. We can also easily reconstruct $h_{G_m}^{(d)}$, the time-warping function for the $m$-th group and $d$-th state variable by solving the initial value problem:

$$ \left\{ \begin{aligned} \dot{g}_{G_m}^{(d)} &= (\overline{w^{(d)}} + \pi_{G_m}^{w^{(d)}}) g_{G_m}^{(d)} \\ \dot{h}_{G_m}^{(d)} &= g_{G_m}^{(d)} \\ (g_{G_m}^{(d)}(0), h_{G_m}^{(d)}(0)) &= (e^{\overline{\beta_0^{(d)}} + \pi_{G_m}^{\beta_0^{(d)}}}, 0). \end{aligned} \right. $$

A similar procedure holds for finding $\overline{h}^{(d)}$, the mean time-warping function at the population scale for the $d$-th state variable.

## 2.3. Statistical inference of the unified model

Now that we have a parametric form for $(w_i^{(d)}, \dot{h}_i^{(d)}(0))$, the estimation of the time-distorted NLME-ODE (2.2) can be done with the addition of population time-warping parameters:

$$ \left\{ \begin{aligned} \overline{\beta} &= (\overline{\beta_0^{(d)}}, \overline{\beta_1^{(d)}}, (\overline{\delta\beta_k^{(d)}})_{k=2,...,K-1})_{d=1,...,D}, \\ \pi^{\beta} &= \left( \left( \pi_{G_m}^{\beta_0^{(d)}}, \pi_{G_m}^{\beta_1^{(d)}}, (\pi_{G_m}^{\delta\beta_k^{(d)}})_{k=2,...,K-1} \right)_{m=1,...,M} \right)_{d=1,...,D}, \\ \Sigma_{b^{\beta}} &= \mathrm{diag} \left( \left( \sigma_{\beta_0^{(d)}}^2, \sigma_{\beta_1^{(d)}}^2, (\sigma_{\delta\beta_k^{(d)}}^2)_{k=2,...,K-1} \right)_{d=1,...,D} \right), \end{aligned} \right. $$

to the structural ones $(\overline{\xi}, \pi^{\xi}, \Sigma_{b^{\xi}})$. From the observation model (1.2), we can then derive the sampling density for the observations which makes likelihood-
based inference a natural choice for time-warping function estimation. Moreover, the related estimator variance automatically provides uncertainty quantification for parameter inference or prediction purposes.

Still, even though we aim to estimate as precisely as possible the non-null functions $w_i^{(d)}$, indicating the presence of perturbations which cannot be accounted by original ODE parameter variations alone, we want to do so while avoiding over-fitting and the subsequent practical identifiability issues coming from it which can appear with the additional estimation of functions $\{h_i\}_{i=1,...,n}$ in our sparse data case. For this, we add a priori distributions on the parameters defining $(\bar{h}, \pi^h, b_i^h)$ in such a way that $h_i(t) \simeq t$. That is, we enforce $X_i \circ h_i$ to stay relatively close to the evolution speed induced by the original ODE (1.1) unless the data strongly advocates otherwise. It means that in case of disagreement between data and model prediction, we prioritize endogenous explanations in terms of structural parameter variations before moving to the external ones indicating the presence of model misspecification. This is in accordance with the interpretability and the confidence we have in the postulated ODE (1.1) which is assumed to be a mostly correct, yet possibly incomplete, summary of previous biological knowledge. For this, we add the prior distributions $\bar{\beta} \sim N(0, \alpha_{\beta}I_{DK})$ and $\pi^{\beta} \sim N(0, \alpha_{\beta}I_{DMK})$ which indeed favor parameter values leading to $w_i^{(d)}(t) \simeq 0$ and $\dot{h}_i^{(d)}(0) = 1$, the level of accepted departure from these reference values being quantified by the hyperparameter $\alpha_{\beta}$, in practice chosen according to data abundance and quality. This is equivalent to starting with the null hypothesis that there is no time-warping. By doing so, temporal distortion is considered when no structural ODE parameter covariates can easily accommodate group specific features. Moreover, to maintain subject specific temporal variations at a reasonable level, we set $\Sigma_{b^\beta}$ to low values chosen similarly as $\alpha_{\beta}$. To account for the specified prior distribution, we choose as estimator of the population parameters $(\bar{\xi}, \pi^\xi, \Sigma_{b^\xi}, \bar{\beta}, \pi^\beta)$ the maximum a posteriori and for the subject specific ones, the empirical Bayes estimates; the variance is estimated via the Fisher Information matrix. All of the numerical methods are implemented using the SAEM algorithm in Monolix [7].

*Remark 2.* More complex parametrizations could have been chosen for $w_i$, ranging from splines to Gaussian processes and Neural Networks. In practice, data sparsity drives us to choose this simple linear approximation to avoid identifiability problems for $(\bar{\beta}, \pi^\beta)$ inference. In particular, the sensitivity of $w_i^{(d)}$ with respect to $\beta_{ik}^{(d)}$ on interval $[t_k, t_{k+1}]$ makes it easier to aggregate
information coming from subjects with different observation timepoints for $(\bar{\beta}, \pi^{\beta})$ estimation.

## 3. numerical experiment

### 3.1. *Experimental setting*

We first test our methodology on simulated data to ascertain whether our procedure can accurately estimate both the parameters specifying the time-warping functions and the structural parameters within a time-warped NLME-ODE, with or without the presence of model error. As a test model, we consider the ODE proposed by [21] describing the evolution of target/infected CD4+T cells and HIV viral load:

$$ 
\begin{aligned}
\dot{\tilde{T}} &= \lambda - \phi \tilde{T} \tilde{V} - d_T \tilde{T} \\
\dot{\tilde{I}} &= \frac{\phi \tilde{T} V}{1+\gamma} - d_I \tilde{I} \\
\dot{V} &= k \tilde{I} - cV.
\end{aligned}
 $$

Here, it is slightly modified to account for the immune response. In this ODE model, the target cells $\tilde{T}$ are produced with a constant rate $\lambda$ and a proportion of them, quantified by the infection rate $\phi$, eventually become infected. Sill, only a fraction of these initially infected cells, quantified by the parameter $\gamma$, will become able to release new virions. It mirrors the cytotoxic immune response action on the latently infected cells before they enter into their active phase. These actively infected cells, denoted $\tilde{I}$, then produce viruses of concentration $V$ with rate $k$. The remaining terms $d_T, d_I$ and $c$ are respectively the death rate of $\tilde{T}$, $\tilde{I}$, and the viral clearance rate.

The goal here is to represent viral load rebound just after the interruption of a long antiretroviral therapy (ART) exposure. This design is also referred as analytical treatment interruption and used for therapeutic vaccine trials, for example in HIV [29]. So, at $t = 0$ we assume there is no infection $(\phi \simeq 0)$, and thus the system is near its steady state i.e $(\tilde{T}(0), \tilde{I}(0), V(0)) = (\lambda/d_T, 0, V_0)$ but still with $V_0 \neq 0$ to mimic the existence of latent viral reservoir [30]. Now, to stay close to the real observation framework considered in the next section, we only use measurement of viral load $V$ for parameter estimation, making the ODE not structurally-identifiable. Thus, to reduce the number of parameters, we focus instead on the evolution of $(T, I, V) =$
$$ (\tilde{\frac{\dot{I}}{\lambda}}, k\tilde{I}, V) \text{ ruled by the ODE:} $$

$$ 
\begin{aligned}
\dot{T} &= 1 - \phi TV - d_T T \\
\dot{I} &= \phi \Lambda TV - d_I I \\
\dot{V} &= I - cV
\end{aligned}
\qquad (3.1) $$

whose definition requires the structural parameters $(\phi, \Lambda, d_T, d_I, c)$ with $\Lambda := \frac{k\lambda}{\gamma+1}$ and initial conditions $(T(0), I(0), V(0)) = (1/d_T, 0, V_0)$. We set the mean parameter to the values given in [31], $\phi = 7.7 \times 10^{-7}$ viral load$^{-1}$.day$^{-1}$, $d_T = 0.05$ day$^{-1}$, $d_I = 0.4$ day$^{-1}$, $c = 23$ day$^{-1}$ as well as $\Lambda = 2.5 \times 10^6$ viral load.day$^{-2}$ (corresponding to $\lambda = 71$, $k = 5 \times 10^4$ and $\gamma = 0.4$) and $V(0) = 10^{-4}$.

From ODE (3.1), we then simulate data mimicking the temporal evolution of a population divided in subgroups receiving different treatments. For this, we separate the subjects into three groups, respectively denoted by $G_0, G_1$ and $G_2, G_0$ being the placebo group and $G_1, G_2$ the ones receiving actual treatments. To account for therapeutic actions on the immune response mechanism acknowledged by the ODE (3.1), we still assume that $\Lambda$ is centered around $\bar{\Lambda} = 2.5 \times 10^6$ but with group and subject specific variation given by:

$$ \ln \Lambda_i = \ln \bar{\Lambda} + \pi_{G_1}^{\ln \Lambda} \mathbb{I}_{i \in G_1} + \pi_{G_2}^{\ln \Lambda} \mathbb{I}_{i \in G_2} + b_i^{\ln \Lambda} \text{ with } b_i^{\ln \Lambda} \sim N(0, \sigma_{\ln \Lambda}^2 = 0.15^2). $$

We add temporal distortion to the ODE (3.1) with time-warping function following the finite basis decomposition given by (2.3) with $K = 2$ basis elements and elapsed time of $t_2 - t_1 = t_1 - t_0 = 15$ days. For all subjects, we simulate model error presence at the mean population level for $(T, I)$ dynamic, if we consider the ODE (3.1) as the correct model, by setting $(\overline{\beta_1^T}, \overline{\delta\beta_2^T}) = (0.25, 0)$, $(\overline{\beta_1^I}, \overline{\delta\beta_2^I}) = (0.1, 0.05)$ and $\dot{h}_i^{(T)}(0) = \dot{h}_i^{(I)}(0) = 1$. We also consider treatment specific time-warping variations on $\overline{\beta_1^I}$ and random effects on $(\beta_{i1}^T, \beta_{i1}^I)$, leading to the subsequent parametrization, as particular case of (2.4):

$$ 
\left\{
\begin{array}{l}
\beta_{i1}^T = \overline{\beta_1^T} + b_i^{\beta_1^T} \text{ with } b_i^{\beta_1^T} \sim N(0, \sigma_{\beta_1^T}^2 = 0.02^2) \\
\beta_{i2}^T = \beta_{i1}^T \\
\beta_{i1}^I = \overline{\beta_1^I} + \pi_{G_1}^{\beta_1^I} \mathbb{I}_{i \in G_1} + \pi_{G_2}^{\beta_1^I} \mathbb{I}_{i \in G_2} + b_i^{\beta_1^I} \text{ with } b_i^{\beta_1^I} \sim N(0, \sigma_{\beta_1^I}^2 = 0.02^2) \\
\beta_{i2}^I = \beta_{i1}^I + \overline{\delta\beta_2^I}.
\end{array}
\right.
$$
We then choose different values for $(\pi_{G_m}^{\ln \Lambda}, \pi_{G_m}^{\beta_1^I})_{m=1,2}$ to reproduce different treatment actions between the groups. For $m=1$, both $\Lambda$ and the time warping functions are changed by choosing $(\pi_{G_1}^{\ln \Lambda}, \pi_{G_1}^{\beta_1^I}) = (-0.5, 0.05)$. For $m=2$, only the time-warping parameter is affected with $(\pi_{G_2}^{\ln \Lambda}, \pi_{G_2}^{\beta_2^I}) = (0, -0.1)$.

We simulate a population of $n=150$ subjects equally partitioned between the groups $m=0, 1, 2$ for which we generate $n_i = 21$ viral load measurements each on the same interval of $[0, 40]$ weeks given by the observation model:

$$Y_{i,j} = \log_{10}(V(t_{i,j})) + \varepsilon_{i,j} \text{ with } \varepsilon_{i,j} \sim N(0, \sigma_{\varepsilon}^2 = 0.15^2),$$

with two measurements a week for the first four weeks followed by one measurement each week afterward. From this synthetic dataset, we estimate the mean parameter values $(\ln \bar{\Lambda}, \bar{\beta}_1^T, \bar{\beta}_1^I, \bar{\delta\beta}_2^I)$, the covariates $(\pi_{G_m}^{\ln \Lambda}, \pi_{G_m}^{\beta_1^I})_{m=1,2}$ and the variance $(\sigma_{\ln \Lambda}^2, \sigma_{\varepsilon}^2)$, the other ones being set to their true values. We also set $\alpha_{\beta} = 0.5$ in prior distributions $(\bar{\beta}_1^T, \bar{\beta}_1^I, \bar{\delta\beta}_2^I) \sim N(0, \alpha_{\beta}I_3)$ and $(\pi_{G_1}^{\beta_1^I}, \pi_{G_2}^{\beta_2^I}) \sim N(0, \alpha_{\beta}I_2)$.

We proceed to $N_{MC} = 10^3$ replicates of such simulated population data along with subsequent parameter estimation. From these replicates, we estimate the mean as well as the relative (i.e. scaled) values of the bias (Rel. Bias), the empirical standard deviation (Emp. Std) and the Mean Square Error (MSE) of our estimator to quantify its practical accuracy. We also compute the coverage rate of the 95% confidence intervals corresponding to the empirical standard deviation (Cov. R.). For a well conditioned estimation problem, the actual coverage rate of the theoretical 95% confidence interval should be close to 0.95.

Our goal here is to quantify our ability to retrieve the parameters associated with the time-warping functions as well as the effect of ignoring, in part or in totality, the temporal distortion on the estimation of ODE parameter covariates. That is why we consider three estimation cases depending on the assumed time-warping model and corresponding to an increasing level of model error presence:

**Case 1** Well-specified case: both assumed models for estimation and data simulation perfectly coincide.
**Case 2** Attenuated time-warping: parameter estimation is made from a model underestimating time-distortion action by stating that $\overline{\beta_1^T} = 0.1$ instead of $\overline{\beta_1^T} = 0.25$.

**Case 3** Missing time-warping: parameter estimation is made from the original ODE (3.1) which ignores time-distortion.

### 3.2. Results

The estimation results are given in Table 1. Regarding the well-specified case 1, the low MSE for all parameters as well as the consistency between the theoretical and actual coverage rate of the confidence intervals indicate a well conditioned estimation problem. In particular, even the time warping parameters ($\overline{\beta_1^T}, \overline{\beta_1^I}$) linked to the dynamic of unobserved state variables $T$ and $I$ are accurately estimated. We only observe a small coverage rate drop for $\overline{\delta\beta_2^I}$ estimation, translating the difficulty to infer this parameter linked to the dynamics of the system once this latter is about to enter into its steady state phase. Now in case 2, where the temporal distortion is wrongly attenuated, a degradation of estimation accuracy appears for all parameters which is particularly marked for the time-warping ones with a significant drop in coverage rates. Still, by comparing case 2 with case 3, we can see that accounting, even in a minimal way, for temporal distortions has a positive effect on inference of structural parameters which is even more degraded in this last setting. In particular, all coverage rates drop to 0 for $\overline{\Lambda}$ and its covariate. Moreover, this is accompanied with significantly negative estimated value for $\pi_{Gp_2}^{\ln\Lambda}$, leading to the false conclusion that both treatments have an action on the mechanism represented by $\Lambda$. This is due to compensating effects, structural parameter values try to account for external features acting on $T$ dynamic missed by the assumed model.

## 4. Real data analysis: siv viral load rebound after treatment cessation

In this section, we present a case study based on pre-clinical data collecting the SIV viral load evolution on non human primates (NHPs) after anti-retroviral therapy interruption [31]. We demonstrate the methodology presented in the previous sections and how it can impact data analysis in practice.
<table><caption>Table 1: Estimation results for Case 1: well defined model; Case 2: attenuated warping and Case 3: no time-warping.</caption><thead><tr><th></th><th>Case</th><th>Mean</th><th>Rel. Bias</th><th>Emp. Std</th><th>MSE</th><th>Cov. R.</th></tr></thead><tbody><tr><td rowspan="3">ln Λ̄ ≃ 14.7</td><td>1</td><td>14.7</td><td>-1.5 × 10<sup>-3</sup></td><td>2.2 × 10<sup>-2</sup></td><td>4.0 × 10<sup>-4</sup></td><td>0.96</td></tr><tr><td>2</td><td>14.8</td><td>1.1 × 10<sup>-2</sup></td><td>2.3 × 10<sup>-2</sup></td><td>6.5 × 10<sup>-4</sup></td><td>0.93</td></tr><tr><td>3</td><td>15.0</td><td>2.0 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-2</sup></td><td>8.0 × 10<sup>-4</sup></td><td>0.00</td></tr><tr><td rowspan="3">π<sub>Gp<sub>1</sub></sub> = -0.5</td><td>1</td><td>-5.1 × 10<sup>-1</sup></td><td>-2.0 × 10<sup>-2</sup></td><td>3.1 × 10<sup>-2</sup></td><td>1.3 × 10<sup>-3</sup></td><td>0.96</td></tr><tr><td>2</td><td>-4.8 × 10<sup>-1</sup></td><td>4.0 × 10<sup>-2</sup></td><td>3.2 × 10<sup>-2</sup></td><td>2.6 × 10<sup>-3</sup></td><td>0.94</td></tr><tr><td>3</td><td>-2.4 × 10<sup>-1</sup></td><td>5.2 × 10<sup>-1</sup></td><td>3.0 × 10<sup>-2</sup></td><td>2.6 × 10<sup>-1</sup></td><td>0.00</td></tr><tr><td rowspan="3">π<sub>Gp<sub>2</sub></sub> = 0</td><td>1</td><td>≈ 0.0</td><td>1.8 × 10<sup>-3</sup></td><td>3.2 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-3</sup></td><td>0.95</td></tr><tr><td>2</td><td>-0.2 × 10<sup>-1</sup></td><td>-0.2 × 10<sup>-1</sup></td><td>3.2 × 10<sup>-2</sup></td><td>1.4 × 10<sup>-3</sup></td><td>0.95</td></tr><tr><td>3</td><td>-2.8 × 10<sup>-1</sup></td><td>-2.8 × 10<sup>-1</sup></td><td>2.1 × 10<sup>-2</sup></td><td>7.8 × 10<sup>-2</sup></td><td>0.00</td></tr><tr><td rowspan="3">β̄<sub>1</sub><sup>T</sup> = 0.25</td><td>1</td><td>2.5 × 10<sup>-1</sup></td><td>1.8 × 10<sup>-2</sup></td><td>5.3 × 10<sup>-2</sup></td><td>3.2 × 10<sup>-3</sup></td><td>0.95</td></tr><tr><td>2</td><td>0.1 (Fixed)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>3</td><td>0 (Fixed)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">β̄<sub>1</sub><sup>I</sup> = 0.1</td><td>1</td><td>1.0 × 10<sup>-1</sup></td><td>8.1 × 10<sup>-3</sup></td><td>3.4 × 10<sup>-2</sup></td><td>1.2 × 10<sup>-3</sup></td><td>0.95</td></tr><tr><td>2</td><td>9.2 × 10<sup>-2</sup></td><td>-4.6 × 10<sup>-2</sup></td><td>3.4 × 10<sup>-2</sup></td><td>4.4 × 10<sup>-2</sup></td><td>0.74</td></tr><tr><td>3</td><td>0 (Fixed)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">π<sub>Gp<sub>1</sub></sub> = 0.05</td><td>1</td><td>5.1 × 10<sup>-2</sup></td><td>2.1 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-1</sup></td><td>1.2 × 10<sup>-2</sup></td><td>0.96</td></tr><tr><td>2</td><td>5.5 × 10<sup>-2</sup></td><td>7.3 × 10<sup>-2</sup></td><td>1.1 × 10<sup>-1</sup></td><td>1.6 × 10<sup>-2</sup></td><td>0.90</td></tr><tr><td>3</td><td>0 (Fixed)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">π<sub>Gp<sub>2</sub></sub> = -0.1</td><td>1</td><td>-1.0 × 10<sup>-1</sup></td><td>1.6 × 10<sup>-3</sup></td><td>5.0 × 10<sup>-2</sup></td><td>2.5 × 10<sup>-3</sup></td><td>0.95</td></tr><tr><td>2</td><td>-9.5 × 10<sup>-2</sup></td><td>-4.1 × 10<sup>-3</sup></td><td>5.1 × 10<sup>-2</sup></td><td>4.3 × 10<sup>-3</sup></td><td>0.87</td></tr><tr><td>3</td><td>0 (Fixed)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">δβ̄<sub>2</sub><sup>I</sup> = 0.05</td><td>1</td><td>4.1 × 10<sup>-2</sup></td><td>-1.7 × 10<sup>-1</sup></td><td>2.1 × 10<sup>-1</sup></td><td>8.0 × 10<sup>-2</sup></td><td>0.91</td></tr><tr><td>2</td><td>4.5 × 10<sup>-2</sup></td><td>-1.0 × 10<sup>-1</sup></td><td>2.2 × 10<sup>-1</sup></td><td>5.5 × 10<sup>-2</sup></td><td>0.94</td></tr><tr><td>3</td><td>0 (Fixed)</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>
### 4.1. Available data

We analyze the measured viral load concentration evolution obtained from studies presented in [22, 23] focusing on SIV positive NHPs. In these pre-clinical trials, NHPs received antiretroviral therapy (ART) during several months after SIV infection. ART is then stopped and viral load is measured at different timepoints afterward. The point of such studies is to quantify the effect of treatments the subjects have received in addition to ART on the viral rebound dynamic. The elapsed time between infection and ART initiation is study dependent, nine weeks for the one presented in Lim et al. (comprising 19 animals and now denoted study 1) and one week in Borducchi et al. (comprising 34 animals and named study 2 afterward). In study 1, the NHPs either received ART alone (control group) or ART+Tlr7 (a molecule involved in the antiviral innate immune system response), on average $n_i = 20$ measurement points are available per animal with final observation made $t = 112$ days after ART cessation. In study 2, the effect of a therapeutic vaccine has also been investigated for which we have access to $n_i = 13$ measures with a follow-up ending at $t = 168$ days after treatment interruption. Once gathered, the data are clustered into four treatment groups: NHPs treated with ART alone, ART+Tlr7, ART+Vaccine or ART+Tlr7+Vaccine. The longitudinal follow-up is presented in Figure 4.1.

### 4.2. Model & parametrization

As in the previous section, we consider the time-warped version of the ODE-based statistical model (3.1) with $d_T = 0.05$, $d_I = 0.4$ and $c = 23$ (in day$^{-1}$) and we are left with $(\phi, V(0), \Lambda)$ as the unknown ODE parameters. To quantify treatment action on newly infected cells and their subsequent viral production, we consider group-specific values for $\Lambda$. In order to account for treatment effects on mechanisms which are not described by the ODE model (3.1) and possibly acting on all state-variables $(T, I, V)$, we add time warping functions. We assume that $\dot{h}_i(0) = 1$ for all subjects, and we choose $K = 2$ elements elapsed with $t_{k+1} - t_k = 5$ days for the parametrization of $(w_i^T, w_i^I, w_i^V)$. We add random effects on $\Lambda_i$ and $(\overline{\beta_1^{(T)}}, \overline{\beta_1^{(I)}}, \overline{\beta_1^{(V)}})$ to end up with the following parametrization, specifying both structural ODE param-
<figure><img src="image_2.png" /><figcaption>Figure 4.1: Raw data available in [22, 23]. Top: study 1 with the *n* = 8 subjects of control group (left) and the *n* = 11 ones of Tlr7 group (right). Middle: study 2 with the *n* = 9 subjects of control group (left) and the *n* = 8 ones of Tlr7 group (right). Bottom: study 2 with the *n* = 9 subjects of vaccine group (left) and the *n* = 8 ones of Tlr7+vaccine group (right).</figcaption></figure>
eters and the time-warping ones as in equation (2.4):

$$
\left\{
\begin{array}{l}
\ln \Lambda_i = \ln \bar{\Lambda} + \sum_m \pi_{G_m}^{\ln \Lambda} 1_{i \in G_m} + b_i^{\ln \Lambda} \text{ with } b_i^{\ln \Lambda} \sim N(0, \sigma_{\ln \Lambda}^2), \\
\beta_{i1}^{(d)} = \bar{\beta}_1^{(d)} + \sum_m \pi_{G_m}^{\beta_1^{(d)}} \mathbb{I}_{i \in G_m} + b_i^{\beta_1^{(d)}} \text{ with } b_i^{\beta_1^{(d)}} \sim N(0, \sigma_{\beta_1^{(d)}}^2), \\
\beta_{i2}^{(d)} = \beta_{i1}^{(d)} + \delta \beta_2^{(d)}
\end{array}
\right.
$$

for (d) ∈ {T, I, V} and G_m ∈ {Tlr7, Vac, Tlr7&Vac} for NHPs treated with Tlr7, therapeutic vaccine and both, the NHPs receiving ART alone being chosen as the control group (i.e. G_0).

The limited number of available subjects per arm drives us to drastically limit the allowed inter-subject variability if we want to infer a common population trend, that is why we choose the small value $\sigma_{\beta_1^{(d)}}^2 = 0.02^2$ for the random-effect variance of time-warping. This data sparsity also leads us to choose the weak value of $\alpha_{\beta} = 0.2$ for the specification of the prior distribution of $\bar{\beta}$ and $\pi^{\beta}$ to isolate group/state-variable specific temporal distortion while avoiding over fitting. For the same reason, we have considered $(\phi, V(0))$ as nuisance parameters to remove practical identifiability issues and we have estimated them via profiled likelihood on the grid $[10^{-5}, 10^{-4}] \times [10^{-4}, 10^{-2}]$. We are left with the joint estimation of:

$$
\left\{
\begin{array}{l}
\left( \ln \bar{\Lambda}, \pi_{\text{Tlr7}}^{\ln \Lambda}, \pi_{\text{Vac}}^{\ln \Lambda}, \pi_{\text{Tlr7\&Vac}}^{\ln \Lambda}, \sigma_{\ln \Lambda}^2 \right), \\
\left( \bar{\beta}_1^{(d)}, \delta \beta_2^{(d)}, \pi_{\text{Tlr7}}^{\beta_1^{(d)}}, \pi_{\text{Vac}}^{\beta_1^{(d)}}, \pi_{\text{Tlr7\&Vac}}^{\beta_1^{(d)}} \right)
\end{array}
\right.
\quad (d) \in \{T,I,V\}
$$

485 done by a posteriori maximization.

4.3. *Inference strategy for treatment action*

Now, to clearly separate the specific effect of each treatment on the mech-
anism represented by Λ from other actions on the dynamics of T, I or V while
acknowledging the impact of model misspecification on estimation accuracy,
we proceed to the following sequential estimation strategy in three steps
based on models of increasing complexity.

492 **Step 1, original model inference** We estimate the ODE (3.1) without
time-warping.
493

494 **Step 2, mean time-warped model inference** We estimate distorted ver-
sions of ODE (3.1) with state-specific time-warping presence tested on
495
any possible combinations within {$T, I, V$}, but without treatment specific covariates on them i.e. with
$$ \left( \pi_{\text{Tlr7}}^{\beta_1^{(d)}}, \pi_{\text{Vac}}^{\beta_1^{(d)}}, \pi_{\text{Tlr7&Vac}}^{\beta_1^{(d)}} \right)_{(d) \in \{T, I, V\}} = 0. $$

**Step 3, treatment-specific time-warped model inference** Addition and estimation of $\left( \pi_{\text{Tlr7}}^{\beta_1^{(d)}}, \pi_{\text{Vac}}^{\beta_1^{(d)}}, \pi_{\text{Tlr7&Vac}}^{\beta_1^{(d)}} \right)$ on the subset of state variables within {$T, I, V$} for which the inclusion of time-warping has led to the most adequate model in step 2.

The ability to compare the adequacy of the different tested models requires a common criterion, we choose the Akaike Information Criterion (AIC) for this task which balances model fidelity to the data and its parsimony; quantified by its number of parameters. For all investigated models, their estimation are done in presence and absence of treatment covariate on $\Lambda$ i.e. when $1/(\pi_{\text{Tlr7}}^{\ln\Lambda}, \pi_{\text{Vac}}^{\ln\Lambda}, \pi_{\text{Tlr7&Vac}}^{\ln\Lambda})$ are estimated or 2/are set to 0. This is here to ensure that a potential treatment effect found on $\Lambda$ is robust to model error presence and is not likely to be explained from compensating effects alone. Step 1 is used as a benchmark to quantify the differences in terms of model adequacy and the conclusions we hold regarding treatment effect on $\Lambda$ between the original ODE and its different time-distorted versions. From there, the goal of step 2 is to precisely determine which state variable specific time-distortions need to be accounted for and which state equations are the most prone to model error and/or external treatment action. Still, step 2 cannot tells us if the required time-distortion is symptomatic of mean model error presence affecting all subjects, no matter the received treatments, or if it indicates treatment specific action. That is the point of step 3 which investigates the need for time-warping covariates.

## 4.4. Results

The model adequacy results given in AIC for models of steps 1, 2 and 3, both in absence and presence of covariates on $\Lambda$, are presented in Table 2. The parameter estimation corresponding to the best retained model of step 1 (Original ODE) and step 3 (Time-warped ODE with treatment specific covariates) are presented in Table 3.

The original ODE estimation (step 1) detects a significant effect of treatment on the immune system mechanism represented by $\Lambda$. Moreover, accounting also for external disturbances (step 2) always improves model adequacy no matter on which state variables time-warping is added. The lower
<table><caption>Table 2: Model adequation results for step 1, step 2 and step 3. The results are given in AIC with the lowest one obtained at each step indicated in bold.</caption><thead><tr><th>Model</th><th>No covariate on Λ</th><th>Covariate on Λ</th></tr></thead><tbody><tr><td>Step 1</td><td></td><td></td></tr><tr><td></td><td>1576</td><td><strong>1547</strong></td></tr><tr><td>Step 2, time warping on:</td><td></td><td></td></tr><tr><td><i>T</i></td><td>1532</td><td>1503</td></tr><tr><td><i>I</i></td><td>1449</td><td>1433</td></tr><tr><td><i>V</i></td><td>1386</td><td>1350</td></tr><tr><td>(<i>T</i>, <i>I</i>)</td><td>1002</td><td><strong>966</strong></td></tr><tr><td>(<i>T</i>, <i>V</i>)</td><td>1024</td><td>995</td></tr><tr><td>(<i>I</i>, <i>V</i>)</td><td>1370</td><td>1346</td></tr><tr><td>(<i>T</i>, <i>I</i>, <i>V</i>)</td><td>1000</td><td>967</td></tr><tr><td>Step 3, time warping covariates on:</td><td></td><td></td></tr><tr><td><i>T</i></td><td>990</td><td>950</td></tr><tr><td><i>I</i></td><td>1006</td><td>959</td></tr><tr><td>(<i>T</i>, <i>I</i>)</td><td>984</td><td><strong>942</strong></td></tr></tbody></table>

AIC is obtained by adding time-warping on *T* and *I* and thus is then used as reference model for step 3. Interestingly, when considering temporal distortion on only one state variable at a time, *V* gives the lowest AIC. Still, once model error is corrected for target/infected cell dynamics, the equation ruling *V* is considered as statistically correct at the population scale. Additionally, no matter the tested structure for the time-distortion, considering covariates on Λ always improves AIC. From step 3, we can see that, in addition to model error, treatment specific features are significantly missed by the original ODE on *T* and *I* dynamic. Thus, even though we proceed to treatment inference from an incomplete model, this procedure allows us to separate what is due to general model misspecification from actual treatment effect in the case of disagreement between the assumed model dynamics and the raw data kinetics.

Now, let us focus on parameter estimation results. No matter the considered models, we derive a treatment effect on Λ for vaccine alone or in combination with Tlr7, but not for Tlr7 alone. The closeness between the estimation results of ($\pi_{Tlr7}^{\ln\Lambda}$, $\pi_{Vac}^{\ln\Lambda}$, $\pi_{Tlr7&Vac}^{\ln\Lambda}$) from the original ODE and its time-warped version supports the hypothesis that there is indeed a vaccine effect, alone or in combination with Tlr7, on the viral production from newly
<table><caption>Table 3: Estimation results for original ODE (Step 1) and ODE with mean time-distortion and treatment effects on (T, I) (Step 3). The values between brackets represent the 95% confidence interval. Bold values are significantly different from 0, either indicating significative treatment effect or model error presence.</caption><thead><tr><th></th><th>Original Model (Step 1)</th><th>Treatment specific time-warped model (Step 3)</th></tr></thead><tbody><tr><td>φ</td><td>1.2 × 10<sup>-5</sup> (Profiled)</td><td>2.8 × 10<sup>-5</sup> (Profiled)</td></tr><tr><td>V(0)</td><td>3.3 × 10<sup>-4</sup> (Profiled)</td><td>9.1 × 10<sup>-4</sup> (Profiled)</td></tr><tr><td>ln Λ</td><td>12.3 [12.2; 12.5]</td><td>11.8 [11.7; 12.0]</td></tr><tr><td>π<sub>Tlr7</sub><sup>ln Λ</sup></td><td>0.0 [-0.1; 0.1]</td><td>0.0 [-0.1; 0.1]</td></tr><tr><td>π<sub>Vac</sub><sup>ln Λ</sup></td><td>-0.2 [-0.3; -0.1]</td><td>-0.2 [-0.3; -0.1]</td></tr><tr><td>π<sub>Tlr7&Vac</sub><sup>ln Λ</sup></td><td>-0.4 [-0.5; -0.3]</td><td>-0.4 [-0.6; -0.2]</td></tr><tr><td>β<sub>1</sub><sup>T</sup></td><td>0 (Fixed)</td><td>-7.3 × 10<sup>-2</sup> [-1.0 × 10<sup>-1</sup>; -4.0 × 10<sup>-2</sup>]</td></tr><tr><td>π<sub>Tlr7</sub><sup>β<sub>1</sub><sup>T</sup></sup></td><td>0 (Fixed)</td><td>2.6 × 10<sup>-2</sup> [1.4 × 10<sup>-2</sup>; 3.8 × 10<sup>-2</sup>]</td></tr><tr><td>π<sub>Vac</sub><sup>β<sub>1</sub><sup>T</sup></sup></td><td>0 (Fixed)</td><td>4.5 × 10<sup>-2</sup> [2.8 × 10<sup>-2</sup>; 6.4 × 10<sup>-2</sup>]</td></tr><tr><td>π<sub>Tlr7&Vac</sub><sup>β<sub>1</sub><sup>T</sup></sup></td><td>0 (Fixed)</td><td>6.4 × 10<sup>-2</sup> [4.5 × 10<sup>-2</sup>; 8.2 × 10<sup>-2</sup>]</td></tr><tr><td>δβ<sub>2</sub><sup>T</sup></td><td>0 (Fixed)</td><td>-6.4 × 10<sup>-2</sup> [-6.5 × 10<sup>-2</sup>; -6.3 × 10<sup>-1</sup>]</td></tr><tr><td>β<sub>1</sub><sup>I</sup></td><td>0 (Fixed)</td><td>8.8 × 10<sup>-3</sup> [-1.5 × 10<sup>-2</sup>; 3.0 × 10<sup>-2</sup>]</td></tr><tr><td>π<sub>Tlr7</sub><sup>β<sub>1</sub><sup>I</sup></sup></td><td>0 (Fixed)</td><td>1.3 × 10<sup>-2</sup> [4.2 × 10<sup>-3</sup>; 2.3 × 10<sup>-2</sup>]</td></tr><tr><td>π<sub>Vac</sub><sup>β<sub>1</sub><sup>I</sup></sup></td><td>0 (Fixed)</td><td>2.4 × 10<sup>-3</sup> [9.1 × 10<sup>-3</sup>; 3.7 × 10<sup>-2</sup>]</td></tr><tr><td>π<sub>Tlr7&Vac</sub><sup>β<sub>1</sub><sup>I</sup></sup></td><td>0 (Fixed)</td><td>2.9 × 10<sup>-2</sup> [4.9 × 10<sup>-3</sup>; 5.3 × 10<sup>-2</sup>]</td></tr><tr><td>δβ<sub>2</sub><sup>I</sup></td><td>0 (Fixed)</td><td>-1.2 × 10<sup>-1</sup> [-1.3 × 10<sup>-1</sup>; -1.1 × 10<sup>-1</sup>]</td></tr><tr><td>-2LL</td><td>1535</td><td>910</td></tr><tr><td>AIC</td><td>1547</td><td>942</td></tr></tbody></table>
infected cells which is not a statistical artifact coming from the estimation of an over-simplified model. Regarding time-warping parameters now, the estimates of $(\overline{\beta_1^T}, \overline{\delta\beta_2^T}, \overline{\delta\beta_2^I})$ are significantly different from 0. This supports the hypothesis of model error presence at the population scale in the equation ruling the evolution of $(T, I)$ in ODE (3.1). Moreover, since we significantly got $(\pi_{Tlr7}^{\beta_1^T}, \pi_{Vac}^{\beta_1^T}, \pi_{Tlr7&Vac}^{\beta_1^T}) \neq 0$ and $(\pi_{Tlr7}^{\beta_1^I}, \pi_{Vac}^{\beta_1^I}, \pi_{Tlr7&Vac}^{\beta_1^I}) \neq 0$ as well, we infer additional effects from the time-warped ODE of all tested treatment strategies on the target/infected cell kinetics compared to the original model. In particular, we infer a therapeutic effect of Tlr7 alone which was originally missed.

From the parameter estimation, we can predict the evolution of $(T, I, V)$ from both the original ODE and its time-warped version, they are presented in figure 4.2. We can also reconstruct the mean time-warping functions $(\overline{h^T}, \overline{h^I})$, i.e. corresponding to the mean parameter values estimated via maximum *a posteriori*, as well as the treatment specific ones $(h_{Tlr7}^T, h_{Vac}^T, h_{Tlr7&Vac}^T)$ and $(h_{Tlr7}^I, h_{Vac}^I, h_{Tlr7&Vac}^I)$ when additionally accounting for covariate estimators this time, they are presented in figure 4.3.

The temporal distortion has a significant effect on the predicted viral load evolution as well as the underlying target/infected cells, even though both models predict that the long term viral load is best controlled with the combination vaccine+Tlr7, followed by vaccine alone, then by Tlr7 alone. The time-warped model predicts a significantly lower viral load and scaled infected cell population at steady state comparing to the ODE (3.1) for all groups. Moreover, the time-warped ODE presents a more lagged transient phase which is consistent with the estimated values of $(\overline{\beta_1^T}, \overline{\delta\beta_2^T}, \overline{\delta\beta_2^I})$ which are all significantly negative. It leads the mean time-warping functions $(\overline{h^T}, \overline{h^I})$ to be below the natural temporal evolution $h(t) = t$ on the restricted observation interval $[0, 20]$ days. This indicates that the ODE (3.1) alone drives $(T, I)$, and consequently $V$, to have a faster evolution, than the ones derived from data at the beginning of infection rebound. Thus, consistent model improvement should lead to a slower evolution of $(T, I)$ at the beginning of treatment interruption.
<figure><img src="image_3.png" /><figcaption>Figure 4.2: Predicted evolution up to 100 days after treatment cessation of scaled target cells *T* (top), scaled infected cells *I* (middle), viral load *V* (bottom) for the original model (left) and its time-warped version (right). The shaded area represents the 95% confidence interval.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4.3: Reconstituted time-warping functions up to 20 days after treatment cessation for T: ($\hat{h}^T$, $h_{Tlr7}^T$, $h_{Vac}^T$, $h_{Tlr7&Vac}^T$) (left) and I: ($\hat{h}^I$, $h_{Tlr7}^I$, $h_{Vac}^I$, $h_{Tlr7&Vac}^I$) (right). The grey dashed lines represents $h(t) = t$. The shaded area represents the 95% confidence interval.</figcaption></figure>
# 5. Conclusion

In this work, we promote a unifying framework incorporating curve registration techniques into the NLME-ODE setting classically considered for clinical/pre-clinical trial data analysis. By adding time-warping at the dynamic scale, we describe the evolution of a given subject via a perturbed version of an initial ODE, seen as a reference model summarizing the available knowledge we have on the biological process of interest. Our approach leads to ODE parameter estimates that are more robust to model misspecification as illustrated in simulated data cases. It also provides model error assessment techniques for each equation ruling the evolution of each state variable given as classic statistical test on parameter values. Additionally, it allows us to investigate the presence of group specific alterations on the unobserved state variables dynamics with descriptive approaches. By relying on a finite basis parametrization for the time-warping functions, we can consider their joint estimation with the ODE parameters as a classic parametric inference problem in an NLME-ODE, naturally solved with likelihood-based procedures.

From this, we can assess the effect of a therapeutic strategy on both 1/mechanisms described by the original ODE but also on 2/pathways ignored by this causal model but represented as temporal distortion using time-warping functions. The benefits of all these features are illustrated by the analysis of longitudinal viral load measurement made on SIV positive NHPs. Our procedure infers the enhancement brought by different therapeutic strategies on an immune mechanism modeled by the underlying ODE-based model while mitigating model error presence action on covariate estimation accuracy. Moreover, it allows us to move beyond the classic mechanistic based inference by quantifying treatment effects missed by the assumed ODE on the kinetics of unobserved cells populations. These results are consistent with the fact that we analyze immunotherapy treatments which primarily act on the immune response host and its ability to prevent infection and/or kills latently infected cells. So, such therapies are expected to have direct influences on target/infected cell dynamics but only indirect actions on viral load evolution. In particular, such treatments are supposed to control the size of the latent viral reservoir made up of long-term latently infected cells [30]. These cell populations being ignored by the ODE (3.1), we cannot properly account for treatment action on this mechanism via covariate specification on some of the structural parameters. This can explain the
618
619
620
621

under-estimation of treatment effect via the original ODE which is, at least
in part, detected and corrected by its time-warped version. Additionally, our
method enlightens a pattern in the apparent model error, providing guidance
for model refinement.

622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643

Still, most of the obtained results in the real data case rely on the reason-
able dimension of both the considered ODE (three states variables) and the
time-warping parametrization (two parameters per state variable for mean
model error, one per state variable and treated group for treatment effect).
Some gaps need to be filled to adapt and automatize the presented inference
procedure to more complex models. Refinement of time-warping functions
would then require multiple testing procedures tailored to their finite basis
decomposition to rigorously assess the presence of model error or group effect.
Alternatively, we can construct tests on $\bar{h}^{(d)}$ or $h_{G_m}^{(d)}$ to mitigate the impor-
tance of parametrization choice but this would demand adapting statistical
tests to infinite dimensional functional spaces. In the presented analysis,
all possible locations for time-warping among the three state variables have
been tested. However, a similar approach for a higher dimensional ODE
would have led to a combinatorial explosion, a problem aggravated by the
presence of group specific temporal distortions. Hence, there is a need for a
proper covariate selection algorithm adapted to our functional parameters in
more complex models. Finally, we may require additional structures for time-
 warping, in particular to construct causality preserving transformations. For
example, if we consider $(X^{(1)}, X^{(2)})$ with $X^{(2)}$ evolution known to be totally
driven by $X^{(1)}$, we may want to enforce that their distorted counterparts
$X^{(1)} \circ h^{(1)}$ and $X^{(2)} \circ h^{(2)}$ respect $h^{(1)}(t) \le h^{(2)}(t)$ i.e. that the evolution of
$X^{(1)}$ precedes the one of $X^{(2)}$.

644
645
646
647
648
649

By overcoming these issues, we can extend this procedure for model error
location and correction in higher dimensional settings in a way that can in-
corporate pre-existing knowledge about the causal relationships linking the
variables of interest. This will provide a valuable device for the construc-
tion of mechanistic models used for complex population longitudinal data
analysis. This is left for future development.

650

Acknowledgements

651
652
653

This work was supported by the Investissements d'Avenir program man-
aged by the ANR under reference ANR-10-LABX-77 Vaccine research insti-
tute.
654

655

656

657

658

659

660

661

John Fricks received partial financial support from the French government in the framework of the France 2030 program IdEx Université de Bordeaux.

The funders have no role in the study design, data collection, data analysis, data interpretation, or writing the report.

We thank Lixoft SAS for their support in using Monolix. Numerical computations were in part 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 (https://www.plafrim.fr).

662

## Data accessibility/Results Reproducibility

663

664

665

An R/Monolix implementation able to simulate the synthetic data and estimate the model parameters to reproduce the results presented in section 3 is accessible in the supplementary materials.

666

667

The real data used in section 4 are accessible via the articles [22, 23] where they are originally presented.

668

## References

669

670

671

[1] A. F. M. Lavielle, A. Samson, F. Mentre, Maximum likelihood estimation of long terms hiv dynamic models and antiviral response., Biometrics 67 (2011) 250-259.

672

673

674

[2] L. Villain, D. Commenges, C. Pasin, M. Prague, R. Thiébaut, Adaptive protocols based on predictions from a mechanistic model of the effect of il7 on cd4 counts, Statistics in medicine 38 (2) (2019) 221-235.

675

676

677

678

[3] M. Alexandre, R. Marlin, M. Prague, S. Coleon, N. Kahlaoui, S. Cardinaud, T. Naninck, B. Delache, M. Surenaud, M. Galhaut, et al., Modelling the response to vaccine in non-human primates to define sars-cov-2 mechanistic correlates of protection, elife 11 (2022) e75427.

679

680

681

682

683

[4] Q. Clairon, M. Prague, D. Planas, T. Bruel, L. Hocqueloux, T. Prazuck, O. Schwartz, R. Thiébaut, 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) (2023) e1011282.
684

685

686

687

[5] C. Pasin, I. Balelli, T. Van Effelterre, V. Bockstal, L. Solforosi, M. Prague, M. Douoguih, R. Thiébaut, Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation, Journal of virology 93 (18) (2019) e00579-19.

688

689

690

691

[6] M. Prague, D. Commenges, J. M. Gran, B. Ledergerber, J. Young, H. Furrer, R. Thiébaut, Dynamic models for estimating the effect of haart on cd4 in observational studies: Application to the aquitaine cohort and the swiss hiv cohort study, Biometrics 73 (1) (2017) 294-304.

692

693

[7] Lixoft, Monolix version 2023r1, Antony, France (2023).
URL http://lixoft.com/products/monolix/

694

695

696

[8] B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Bencourt, M. Brubaker, J. Guo, P. Li, A. Riddell, Stan: A probabilistic programming language, Journal of statistical software 76 (2017) 1-32.

697

698

699

700

701

[9] 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 (2020) 2020-11.

702

703

[10] J. Brynjarsdottir, A. O'Hagan, Learning about physical parameters: The importance of model discrepancy, Inverse Problems 30 (2014) 24.

704

705

[11] P. Kirk, D. Silk, M. Michael, Reverse engineering under uncertainty, in: Uncertainty in Biology, Springer, 2016, pp. 15-32.

706

707

708

[12] Y. Huang, G. Dagne, A bayesian approach to joint mixed-effects models with a skew normal distribution and measurement errors in covariates, Biometrics 67 (2011) 260-269.

709

710

[13] J. Ramsay, B. Silverman, Functional Data Analysis, Springer Series in Statistics, Springer, 2005.

711

712

[14] J. S. Morris, Functional regression, Annual Review of Statistics and Its Application 2 (1) (2015) 321-359.

713

714

[15] H. L. Shang, A survey of functional principal component analysis, AStA Advances in Statistical Analysis 98 (2) (2014) 121-142.
715

716

[16] J. Jacques, C. Preda, Functional data clustering: a survey, Advances in Data Analysis and Classification 8 (3) (2014) 231–255.

717

718

[17] J. O. Ramsay, X. Li, Curve registration, Journal of the Royal Statistical Society Series B: Statistical Methodology 60 (2) (1998) 351–363.

719

720

721

[18] J. S. Marron, J. O. Ramsay, L. M. Sangalli, A. Srivastava, Functional data analysis of amplitude and phase variation, Statistical Science (2015) 468–484.

722

723

724

725

[19] J. Matuk, K. Bharath, O. Chkrebtii, S. Kurtek, Bayesian framework for simultaneous registration and estimation of noisy, sparse, and fragmented functional data, Journal of the American Statistical Association 117 (540) (2022) 1964–1980.

726

727

728

[20] F. Yao, H.-G. Müller, J.-L. Wang, Functional data analysis for sparse longitudinal data, Journal of the American statistical association 100 (470) (2005) 577–590.

729

730

[21] A. S. Perelson, R. M. Ribeiro, Modeling the within-host dynamics of hiv infection, BMC biology 11 (2013) 1–10.

731

732

733

734

[22] E. N. Borducchi, C. Cabral, K. E. Stephenson, J. Liu, P. Abbink, D. Nganga, J. P. Nkolola, A. L. Brinkman, L. Peter, B. C. Lee, et al., Ad26/mva therapeutic vaccination with tlr7 stimulation in siv-infected rhesus monkeys, Nature 540 (7632) (2016) 284–287.

735

736

737

738

739

[23] S.-Y. Lim, C. E. Osuna, P. T. Hraber, J. Hesselgesser, J. M. Gerold, T. L. Barnes, S. Sanisetty, M. S. Seaman, M. G. Lewis, R. Geleziunas, et al., Tlr7 agonists induce transient viremia and reduce the viral reservoir in siv-infected rhesus macaques on antiretroviral therapy, Science translational medicine 10 (439) (2018) eaao4521.

740

741

742

[24] P. Chassat, Functional and shape data analysis under the frenet-serret framework: Application to sign language motion trajectories analysis, Ph.D. thesis, Université Paris-Saclay (2024).

743

744

745

[25] J. Park, N. Brunel, P. Chassat, Curvature and torsion estimation of 3d functional data: A geometric approach to build the mean shape under the frenet serret framework, arXiv preprint arXiv:2203.02398 (2022).
746

747

[26] A. Srivastava, E. P. Klassen, Functional and shape data analysis, Vol. 1, Springer, 2016.

748

749

750

[27] M. C. Kennedy, A. O'Hagan, Bayesian calibration of computer models, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63 (3) (2001) 425-464.

751

752

753

[28] G. Hooker, S. P. Ellner, et al., Goodness of fit in nonlinear dynamics: misspecified rates or misspecified states?, The Annals of Applied Statistics 9 (2) (2015) 754-776.

754

755

756

757

[29] M. Alexandre, M. Prague, R. Thiébaut, Between-group comparison of area under the curve in clinical trials with censored follow-up: Application to hiv therapeutic vaccines, Statistical Methods in Medical Research 30 (9) (2021) 2130-2147.

758

759

760

[30] J. D. Siliciano, R. F. Siliciano, The latent reservoir for hiv-1 in resting cd4+ t cells: a barrier to cure, Current Opinion in HIV and AIDS 1 (2) (2006) 121-128.

761

762

763

[31] M. Prague, J. M. Gerold, I. Balelli, C. Pasin, J. Z. Li, D. H. Barouch, J. B. Whitney, A. L. Hill, Viral rebound kinetics following single and combination immunotherapy for hiv/siv, BioRxiv (2019) 700401.