Olivier Brousse

Gilles Sassatelli

Thierry Gil

Michel Robert

François Grize
email: francois.grize@unil.ch

Eduardo Sanchez

Andrés Upegui

Yann Thoma

The Perplexus Programming Framework: Combining Bio-inspiration and Agent-Oriented Programming for the Simulation of Large Scale Complex Systems

Keywords: Bio-inspired, Scalability, Distributed platform, Ad-hoc Network, Multi-Agent Framework

This paper presents a flexible agent-oriented programming framework that provides native support for bio-inspired mechanisms. This solution, developed within the Perplexus IST European project [IST-034632] 1 aims at providing a means for the specification of applications running on a distributed and pervasive network of mobile nodes. In such applications, the deployed systems may face time-changing environments and bio-inspiration may prove useful bringing self-adaptability to the system. The presented framework features are demonstrated on a proof-of-concept application made of simple robots that autonomously improve their behaviour over time.

Introduction

The Perplexus Project aims at developing a scalable and ubiquitous platform endowed with bio-inspired features to simulate complex phenomena.

The developed platform is composed of a cluster of communicating nodes called UBIDULES (UBIquitous moDULES). These nodes have a microprocessor and a specific reconfigurable integrated circuit; both providing support to distributed bio-inspired mechanisms. Furthermore, application-specific hardware sensors and actuators may be added to any Ubidule for application-specific purposes (environmental interactions). The Perplexus Ubidules integrate at least one bio-inspired chip and a XScale PXA270 microprocessor that is in charge of networking, sensors, actuators and more generally software aspects. Figure 1 gives an overview of the obtained platform which is represented in form of a network of Ubidules. Each Ubidule is in turn represented as a unit made of a hardware partition and a software partition. An embedded Linux operating system runs on the microprocessors and handles most high-level functionalities.

Bio-inspired Mechanisms

Bio-inspiration has found applications in many areas ranging from finance to integrated circuits designs. Many theories exist regarding the emergence of life and evolution of biologic lifeforms; however it is often considered that three main mechanisms exist:

1. Phylogenesis deals with the evolution of a set of species. Evolution gears species towards a better adaptation of individuals to their environment; genetic algorithms are inspired from this very principle of life. 2. Ontogenesis describes the origin and the development of an organism from the fertilized egg to its mature form. Biological processes like healing and fault tolerance are qualified of ontogenetic. 3. Epigenesis refers to features that are not related to the underlying DNA sequence of an organism. Learning as of performed by Artificial Neural Networks (ANN) is a process which scope remains limited to an individual lifetime and therefore is Epigenetic.

Application Case-Study

Among the many challenges that form this project, this paper puts focus on the work realized at the platform level for enabling the implementation of distributed applications that take advantage of the bio-inspired features of the platform. This work relies on an agent-based programming methodology and a programming framework made of resident agents which provide support for the three bio-inspired mechanisms presented previously. Three applications are proposed within the project. In this paper we focus on a so-called embodied application (i.e. robotic application) that takes advantage of learning and evolution capabilities of the platform. This case-study application is intended to be a proof of concept of the proposed approach; it uses the pervasive capabilities of the Perplexus platform for it relies on a network of robots that communicate wirelessly. Each robot is equipped with a prototype Ubidule board that handles the bio-inspired mechanisms such as learning, thanks to a software artificial neural network, and evolution that helps improving robot behaviour over time.

The FIPA Multi-agent System

The distributiveness and capabilities of the Perplexus platform greatly rely on the Agent-oriented programming (AOP) style. AOP derives from the initial theory of agent orientation which was first proposed by Yoav Shoham [START_REF] Shoham | Agent Oriented Programming[END_REF]. Agentorientation was initially defined for promoting a social view of computing and finds natural applications in areas such as Artificial Intelligence or social behaviours modeling. This programming format infers object-oriented programming by endowing objects with additional characteristics; they are viewed as entities which exhibit behaviours, capabilities and are entitled to take decisions.

In an AOP systems Agents exclusively communicate with each other through typed messages. The type of these messages vary very much depending on the nature of the exchanged information. Agents may be informing, requesting, offering, accepting, and rejecting requests, services or any other type of information. AOP furthermore sets constraints on the parameters defining the state of the Agent (beliefs, commitments and choices). These constraints actually define the Agent Oriented computational system which is then viewed as a set of communicating software modules that exhibit a certain degree of awareness. These characteristics naturally geared the PERPLEXUS modeling framework towards AOP which fits perfectly the objectives of the platform.

There exists various Multi-Agent Platforms that allow to develop agent based applications, such as JADE [START_REF] Bellifemine | Developing Multi-Agent Systems with JADE[END_REF], JXTA, FIPA OS, JAX or MADKIT. JADE is a Java environment and is therefore portable, FIPA 2 compliant and exists in a lightweight version called LEAP for Light Extensible Agent Platform [START_REF] Lawrence | LEAP into Ad-Hoc Networks[END_REF] that suits our hardware restrictions (PXA270) as it is running with Sun Microsystem J2ME (Java2 Micro Edition) JAVA Virtual Machine. Agents in a JADE Framework "live" in containers. Both containers and Agent Platforms (APs) may span several hardware hosts. Agents existing inside a given container may very well communicate with agents of other containers and even platforms, through a TCP/IP network using the Message Transport Protocols (MTP).

Perplexus Agent Framework

The here presented Bio-inspired Agent Framework is built on top of a JADE agent platform that has been ported to the Ubidule sotware environment. Figure 2 shows the proposed framework; shaded components correspond to platform specific features while others are part of the JADE environment. Hence, the Directory 1. Agent-level, that facilitate various platform-level operations that may prove useful or necessary for some applications. Among the proposed services, event collection, peer discovery, timestamping or semaphore resource access find natural use in many cases. 2. Service level, that take care of ensuring communication reliability in close interaction with the low-level layer discussed below. Provided services aims at answering various needs such as keeping an up-to-date list of accessible nodes in the platform, resolving peers names (nameservice) or gathering information from the physical network topology such as identifying neighbours, etc.

The Ubidules forming an ad-hoc network of mobile nodes, some dedicated mechanisms are mandatory for ensuring network reliability. These mechanisms are integrated in the lower-level layer that implements the Optimized Source Link State Routing (OLSR [START_REF] Jacquet | Optimized link State Routing protocoll for Ad-Hoc Networks[END_REF]). This adaptive routing protocol dynamically changes routes for maintaining the logical network structure that allow any node to contact any other (multi-hop routing). The network agent makes exclusive use of some services of OLSR for retrieving physical network information; for some of its services such as neighbours discovery.

Validation Application

In order to prove the reliability of the platform, a simple validation application based on a race of robotic toys was developed (Fig. 3). Robots which participate use learning (Epigenesis) for improving their performance which is expressed as lap time. Furthermore, the robot population is evolved (Phylogenesis) over multiple generations in order to obtain a species that is better adapted to the problem. The individuals have an average lifetime that is set by the application; whenever an individual dies a new one resulting from the Phylogenetic process is grown (Ontogenesis) onto the same robot. Each robot is equipped with four infrared sensors required for the perception of its environment; these sensors are well-arranged on the four sides of the robot and compute the distance with the nearest obstacle. The robots are capable of executing four basic motions that can be combined: moving forward (MF), moving backward (MB), turning right (TR) or turning left (TL); the Interface agent I provides methods to access to these sensors and actuators. Later, an Ubicom agent U will be implemented in order to provide methods dedicated to the management of the bio-inspired circuit. In our application, the population is composed of three robots, the genome (Fig. 3) of a robot is composed of the values of four thresholds which are used to convert the analog information issued by the infrared sensors into a binary word representation; the data provided by the infrared sensors are compared to the values of the threshold contained in a genome, if a value is lower than the specified threshold, then an obstacle is detected otherwise it is not. The fitness function reports the quality of the expected solution; in our application the fitness is represented by the lap time. Figure 3 shows a schematic representation of the robot controller we have implemented. The output of the system which issues the motion commands (MF, MB, TR, TL) are computed by a software ANN which receives as inputs the four boolean informations of the proximity sensors. This controller represents the behaviour of the Epigenetic agent of our multi agent platform. The genomic description of an individual is made of an array containing all threshold values that are used by the sensors (expresses sensorial abilities) and the specificities of the ANN. In our application, the genome is coded as a 48-bit string as described in Figure 3. A demonstration video is available online at: http://www.lirmm.fr/~brousse/Ubibots.

Conclusion

This paper presents some of the work realized within the confines of the Perplexus Project. These contributions provide a reliable middleware to support both bio-inspiration and agent-oriented programming for distributed pervasive platforms. The proposed framework has been designed with broader application fields in mind and should prove appropriate for many sensor-network applications where adaptability brings advantages. This software environment has been tested through a POE (Phylogenesis, Ontogenesis, Epigenesis) robotic application that proves the reliability of the platform with satisfying result concerning communication and basic POE computation. Future work relies on enabling transparent migration of some software agents (JADE agents) to the bio-inspired integrated circuit in order to increase the performances of the platform.

Fig. 1 .

 1 Fig. 1. Perplexus platform overview

2Fig. 2 .

 2 Fig. 2. Bio-inspired Framework

Fig. 3 .

 3 Fig. 3. Sofware application example

[IST-2006-034632] STREP IST project web site: http://www.perplexus.org

Acknowledgements

This project is funded by the Future and Emerging Technologies programme IST-STREP of the European Community, under grant IST-034632 (PERPLEXUS). The information provided is the sole responsibility of the authors and does not reflect the Communitys opinion. The Community is not responsible for any use that might be made of data appearing in this publication.


Quentin Clairon
email: quentin.clairon@u-bordeaux.fr

Chloé Pasin

Irene Balelli

Rodolphe Thiébaut

Mélanie Prague

Parameter estimation in nonlinear mixed effect models based on ordinary differential equations: An optimal control approach

Keywords: Dynamic population models, Ordinary differential equations, Optimal control theory, Mechanistic models, Nonlinear mixed effects models, Clinical trial analysis

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.

Parameter estimation in nonlinear mixed effect models based on ordinary differential equations: an optimal control approach 1 Introduction

Ordinary differential equation (ODE) models are standard in population dynamics, epidemiology, virology, pharmacokinetics, or genetic regulatory network analysis since they can describe the main mechanisms of interaction between different biological components of complex systems and their evolution over time and because they also provide reasonable approximations to stochastic dynamics [START_REF] Perelson | Hiv-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time[END_REF]; [START_REF] Engl | Inverse problems in systems biology[END_REF]; Villain et al ( 2019)).

For experimental designs with a large number of subjects and a limited number of individual measurements, nonlinear mixed-effects models may be more relevant than single-subject models, since they allow to collect informations from the entire population while accounting for variability among individuals. For instance, clinical trials and pharmacokinetics/pharmacodynamics studies often fall into this category (M. [START_REF] Lavielle | Maximum likelihood estimation of long terms hiv dynamic models and antiviral response[END_REF]; [START_REF] Guedj | Maximum likelihood estimation in dynamical models of hiv[END_REF]). Formally, we are interested in a population where the dynamic of each subject i ∈ 1, n is modeled by the d-dimensional ODE: ẋi (t) = f θ,bi (t, x i (t))

x i (0) = x i,0 (1)
where f is a d-dimensional vector field, θ is a p-dimensional parameter, b i ∼ N (0, Ψ) is a q-dimensional random effect where Ψ is a variancecovariance matrix, x i,0 is the initial condition for subject i. We denote X θ,bi,xi,0 the solution of (1) for a given set (θ, b i , x i,0 ). In (1), we can also incorporate covariate functions z i which are omitted here for the purpose of clarity.

Our goal is to estimate the true population parameters (θ * , Ψ * ) as well as the true subject specific realizations {b * i } i∈ 1, n from partial and noisy observations coming from n subjects and described by the following observational model:

y ij = CX θ * ,b * i ,x * i,0 (t ij ) + ϵ ij .
For the i-th subject, we denote t ij its j-th measurement time-point on the observation interval [0, T ] and n i its total number of available measurements. Here C is a d o ×d sized observation matrix emphasizing the potentially partially observed nature of the process and ϵ ij ∼ σ * × N (0, I d o ) is the measurement error. The vector y i = {y ij } j∈ 1, ni corresponds to the set of observations available for the i-th subject and y = {y i } i∈ 1, n is the set of all observations in the population. We also assume that only a subset x k * i,0 of x * i,0 is perfectly known, the other ones, denoted x u * i,0 , being unknown and they are ordered as follows x i,0 =

x u i,0

T , x k i,0

T T for the sake of clarity. Nonetheless, pre-existing information can be available for x u * i,0 under the form of a priori distribution with a possibly parameter dependent density P x u i,0 | θ, Ψ, b i . The same holds for (θ, Ψ) for which a priori distribution P (θ, Ψ) can be available.

Our problem belongs to the class of parameter estimation problem in nonlinear mixed effect models based on ODEs (NLME-ODEs). In this context, frequentist methods based on likelihood maximization (via different numerical procedures: Laplace approximation [START_REF] Pinheiro | Approximations to the loglikelihood function in the nonlinear mixed effects model[END_REF], Gaussian quadrature [START_REF] Guedj | Maximum likelihood estimation in dynamical models of hiv[END_REF] or SAEM [START_REF] Comets | Parameter estimation in nonlinear mixed effect models using saemix, an r implementation of the saem algorithm[END_REF]; Lavielle and Mentré (2007)) and Bayesian ones aiming to reconstruct the a posteriori distribution or to derive the maximum a posteriori estimator (via MCMC algorithms [START_REF] Lunn | Winbugs -a bayesian modelling framework: Concepts, structure and extensibility[END_REF]; [START_REF] Huang | A bayesian approach to joint mixed-effects models with a skew normal distribution and measurement errors in covariates[END_REF], importance sampling [START_REF] Raftery | Estimating and projecting trends in hiv/aids generalized epidemics using incremental mixture importance sampling[END_REF], approximation of the asymptotic posterior distribution [START_REF] Prague | Nimrod: A program for inference via a normal approximation of the posterior in models with random effects based on ordinary differential equations[END_REF]) have been proposed. In particular, dedicated methods/softwares using the structure of ODE models have been implemented to increase numerical stability and speed up convergence rate [START_REF] Tornoe | Non-linear mixed-effects pharmacokinetic/pharmacodynamic modelling in nlme using differential equations[END_REF]), to reduce the computational time [START_REF] Donnet | Estimation of parameters in incomplete data models defined by dynamical systems[END_REF]) or to avoid the repeated model integration and estimation of initial conditions (Wang et al (2014)). However, all the preceding methods face similar pitfalls due to specific features of population models based on ODEs (with the exception of Wang et al ( 2014)):

1. They do not take into account the presence of model misspecification, a common feature of ODE models used in biology. Indeed, the ODE modeling process suffers from model inadequacy, understood as the discrepancy between the mean model response and the real process, as well as residual variability subject to specific stochastic perturbations or missing elements that disappear by averaging over the entire population [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. As examples of the causes of model inadequacies, one can think of the ODE models used in epidemiology and virology, which are derived by approximations in which, for example, interactions are modeled by pairwise products, while higher order terms and/or the influence of unknown/unmeasured external factors are neglected. As for residual variability, recall that biological processes are often stochastic and the justification for deterministic modeling lies in the approximation of stochastic processes [START_REF] Kurtz | Strong approximation theorems for density dependent markov chains[END_REF]; [START_REF] Kampen | Stochastic Process in Physics and Chemistry[END_REF]). Moreover, in the context of NLME-ODEs, new sources of model uncertainties emerge. Firstly, error measurement in covariates z i can lead to use a proxy function z i instead of z i [START_REF] Huang | A bayesian approach to joint mixed-effects models with a skew normal distribution and measurement errors in covariates[END_REF]). Secondly, the sequential nature of most inference methods leads to estimate {b * i } i∈ 1, n based on an approximation θ instead of θ * . Thus, the structure of mixed-effect models spreads measurement uncertainty into the mechanistic model structure during the estimation. It turns classical statistical uncertainties into model error causes. Estimation of θ * , Ψ * and Parameter estimation in NLME-ODEs: an optimal control approach {b * i } i∈ 1, n must be performed in the presence of the model error, although it is known to dramatically affect the accuracy of methods that do not take it into account [START_REF] Brynjarsdottir | Learning about physical parameters: The importance of model discrepancy[END_REF]). 2. They have to estimate or make assumptions on the true unknown initial conditions x u * i,0 . In ODE models, the initial conditions are generally nuisance parameters in the sense that inferring their values does not bring answers to the scientific questions which motivate the model construction but is necessary for the estimation of the relevant parameters. For example, partially observed compartmental models used in pharmacokinetics/pharmacodynamics often involve unknown initial conditions which needs to be inferred to estimate the transmission rates between compartments, which are the true parameters of interest. Unknown initial conditions imply either assumptions on x u * i,0 values (M. Lavielle and Mentre ( 2011)), another potential cause of model misspecifications, or their estimation [START_REF] Huang | Modeling long-term longitudinal hiv dynamics with application to an aids clinical study[END_REF]). This latter case requires a priori knowledge to derive P x u i,0 | θ, Ψ, b i expression and simultaneous inference of (b * i , x u * i,0 ) as subject specific parameters. This increases the complexity of the related optimization problem and can degrade estimation accuracy. 3. They can face accuracy degradation when the inverse problem of parameter estimation is ill-posed [START_REF] Engl | Inverse problems in systems biology[END_REF]) due to practical identifiability issues. Ill-posedness in ODE models is often due to the geometry induced by the mapping (θ, b i , x i,0 ) -→ CX θ,bi,xi,0 , where there can be a small number of relevant directions of variation skewed from the original parameter axes [START_REF] Gutenkunst | Universally sloppy parameter sensitivities in systems biology models[END_REF]). This problem, called sloppiness, often appears in ODE models used in biology [START_REF] Gutenkunst | Universally sloppy parameter sensitivities in systems biology models[END_REF]; [START_REF] Leary | Computational models in the age of large datasets[END_REF]) and leads to an ill-conditioned Fisher Information Matrix. For maximum likelihood estimators this is a cause of high variance due to the Cramér-Rao bound. For Bayesian inference, it leads to a nearly singular asymptotic posterior distribution because of Bernsteinvon Mises theorem (see [START_REF] Campbell | Bayesian collocation tempering and generalized profiling for estimation of parameters from differential equation models[END_REF] for the computational induced problems). Although this problem is in part mitigated in NLME-ODEs which merge different subjects for estimating (θ * , Ψ * ) and use distribution of b i | Ψ as prior at the subject level [START_REF] Lavielle | What do we mean by identifiability in mixed effects models? Journal of pharmacokinetics and pharmacodynamics Lavielle M, Mentré F[END_REF]), estimation accuracy can benefit from the use of regularization techniques. These specific features of ODE-based population models limit the amount of information classic approaches can extract for estimation purposes from observations no matter their qualities or abundances. This advocates for the development of new estimation procedures. Approximate methods (Varah (1982); [START_REF] Ramsay | Parameter estimation for differential equations: A generalized smoothing approach[END_REF]) have already proven to be useful for ODE models to face such issues. They rely on an approximation of the solution of the original ODE (1) which is expected to have a smoother dependence with respect to the parameters and to relax the constraint imposed by the model during the estimation process. The interest of such approximations is twofold. Firstly, they produce estimators with a better conditioned variance matrix comparing to classic likelihood based approaches. Secondly, they reduce the effect of model error on estimator accuracy. Also, some of these approximations bypass the need to estimate initial conditions [START_REF] Ramsay | Parameter estimation for differential equations: A generalized smoothing approach[END_REF]; [START_REF] Clairon | A regularization method for the parameter estimation problem in ordinary differential equations via discrete optimal control theory[END_REF]). Still, these methods are currently limited to cases where observations are coming from one subject.

In this work, we develop a new estimation method adapted to NLME-ODEs integrating such approximations to mitigate the effect of model misspecification and poorly identifiable parameters on estimation accuracy, while avoiding the need to estimate x u * i,0 as additional subject specific parameters. We propose here a hierarchical profiling approach taking the form of a nested estimation procedure. The unknowns initial conditions x u * i,0 i∈ 1,n are seen as nuisance parameters for {b * i } i∈ 1,n estimation, which are in turn considered as nuisance parameters for population parameter (θ * , Ψ * , σ * ) estimation. This lead to the construction of outer, middle and inner criteria for the estimation of (θ * , Ψ * , σ * ), {b * i } i∈ 1,n and x u * i,0 i∈ 1,n respectively. The inner criteria is designed to incorporate P x u i,0 | θ, Ψ, b i if an expression is proposed for it but can also be defined without if no prior information exists for x u * i,0 i∈ 1,n . Also, this criterion accounts for model error presence by assuming that the actual dynamic of each subject is better described by a perturbed version of the ODE (1). This added perturbation aims to capture different sources of errors at the subject level [START_REF] Brynjarsdottir | Learning about physical parameters: The importance of model discrepancy[END_REF]; Tuo and Wu (2015)). We control the magnitude of the acceptable perturbations by defining the inner criteria through a cost function balancing the two contrary objectives of fidelity to the observations and to the original model: to this end, we introduce a model discrepancy penalization term. The practical computation of the chosen perturbations requires to solve optimal control problems [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF]) known as tracking problems. This is done using a method inspired by [START_REF] Cimen | Global optimal feedback control for general nonlinear systems with nonquadratic performance criteria[END_REF] which has the advantage to automatically provide an estimator for x u * i,0 with no additional computational costs. This is the key element to efficiently profile on initial conditions during b * i estimation. In section 2, we present the inner, middle and outer criteria used to define our estimator. In section 3, we compare our approach with classic maximum likelihood in simulations. Then, we proceed to the real data analysis coming from clinical studies and a model of the antibody concentration dynamics following immunization with an Ebola vaccine in East African participants [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]). Section 5 concludes and discuss future extensions of the method.

2 Estimator construction: definition of the inner/middle/outer criteria

From now on, we use the Choleski decomposition σ 2 Ψ -1 = △ T △ and the parametrization ϕ := (θ, ∆, σ) instead of (θ, Ψ, σ) to enforce positiveness Parameter estimation in NLME-ODEs: an optimal control approach and symmetry of Ψ and denote in a summarized way the set of all population parameters. The norm ∥.∥ 2 denotes the classic Euclidean one defined by ∥b∥ 2 = √ b T b. Similarly as in the Expectation-Maximization (EM) algorithm, we estimate the population and individual parameters via a nested procedure:

Estimator ϕ obtained by minimization of an outer criterion F based on an approximation of min b min x u 0 (-ln P [ϕ, b, x u 0 | y]), the log joint-distribution of (ϕ, b, x u 0 ) sequentially profiled on b := {b i } i∈ 1,n and x u 0 := x u 0,i i∈ 1,n , which are respectively the set of all random effects and unknown initial conditions among all subjects. Estimator b i := b i (ϕ) obtained for each subject i by minimization of a middle criterion G i based on an approximation of min x u 0,i -ln P(y i , b i , x u 0,i | ϕ) , the log joint-distribution of the data, the random effects and unknown initial conditions profiled on the latter. Estimator x u 0,i := x u 0,i (ϕ, b i ) obtained for each subject i by minimization of an inner criterion H i based on an approximation of -ln P(y i , x u 0,i | ϕ, b i ), the log joint-distribution of the data and unknown initial conditions.

Our estimation procedure can be expressed in a pseudo-algorithmic way.

1/ Outer criteria minimization:

ϕ = arg min ϕ F (ϕ) := arg min ϕ -2 ln P(ϕ, b, x u 0 | y),
for a given subject i and ϕ value:

2/ middle criteria minimization:

b i (ϕ) = arg min bi G i (b i | ϕ) := arg min bi -2 ln P(y i , b i , x u 0,i | ϕ),
for a given b i value:

3/ inner criteria minimization:

x u 0,i (ϕ, b i ) = arg min x u 0,i H i (x 0,i | ϕ, b i ) := arg min x u 0,i -2 ln P(y i , x u 0,i | ϕ, b i ).
In the following sections, we derive the expressions of F, G i and H i starting with H i since each criterion construction rely on lower level ones. Finally, despite that the following formal presentation of criteria are made for any P x u i,0 | ϕ, b i expressions, we have to restrict ourselves to uniform, normal and log-normal densities in practice to use our numerical procedures.

Inner criterion

In this section, we construct the criteria H i used to estimate x u * 0,i for a given (ϕ, b i ) value. A classic procedure would lead to jointly estimate b * i , x u * 0,i by maximization of the log joint-likelihood function of the data and (b i , x u 0,i ). However for each subject, we want to:

1. profile on x u * 0,i during random effects estimation to limit b * i estimation degradation due to presence of nuisance parameters, 2. use prior knowledge given by P x u i,0 | ϕ, b i if available, 3. allow an acceptable deviation from the assumed model at the population level to take into account possible model misspecifications. To solve the first and second point, we define our estimator:

1. as the maximizer of the joint conditional likelihood P(

y i , x u 0,i | ϕ, b i ) if P x u i,0 | ϕ, b i is available, 2. otherwise as the maximizer of P(y i | ϕ, b i , x u 0,i ). Since P(y i , x u 0,i | ϕ, b i ) = P(y i | ϕ, b i , x u 0,i )P x u i,0 | ϕ, b i , we have arg max x u 0,i P(y i | ϕ, b i , x u 0,i ) = arg max x u 0,i P(y i , x u 0,i | ϕ, b i ) if P x u i,0 | ϕ, b i is constant.
So, the estimation criteria in absence of prior information is equivalent to choosing a uniform prior over x u 0,i space and constitute only a particular case. We will thus focus on P(y i , x u 0,i | ϕ, b i ) from now one. We have:

P(y i , x u 0,i | ϕ, b i ) = j P(y ij | ϕ, b i , x u 0,i )P x u i,0 | ϕ, b i = j (2π) -d o /2 σ -d o e -0.5∥CX θ,b i ,x 0,i (tij )-yij ∥ 2 2 /σ 2 P x u i,0 | ϕ, b i ,
from which we derive the joint likelihood estimator:

x u 0,i (ϕ, b i ) = arg min x u 0,i -2 ln P(y i , x u 0,i | ϕ, b i , x u 0,i ) = arg min x u 0,i 1 σ 2 j CX θ,bi,x0,i (t ij ) -y ij 2 2 -2 ln P x u i,0 | ϕ, b i .
We also want to allow the presence of perturbations at the subject scale comparing to the original model defined at the population level. For this, we assume the regression function is no longer X θ,bi,x0,i , but rather X θ,bi,xi,0,ui , the solution of:

ẋi (t) = f θ,bi (t, x i (t)) + Bu i (t) x i (0) = x i,0 . (2)
This perturbed ODE has been obtained by the addition of the forcing term t → Bu i (t) to ODE (1) with B a d×d u matrix and u i a function in L 2 [0, T ] , R du representing the perturbation. However, to ensure the possible perturbations remain small, we replace the data fitting criterion

j CX θ,bi,x0,i (t ij ) -y ij 2 2 by min ui C i (x u i,0 , u i | θ, b i , U ),
where

C i (x u i,0 , u i | θ, b i , U ) = j CX θ,bi,x0,i,ui (t ij ) -y ij 2 2 + ∥u i ∥ 2 U,L 2 ,
Parameter estimation in NLME-ODEs: an optimal control approach

and ∥u i ∥ 2 U,L 2 = T 0 u i (t) T U u i (t)dt is the weighted Euclidean norm.
Here, the magnitude of the allowed perturbations is controlled by a positive definite and symmetric weighting matrix U . Finally, we obtain:

x u 0,i (ϕ, b i ) := arg min x u 0,i H i (x u 0,i | ϕ, b i ) (3)
where

H i (x u 0,i | ϕ, b i ) = min x u 0,i 1 σ 2 min ui C i (x u i,0 , u i | θ, b i , U ) -2 ln P x u i,0 | ϕ, b i . Computing H i (x u 0,i | ϕ, b i ) requires to solve the infinite dimensional opti- mization problem min ui C i (x u i,0 , u i | θ, b i , U ) in L 2 [0, T ] , R du .
This problem belongs to the field of optimal control theory for which dedicated approaches have been developed [START_REF] Sontag | Mathematical Control Theory: Deterministic finitedimensional systems[END_REF]; [START_REF] Aliyu | Nonlinear H-Infinity Control, Hamiltonian Systems and Hamilton-Jacobi Equations[END_REF]; [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF]). Here we use the same method as in [START_REF] Clairon | A regularization method for the parameter estimation problem in ordinary differential equations via discrete optimal control theory[END_REF] which is detailed in Appendix A. All it requires from the user is to specify a pseudo-linear representation of ODE (1), i.e., a possibly state-dependent matrix A θ,bi (t, x i (t)) and state-independent vector r θ,bi (t) such that:

f θ,bi (t, x i (t)) = A θ,bi (t, x i (t)) x i (t) + r θ,bi (t). ( 4
)
This formulation is crucial for solving the optimal control problem in a computationally efficient way. Linear models already fit in this formalism with A θ,bi (t) := A θ,bi (t, x i (t)). For nonlinear models, the pseudo-linear representation is not unique but always exists (in order to exploit this non-uniqueness as an additional degree of freedom, see [START_REF] Cimen | State-dependent riccati equation (sdre) control: A survey[END_REF] section 6). This method presents the advantage of formulating min ui C i (x u i,0 , u i | θ, b i , U ) as a quadratic form (or a sequence of quadratic forms) with respect to x u 0,i . Thus, if we choose a uniform, normal or log-normal law for P

x u i,0 | ϕ, b i , arg min x u 0,i H i (x u 0,i | ϕ, b i )
has a closed form expression (approximated for log-normal), and obtaining x u 0,i (ϕ, b i ) does not add any computational complexity comparing to

min ui C i (x u i,0 , u i | θ, b i , U ). For a given x u 0,i , the perturbation u i corresponding to the solution of min ui C i (x u i,0 , u i | θ, b i , U
) is named optimal control and denoted u i,ϕ,bi,x u 0,i . In particular, we denote u i,θ,bi := u i,θ,bi, x0,i the optimal control corresponding to the initial condition estimator x 0

,i = x u 0,i (ϕ, b i ) T , x k i,0 T T
. The solution of (2) corresponding to the optimal control u i := u i,θ,bi is denoted X θ,bi and named optimal trajectory: this will be considered as the regression function for the i-th subject. X θ,bi is thus defined as solution of ODE (2) which needs the smallest perturbation in order to get close to the observations. In particular, X θ,bi and u i,θ,bi are respectively the subject specific state variable and perturbation such that:

H i ( x u 0,i (ϕ, b i ) | ϕ, b i ) = 1 σ 2 j CX θ,bi (t ij ) -y ij 2 2 + ∥u i,θ,bi ∥ 2 U,L 2 -2 ln P( x u 0,i (ϕ, b i ) | ϕ, b i ).
(5)

Again, formal expressions can be derived for both u i,θ,bi and x u 0,i (ϕ, b i ), but they present no interest for the sake of explanation and are left in Appendix A.

Middle criterion

To construct an estimator b i of the random effects, we rely on an approximation of ln P(y i , b i , x u 0,i | ϕ) profiled on the unknown initial conditions. Since

P(y i , b i , x u 0,i | ϕ) = P(y i | ϕ, b i , x u 0,i )P(b i , x u 0,i | ϕ) = P(y i | ϕ, b i , x u 0,i )P(x u 0,i | ϕ, b i )P(b i | ϕ), with P(b i | ϕ) the density of b i ∼ N (0, σ 2 △ T △ -1
), we can define as estimator:

b i (ϕ) = arg min bi min x u 0,i -2 ln P(y i , b i , x u 0,i | ϕ) = arg min bi min x u 0,i 1 σ 2 j CX θ,bi,x0,i (t ij ) -y ij 2 2 -2 ln P x u i,0 | ϕ, b i + ∥∆bi∥ 2 2 σ 2 .
Still, we use the same relaxation & penalization scheme as in the previous section to account for model error presence for b * i estimation. We replace again the term j CX θ,bi,x0,i (t ij ) -y ij 2 2 by min ui C i (x u i,0 , u i | θ, b i , U ) in the previous criteria and we end up with the following estimator:

b i (ϕ) := arg min bi G i (b i | ϕ) (6)
where:

G i (b i | ϕ) = H i ( x u 0,i (ϕ, b i ) | ϕ, b i ) + ∥∆b i ∥ 2 2 σ 2 . ( 7
)

Outer criterion

F general expression

We focus in this section on population parameter estimation. Classic maximum likelihood based approaches generally consider as estimator:

ϕ := E b1 [P(ϕ, b 1 | y)].
That is, they get rid of the unknown subject specific parameters by taking the mean value of

P(ϕ, b 1 | y) where b 1 ∼ N (0, σ 2 △ T △ -1 ).
This generally requires the numerical approximation of integrals of possibly high dimensions (the same as b 1 ), a source of approximation and computational issues [START_REF] Pinheiro | Approximations to the loglikelihood function in the nonlinear mixed effects model[END_REF]).

To avoid this, we consider the random Parameter estimation in NLME-ODEs: an optimal control approach effects as nuisance parameters and rely on a classic profiling approach for ϕ * estimation [START_REF] Murphy | On profile likelihood[END_REF]). Instead of taking the mean, we rely on the profiled joint distribution sequentially with respect to b := {b i } i∈ 1,n and x u 0 = x u 0,i i∈ 1,n , or equivalently min b min x u 0 (-2 ln P(ϕ, b, x u 0 | y)). Bayes formula gives us P(ϕ, b, x u 0 | y) ∝ P(y, b, x u 0 | ϕ)P(ϕ) and we get P(ϕ, b, x u 0 | y) ∝ i P(y i , b i , x u 0,i | ϕ) P(ϕ) by conditional independence of subject by subject observations and subject specific parameters. It follows that

min b min x u 0 (-2 ln P [ϕ, b, x u 0 | y]) ∝ i min bi min x u 0,i -2 ln P(y i , b i , x u 0,i | ϕ) -2 ln P(ϕ),
from which we derive the estimator

ϕ = arg min ϕ i min bi min x u 0,i 1 σ 2 j CX θ,bi,x i 0 (t ij ) -y ij 2 2 -2 ln P(x u 0,i | ϕ, b i ) + ∥∆bi∥ 2 2 σ 2 + (d o i n i + nq) ln σ 2 -2n ln|△| -2 ln P(ϕ)
by using the exact expression of ln P(y i , b i , x u 0,i | ϕ) (computational details are recalled in Appendix B). In order to account for the presence of model error and limit its effect on estimation, we replace in the last expression the classic profiled likelihood estimator for b * i and x u * 0,i by b i (ϕ) and x u 0,i (ϕ, b i ) respectively and X θ,bi,x0,i by X θ,bi . This leads us to the following population parameter estimator:

ϕ := arg min ϕ F (ϕ) , (8)
where:

F (ϕ) = 1 σ 2 i j CX θ, bi(ϕ) (t ij ) -y ij 2 2 + ∆ b i (ϕ) 2 2 -2 ln P( x u 0,i (ϕ, b i (ϕ)) | ϕ, b i (ϕ)) + (d o i n i + nq) ln σ 2 -2n ln|△| -2 ln P(ϕ) . (9) 2.3.2 F profiling on σ for uniform x u 0,i distribution If P(x u 0,i | ϕ, b i ) is constant then x u 0,i (ϕ, b i ) and b i (ϕ) do not depend on σ i.e. x u 0,i (ϕ, b i ) = x u 0,i (θ, b i ) and b i (ϕ) = b i (θ, ∆
) and consequentially neither does X θ, bi(θ,∆) . So, for each (θ, ∆), the maximizer in σ 2 of F has a closed form expression:

σ 2 (θ, ∆) = 1 (d o i n i + qn) i j CX θ, bi(θ,∆) (t ij ) -y ij 2 2 + ∆ b i (θ, ∆) 2 2 . ( 10
)
By using σ 2 (θ, ∆) expression, we get min

σ 2 F (θ, ∆, σ | y) = F [θ, ∆ | y] where: F [θ, ∆ | y] = d o i n i + qn ln σ 2 (θ, ∆) -2n ln|△| -2 ln P(ϕ).
Thus, we can profile F on σ 2 and define our estimator as:

θ, ∆ = arg min (θ,∆) F [θ, ∆ | y] . (11)
An estimator of σ * is obtained from there by computing σ 2 θ, ∆ , given by equation ( 10). The details of F derivation are left in appendix B.

Asymptotic Variance-Covariance matrix estimator for θ, ∆

We derive an estimator of the asymptotic variance of θ, ∆ . Here we restrict to the case described in section 2.3.2 when a uniform distribution is chosen for x u i,0 and the outer criterion is profiled on σ. The general case can be considered similarly, but we withdraw it for the sake of clarity since it is not used in following simulation works. We highlight that in practice the matrix ∆ is parametrized by a vector δ of dimension q ′ , i.e △ := △(δ) and we give here a variance estimator of θ, δ . From this, the variance of ∆ can be obtained using classic delta-methods (see van der Vaart (1998) chapter 3). Conditions on model structural identifiability and regularity are required to derive the existence of this asymptotic variance, we precise such sufficient conditions in appendix D.

Theorem 1 There is a model dependent lower bound λ such that if ∥U ∥ 2 > λ then the estimator θ, δ converges almost surely to a constant value θ, δ such that:

√ n( θ -θ, δ -δ) ⇝ N 0, A(θ, δ) -1 B(θ, δ) A(θ, δ) -1 T , where A(θ, δ) = limn 1 n n i=1 ∂ J(θ,δ,y i ) ∂(θ,δ) , B(θ, δ) = limn 1 n i J(θ, δ, y i ) J(θ, δ, y i ) T and the vector valued function J(θ, δ, y i ) = J θ (θ, δ, y i ) J δ (θ, δ, y i )
is given by:

J θ (θ, δ, y i ) = d dθ h( b(θ, ∆(δ)), θ, ∆(δ), y i ) J δ (θ, δ, y i ) = d dδ h( b i (θ, ∆(δ)), θ, ∆(δ), y i ) - 2 d o E[n1]+q T r △(δ) -1 ∂△(δ) ∂δ k h( b i (θ, ∆(δ)), θ, ∆(δ), y i ) , where h(b i , θ, ∆, y i ) = ∥∆b i ∥ 2 2 + j CX θ,bi (t ij ) -y ij 2 2 .
Parameter estimation in NLME-ODEs: an optimal control approach

The proof is left in appendix D. The practical interest of this theorem is to give an estimator of Variance-Covariance

V ( θ, δ) ≃ A( θ, δ) -1 B( θ, δ) A( θ, δ) -1 T /n with the matrices A( θ, δ) = -1 n n i=1 ∂J( θ, δ,y i ) ∂(θ,δ)
and B = 1 n n i=1 J( θ, δ, y i )J( θ, δ, y i ) T . The (p + q) components of the vector valued function J for 1 ≤ k ≤ p are given by :

J k (θ, δ, y i ) = d dθ k h( b(θ, ∆(δ)), θ, ∆(δ), y i ),
and for p + 1 ≤ k ≤ p + q by

J k (θ, δ, y i ) = d dδ k h( b i (θ, ∆(δ)), θ, ∆(δ), y i ) - 2n d o i ni+qn T r △(δ) -1 ∂△(δ) ∂δ k h( b i (θ, ∆(δ)), θ, ∆(δ), y i ).
Now that we have proven the existence of the variance-covariance matrix V (θ, δ) such that δ -δ ⇝ N 0, V (θ, δ) , we can use the Delta method to derive the asymptotic normality of the original matrix

Ψ δ = σ 2 ∆( δ) T ∆( δ) -1
as well as an estimator of its asymptotic variance. In the case of a diagonal matrix Ψ, composed of the elements Ψ 2 1 , . . . Ψ 2 q and of the parametrization △(δ) = Diag( e δ l l∈ 1,q ) used in section 3, we derive:

   Ψ 1 ( δ) . . . Ψ q ( δ)   -    Ψ 1 (δ * ) . . . Ψ q (δ * )    ⇝ N   0, σ 2    e -δ * 1 0 0 0 . . . 0 0 0 e -δ * q    V (θ * , δ * )    e -δ * 1 0 0 0 . . . 0 0 0 e -δ * q       .
Remark 1 The previous theorem 1 states that we retrieve a parametric convergence rate. Thus, we avoid the pitfall described in [START_REF] Sartori | Modified profile likelihood in models with stratum nuisance parameters[END_REF] for profiled methods in presence of a number of nuisance parameters increasing with the number of subjects (or strata to resume [START_REF] Sartori | Modified profile likelihood in models with stratum nuisance parameters[END_REF] terminology) potentially leading to bias accumulation for score functions among subjects. The i.i.d structure of random effects allows us to rely on central limit theorem to avoid this accumulation phenomenon.

Results on simulated data

We compare the accuracy of our approach with maximum likelihood (ML) in different models and experimental designs reflecting the problems exposed in the introduction, that is estimation in 1/presence of model error, 2/partially observed framework with unknown initial conditions and 3/presence of poorly identifiable parameters. We proceed to Monte-Carlo simulations based on N M C = 100 runs. At each run, we generate n i observations coming from n subjects on an observation interval [0, T ] with Gaussian measurement noise of standard deviation σ * . From these data, we estimate θ * , Ψ * and b * i with both estimation methods. We quantify the accuracy of each estimator ψ p of the population parameters estimate ψ = θ, Ψ via Monte-Carlo computation of the bias Bias( ψ

p ) = E ψ p -ψ * p , the empirical variance V e ( ψ p ) = E E ψ p -ψ * p 2
, the mean squared error M SE( ψ p ) = Bias( ψ p ) 2 + V e ( ψ p ), the estimated variance V ψ p , as well as the coverage rate of the 95%-confidence interval derived from it. This coverage rate, denoted CR in the following results, corresponds to the frequency at which the interval ψ p ± z 0.975 V ψ p contains ψ * p with z 0.975 the 0.975-quantile of the centered Gaussian law. We compute the previous quantities for the normalized values ψ norm p := ψp ψ * p to make relevant comparisons among parameters with different order of magnitude. For b * i , we estimate the mean squared error

M SE( b i ) = E b * i -b i 2 2
. For each subsequent examples, we give the results

for n = 50 and present in appendix C the case n = 20 to analyze the evolution of each estimator accuracy with respect to data sparsity.

In the following, we use the superscript M L to denote the ML estimator. For the fairness of comparison with ML, we choose a non-informative prior i.e. ln P [θ, ∆] = 0 for our method throughout this section. Also, we do not use a distribution for x u i,0 for our approach. For ML which requires it, we will use the right parametric form for P x u i,0 | ϕ, b i . If the ODE (1) has an analytical solution, the ML estimator is computed via SAEM algorithm (SAEMIX package [START_REF] Comets | Parameter estimation in nonlinear mixed effect models using saemix, an r implementation of the saem algorithm[END_REF]). Otherwise, it is done via a restricted likelihood method dedicated to ODE models implemented in the nlmeODE package (Tornoe et al ( 2004)). For our method, we need to select U balancing model and data fidelity in the inner and middle criteria (5)-( 7). We use the method presented in G. [START_REF] Hooker | Parameterizing state-space models for infectious disease dynamics by generalized profiling: measles in ontario[END_REF] to compute EP i (U ), the prediction error for the subject i corresponding to the estimators θ U , b i,U i∈ 1, n obtained for a given matrix U . From this, we compute EP(U ) = i EP i (U ) the global prediction error for the whole population. We test a trial of weighing matrices {U l } l∈ 1,L and retain the one minimizing EP and denote θ, Ψ, b i, i∈ 1, n the corresponding estimator. For solving the optimization problems required for computing our criteria, we use the Nelder-Mead algorithm implemented in the optimr package [START_REF] Nash | Using and extending the optimr package[END_REF]). All optimization algorithms used here require a starting guess value. We start from the true parameter value for each of them. By doing so, we aim to keep distinct two problems: 1)the numerical stability of the estimation procedures, 2)the intrinsic accuracy of the different estimators. These two problems are correlated, but we aim to adress only the latter which corresponds to the issues raised in introduction. Still, we check on preliminary analysis that local minima presence was not an issue in the neighborhood of (θ * , △ * ) by testing different starting points for all methods. 12) and corresponding observations. Right: Solution of ( 12) and a realization of ( 13) for the same parameter values.

No problem appears for our method and SAEMIX. A negligible number of non convergence cases appear for nlmeODE which have been discarded thanks to the convergence criteria embedded in the package.

Application 1 -Partially observed linear model

We consider the population model where each subject i follows the ODE:

   Ẋ1,i = ϕ 2,i X 2,i -ϕ 1,i X 1,i Ẋ2,i = -ϕ 2,i X 2,i (X 1,i (0), X 2,i (0)) = (x 1,0 , x 2,0,i ) (12)
with the following parametrization: log(ϕ 1,i ) = θ 1 + b i and log(ϕ 2,i ) = θ 2 where b i ∼ N (0, Ψ). The true population parameter values are θ * = (θ * 1 , θ * 2 ) = (log (0.5) , log (2)) and Ψ * = 0.5 2 and we are in a partially observed framework where only X 1,i is accessible. The true initial conditions are distributed with x * 1,0,i ∼ N (2, 0.5) and x * 2,0,i ∼ N (3, 1). An analytic solution exists for ODE (12). In particular the first component is given by X 1,i (t) = e -ϕ1,it (x 1,0 + x2,0ϕ2,i ϕ1,i-ϕ2,i (e (ϕ1,i-ϕ2,i)t -1)) and will be used for estimation with the SAEMIX package. We generate n i = 11 longitudinal observations per subject on [0, T ] = [0, 10] with measurement noise of standard deviation σ = 0.05. An example of sampled observations and corresponding solutions are plotted in figure 1. We want to investigate the impact of initial condition, especially the unobserved one x * 2,0,i , on the ML estimator accuracy. Indeed, our method does not need to estimate x * 2,0,i and thus no additional difficulties appear in this partially observed framework. For the ML, however, it is a nuisance subjectspecific parameter that should be estimated and for which no observations are available. For this, we compute θ M L x0 , θ M L x0,2 and θ M L the ML estimator respectively when: 1) both initial conditions are perfectly known, 2) x * 1,0,i is replaced by the measured value, 3) in addition, x * 2,0,i has to be estimated. 12). The different subscripts stand for the following estimation scenarios: 1) x 0 when both initial conditions are set to x * 0,1 , x * 0,2 , 2) x 0,2 when x 0,i is set to y i,0 and x 0,2 to x * 0,2 , 3) absence of subscript when x 0,i is set to y i,0 and x 0,2 is estimated. Results from our method are in bold.

Well-specified case

We used the exact model described in Section 3.1 for the estimation procedure. Thus we are in a completely well-specified setting, with all mechanisms modeled. We present the estimation results in table 1 -left side. For ML, the results are goods in terms of accuracy and consistent in terms of asymptotic confidence interval coverage rate when both initial conditions are known: 95% for θ 1 and θ 2 , which is consistent with theoretical results. However, there is a significant drop in accuracy when x * 2,0,i has to be estimated. In particular, the coverage rate drops to 86% and 80% for θ 1 and θ 2 respectively. Interestingly, ML inaccuracy is driven by bias and under-estimated variance when initial conditions are not known (as shown by a greater V e than V ). In this case our method provides a relevant alternative: it gives accurate estimations with a good coverage rate for all parameters while avoiding the estimation of x * 2,0,i . Variances are properly estimated compared to empirical variances. Estimation of individual random effects is also more accurate with our method, with a MSE for b i 10 times smaller compared to ML with unknown initial conditions.

Misspecified case in presence of model error at the subject level

To mimic the presence of misspecification, we now generate the observations from the hypoelliptic stochastic model:

   dX 1,i = ϕ 2,i X 2,i dt -ϕ 1,i X 1,i dt dX 2,i = -ϕ 2,i X 2,i dt + αdB t (X 1,i (0), X 2,i (0)) = (x 1,0 , x 2,0,i ) (13)
with B t a Wiener process and α = 0.1 the diffusion coefficient. For the sake of comparison, a solution of ( 12) and a realization of its perturbed counterpart given by ( 13) are plotted in figure 1. This framework where stochasticity only Parameter estimation in NLME-ODEs: an optimal control approach affects the unmeasured compartment is known to be problematic for parameter estimation and inference procedures are yet to be developed for sparse sampling case. From figure 1 it is easy to see that the diffusion α will be hard to estimate when we only have observations for X 1,i . Thus, we still estimate the parameters from the model ( 12) which is now seen as a deterministic approximation of the true stochastic process. Still, it is expected that our method will mitigate the effect of stochasticity on the estimation accuracy by taking into account model misspecification. Results are presented in table 1 -right side. The differences between the two methods are similar to the previous well-specified case with an additional loss of accuracy coming from model error for both estimators. However, the misspecification effect for ML is more pronounced comparing to our method which manages to limit the damages done. This illustrates the benefits of taking into account model uncertainty for estimation, in particular here when model error occurs in the unobserved compartment, a situation in which classic statistical criteria for model assessment based on a data fitting criterion are difficult to use.

Application 2 -Partially observed nonlinear model

We consider the model presented in De [START_REF] Gaetano | Mathematical modelling of the intravenous glucose tolerance test[END_REF] for the analysis of glucose and insulin regulation:

   Ġi = S G (G B -G i ) -X i G i İi = γt(G i -h) -m i (I i -I B ) Ẋi = -p 2 (X i + S I (I i -I B )).
(14)

The ODE system (14) rules the behavior of circulating glucose G i and insulin I i in blood as well as insulin X i present in interstitial fluid. We are in a partially observed case where only G i and I i are measured. The values of parameters (p 2 , γ, h, G B , I B ) are set to (-4.93, -6.85, 4.14, 100, 100)

* i,0 = G * 0,i , I * 0,i , X * 0,i are distributed according to ln(x * i,0 ) ∼ N (l x * 0 , Ψ l x * 0
) with l x * 0 = (5.52, 4.88, -7) and Ψ l x * 0 = 0.17 2 , 0.1 2 , 10 -4 . We generate m i = 5 observations on [0, T ] = [0, 180] with Gaussian measurement noise of standard deviation σ * = 3. As in the previous example, we investigate the impact of unknown initial conditions on the estimators accuracy. We are particularly interested in the joint estimation of θ S I , which appears only in the equation ruling the unobserved state variable X i , and x * 0,i required for each subject by ML. For this, we distinguish two cases, 1) when θ S I is known, 2) when θ S I has to be estimated and we denote respectively θ Si and θ the corresponding estimators. Finally, since the model is nonlinear we have to specify a pseudo-linear representation of the

Well-specified Misspecified MSE Bias V e V CR MSE bi MSE Bias V e V CR MSE bi θS G θ M L S i
5e-5 2e-3 4e-5 9e-6 0.95 6e-5 3e-3 6e-5 2e-5 0.85

θ M L
2e-3 0.03 1e-3 8e-5 0.85 2e-3 3e-3 1e-3 2e-4 0.54 θS i 1e-5 4e-4 1e-5 8e-6 0.95 2e-5 -2e-5 2e-5 2e-5 0.93 θ 2e-4 -6e-4 2e-4 2e-4 0.96 3e-4 -1e-3 3e-4 4e-4 0.93

θS I θ M L S i
known known Table 2 Results of estimation for model ( 14). The different subscripts stand for the following estimation scenarios: 1) S i when S i is set to S * i , 2) absence of subscript when S i is estimated. Results from our method are in bold.

θ M L 2e-3 0.
vector field as in (4):

A θ,bi (t, G i , I i , X i ) =   -S G 0 -G i γt -m i 0 0 -p 2 S I -p 2   , r θ,bi (t) =   S G G B -γth + m i I B p 2 S I I B   .

Well-specified case

We present the estimation results in table 2 -left side. Our method has small bias and achieve good coverages in all cases. We obtain smaller MSE than ML and avoid the drop in coverage rate of the confidence interval in the case of θ * S I estimation, which is often needed in practice. The difference between the two estimators behaviors is explained by the fact that they are defined through the construction of two different optimization problems. At the population level, our approach leads to minimize a cost function depending on a 4-dimensional parameter whereas ML, due to its need to estimate x * i,0 , considers a 10-dimensional one. Thus, the parameter spaces explored by each method to look for the minimum are very different.

Misspecified case in presence of model error at the subject level

To mimic misspecification presence, we generate the observations from the stochastic model:

   dG i = (S G (G B -G i ) -X i G i ) dt + α 1 dB 1,t dI i = (γt(G i -h) -m i (I i -I B ))dt + α 2 dB 2,t dX i = (-p 2 (X i + S I (I i -I B ))) dt + α 3 dB 3,t , (15)
where the B i,t are Wiener processes and (α 1 , α 2 , α 3 ) = 2, 2, 2 × 10 -4 their diffusion coefficients. We present the estimation results in table 2 -right side.

Parameters Biological interpretation

Values δL long-lived B-cells declining rate log(2)/(364 × 6)

θ * θ * δ S
Mean log-value for δS , the short-lived cells declining rate log(log(2)/1.2) ≃ -0.54

θ * ϕ S
Mean log-value for ϕS , the antibodies influx from short-lived cells log(2755) ≃ 7.92

θ * ϕ L
Mean log-value for ϕL, the antibodies influx from long-lived cells log(16) ≃ 2.78

θ * δ Ab
Mean log-value for δAb, the antibodies declining rate log(log(2)/24) ≃ -3.54

Ψ * Ψ * ϕ S
Inter individual variance for log(ϕS,i)

0.92 2 Ψ * ϕ L
Inter individual variance for log(ϕL,i)

0.85 2 Ψ * δ Ab
Inter individual variance for log(δAb,i)

0.3 2

Table 3 Biological interpretation and parameter values

For ML, the drop in coverage rate for θ * S G and θ * S I is even more striking when θ * S I needs to be estimated. This is explained by the effect of model misspecification which bias and the fact that ML does not take into account this new source of uncertainty which leads to under-estimation of variance and too narrow confidence intervals. Our method achieved small biais, nominal coverages and small MSE for random effects.

Application 3 -Antibody concentration evolution model

We consider the model presented in [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF] to analyze the antibody concentration, denoted A i , generated by two populations of antibody secreting cells: the short lived, denoted S i , and the long-lived, denoted L i :

       Ṡi = -δ S S i Li = -δ L L i Ȧi = ϑ S,i S i + ϑ L,i L i -δ Ab,i A i (S i (0), L i (0), A i (0)) = (S 0,i , L 0,i , A 0,i ) . ( 16
)
This model is used to quantify the humoral response on different populations after an Ebola vaccine injection with a 2 doses regimen seven days after the second injection when the antibody secreting cells enter in a decreasing phase. These cells being unobserved, the preceding equation can be simplified to focus on antibody concentration evolution:

Ȧi = ϕ S,i e -δ S t + ϕ L,i e -δ L t -δ Ab,i A i ( 17
)
with ϕ S,i := ϑ S,i S 0,i and ϕ L,i := ϑ L,i L 0,i . This equation has an analytic solution provided in 17 which will be used for ML. We consider the following parametrization: log(δ S ) = θ δ S , log(ϕ

S,i ) = θ ϕ S + b ϕ S ,i , log(ϕ L,i ) = θ ϕ L + b ϕ L ,i and log(δ Ab,i ) = θ δ Ab + b δ Ab ,i .
The true parameter values are presented in table 3. According to [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF], δ L was non-identifiable based on the available data and only a lower bound has been derived for it via profiled likelihood. So, to make fair comparisons between our approach and maximum likelihood, we do not estimate it. Regarding population parameters, we are particularly interested in the behavior of estimation methods for θ δ S and θ ϕ S . Indeed, a parameter sensitivity analysis shows the symmetric role of θ δ S and θ ϕ S on the , 2) absence of subscript when θ δ S is estimated. Results from our method are in bold.

ODE solution (see [START_REF] Balelli | A model for maintenance and reactivation of the immune response after vaccination against ebola virus[END_REF]). Thus, they are likely to face practical identifiability problems. To investigate this effect, we estimate the parameters when θ * δ S 1) is known (the corresponding estimators will be denoted with the subscript δ S ), or 2) has to be estimated as well.

Well-specified case

We generate n i = 11 longitudinal observations on the interval [0, T ] = [0, 364] with measurement noise of standard deviation σ * = 100. For each subject i, the initial condition has been generated according to A * 0,i ∼ N (A 0 , σ 2 A0 ) with A 0 = 500 and σ A0 = 260 to reflect the dispersion observed in data presented in [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]. We present the estimation results in table 4 -left side.

Our method gives an improved estimation with a dramatically reduced variance for θ * δ S comparing to ML, as well as an improved estimate for the {b * i } i∈ 1, n in all cases. We assume that is due to the committed estimation error for θ * which causes model error during {b * i } i∈ 1, n estimation, not accounted for by ML. This in turn explains why variance Ψ * is better estimated with our approach. In this mixed-effect context, this cause of model error is systematically present and claims for the use of estimation methods taking it into account when subject specific parameters are critical for the practitioner.

Misspecified case in presence of model error at the subject level

The data are generated with a stochastic perturbed version of ODE ( 17):

dA i = ϕ S,i e -δ S t + ϕ L,i e -δ L t -δ Ab,i A i dt + αdB t (18)
where B t is a Wiener process and α = 10 its diffusion coefficient. The value for α has been chosen big enough to produce significantly perturbed trajectories but small enough to ensure that ODE ( 17) is still a relevant approximation for estimation purpose. The results are presented in table 4 -right side. Our method outperforms the ML for θ * δ S as well as for {b * i } i∈[1, n] estimation and their variances. However, we acknowledge that this last simulation setting is challenging even for our approach with confidence coverage around 90% for most of parameters, below the theoretical rate of 95%.

Real data analysis

We use the presented estimation approach to address the same problem as [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]. This real data example is similar to the synthetic scenario performed in Section 3.3. In brief, we use data from a phase I trial in East Africa evaluating the effect of an heterologous anti-Ebola vaccine strategy in which Ad26.ZEBOV first and then MVA-BN-Filo with a delay of 28 days between the two doses. We consider a population of n=28 individuals, with in average 5 measurements per subject. In order to ensure a fair comparison, we adopt a Bayesian framework for θ = (θ δ S , θ ϕ S , θ ϕ L , θ δ Ab ) and used the same prior distribution as in the original paper:

π(θ) ∼ N         -1 0 0 -4.1     ,     25 0 0 0 0 100 0 0 0 0 100 0 0 0 0 1         .
We set our mesh-size to get 200 discretization points for each subject and we use U = 10 i.e., a value lower than in the simulated data case because of the presence of model error. We also proceed to the log-transformation of the data to stabilize the measurement noise variance. Using the transformation A i (t) := log 10 A i (t) in Equation ( 17), this drives us to use the following nonlinear model:

˙ A i (t) = 1 ln(10) ϕ S,i e -δ S t + ϕ L,i e -δ L t 10 -Ai(t) - δ Ab,i ln(10) . ( 19
)
We choose A θ,bi (t, x, z i (t)) = 1 ln(10) ϕ S,i e -δ S t + ϕ L,i e -δ L t 10 -x

x and r θ,bi (t, z i (t)) = -δ Ab,i ln(10) for the pseudo-linear formulation of the model. In Table 5, we compare our estimations with those presented in [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF] obtained using the NIMROD software [START_REF] Prague | Nimrod: A program for inference via a normal approximation of the posterior in models with random effects based on ordinary differential equations[END_REF] produce estimations with overlapping confidence intervals for θ supporting the previous published results in term of antibodies concentrations dynamics over time. Still, significant differences appear for (Ψ ϕ S , Ψ ϕ L , Ψ δ Ab ) estimation with lower dispersion of random effects in the optimal control approach. This is explained by the fact that a part of the variability is now carried out by subject-specific perturbations u i, θ,bi( θ) i∈ 1, n

. Figure 2 allows to visually check that the individual fits are comparable between the two approaches. Finally, our method can be used to assess the model adequacy via the temporal evolution analysis of u i, θ,bi( θ) i∈ 1, n estimated as byproducts of our method. In Section 2.1, we have also indicated that perturbations u i,θ,bi,x0,i can be computed for an arbitrary set (θ, b i , x 0,i ). In particular, we estimate u i, θ P , bi P ,y0,i i∈ 1, n , the committed error corresponding to ( θ P , b i P

), the population and subject specific estimators obtained in [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]. In Figure 3, we plot both perturbation sets. Our method leads to residual perturbations of smaller magnitudes and narrower confidence intervals. This means that our approach produces an estimation which minimizes the committed model error for each subject comparing to a method based only on a data fitting criterion as in [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]. Moreover, by reducing the size of the confidence interval for estimated perturbations, we conclude to a mean perturbation among the population which is statistically different from zero at the beginning of observation interval. This may indicates presence of model misspecification.

Conclusion

In this paper, we propose an estimation method that addresses problems encountered by classical approaches in NLME-ODE models. We identify the following shortcomings for exact methods such as likelihood-based inference: their difficulty in the presence of model misspecification, their need to estimate initial conditions as regular random effects, and their dramatic performance degradation in the presence of poorly identifiable parameters. We propose here a method based on control theory that accounts for the presence of potential model uncertainty at the subject level and that can be easily profiled on initial conditions. Simulations with and without model error illustrate the advantages Parameter estimation in NLME-ODEs: an optimal control approach Dashed lines: fitted ODE solutions from [START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]. Solid line: optimal trajectories X θ, b i obtained with optimal control approach. Shaded area are the 95% confidence intervals.

of regularization techniques for estimating poorly identifiable parameters, subject-specific parameters, and their variances in NLME-ODEs. In addition, bypassing the estimation of initial conditions represents a clear advantage for partially observed systems comparing to likelihood based approaches, as emphasized in the simulations. Still, this benefit in term of estimation accuracy comes with a computational price. On a server (see https://plafrim-users.gitlabpages.inria.fr/doc/ for more server details) with the parallelization package Snow in R language, it takes approximately 10-15 minutes to obtain an estimation for the twodimensional linear model, 30 minutes for the insulin model and 3-4 hour for the antibody concentration evolution one, whereas it was a matter of minutes for the other approaches. Nevertheless, the use of compiled languages and proper Parameter estimation in NLME-ODEs: an optimal control approach parallelization could reduce the computation time. Moreover, we have willingly separated the formal definition of the optimal control problem required by our method and the numerical procedure used to solve it, in case it may exists better suited approaches for this specific control problem. Right now, our current strategy allows us to profile on initial conditions, therefore looking for another numerical procedure is beyond the scope of this paper. Finally, the qualitative assessment of model misspecification exposed in section 4 can be made more rigorous. In a one subject setting, the estimation of a perturbation term at the derivative level via non-parametric procedures to test model error presence has been already explored [START_REF] Hooker | Goodness of fit in nonlinear dynamics: misspecified rates or misspecified states?[END_REF]; [START_REF] Engelhardt | A bayesian approach to estimating hidden variables as well as missing and wrong molecular interactions Parameter estimation in NLME-ODEs: an optimal control approach in ordinary differential equation-based mathematical models[END_REF]). Comparing to statistical methods solely based on data fitting criteria, they generally produce more sensitive statistical tests and can explore misspecification presence even for unobserved state-variables. Our control based approach can extend such tests to a population framework, while avoiding issues due to hyperparameter selection required for non-parametric statistical methods which can appear for a growing number of subjects. For example, to stay in a Bayesian setting, we can specify a prior distribution for the controls and then compare it with the obtained posterior once the inference is made. This would lead to a semi-parametric inference problem for which optimal control based approach has already been proven useful (see [START_REF] Clairon | A regularization method for the parameter estimation problem in ordinary differential equations via discrete optimal control theory[END_REF]). This is a subject for further work. Supplementary information. A supplementary file containing the appendixes and proof of theorem 1 is available alongside this article.

Fig. 1

 1 Fig.1Left: Examples of solutions of (12) and corresponding observations. Right: Solution of (12) and a realization of (13) for the same parameter values.

  and we aim to estimate θ = (θ S G , θ S I , θ m ), linked to the original model via the parametrization: log(S G ) = θ S G , log(S I ) = θ S I and log(m i ) = θ m + b i where b i ∼ N (0, Ψ). The true population parameter values are θ * = (-3.89, -7.09, -1.81) and Ψ * = 0.26 2 . The true subject-specific initial conditions x

Fig. 2

 2 Fig.2Examples of fitted trajectories for both methods for four different random subjects. Dashed lines: fitted ODE solutions from[START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF]. Solid line: optimal trajectories X θ, b i obtained with optimal control approach. Shaded area are the 95% confidence intervals.

  and 95% CI for OCA estimation

Fig. 3

 3 Fig. 3 1) Up: Estimated residual controls for each subject, 2) bottom: mean optimal control and 95% confidence interval for the optimal controls a) left: u i, θ P , b i P ,y 0,i obtained from parameter estimation in Pasin et al (2019), b) right: u i, θ,b i ( θ) obtained from our estimation.

Table 1

 1 Results of estimation for model (

			Well-specified					Misspecified
			MSE	Bias	V e	V	CR	MSE bi	MSE	Bias	V e	V	CR	MSE bi
	θ1	θ M L x 0 θ M L x 0,2 θ M L	0.01 0.01 0.04	0.01 0.01 -0.04	0.01 0.01 0.04	0.01 0.01 0.01	0.95 0.94 0.86		0.01 0.01 0.05	4e-4 -3e-4 0.02	0.01 0.01 0.05	0.01 1e-4 0.01	0.91 0.89 0.81
		θ	5e-3	8e-3	8e-3 1e-2 0.97		0.01	-8e-3	7e-3	0.05	0.97
	θ2	θ M L x 0 θ M L x 0,2 θ M L	4e-5 6e-5 4e-3	1e-3 1e-3 -0.01	4e-5 6e-5 3e-3	4e-5 8e-5 1e-4	0.95 0.94 0.80		1e-4 1e-4 4e-3	-1e-3 -1e-3 -2e-3	1e-4 2e-4 4e-3	1e-4 0.01 2e-4	0.83 0.82 0.63
		θ	5e-5	2e-3	4e-5 4e-5 0.93		1e-4	2e-5	1e-4 1e-4 0.92
	Ψ	θ M L x 0 θ M L x 0,2 θ M L	0.01 0.02 0.05	-0.03 -0.03 0.17	0.01 0.01 0.02	7e-3 7e-3 0.02	1 1 1	0.01 0.01 0.10	0.01 0.01 0.09	-0.003 -0.005 0.21	0.01 0.01 0.04	0.01 0.01 0.03	1 1 1	0.01 0.01 0.12
		θ	0.01	-0.01	0.01	0.01	0.92	0.01	0.02	-0.02	0.02	0.01	0.90	0.01

Table 4

 4 Results of estimation for model(17). The different subscripts stand for the following estimation scenarios: 1) δ S when θ δ S is set to θ *

			Well-specified					Misspecified
			MSE	Bias	V e	V	CR	MSE bi	MSE	Bias	V e	V	CR	MSE bi
		θ M L δ S	known						known
	θδ S	θML θδ S	2.13 known	0.78	1.51	70.64 0.92		3.88 known	1.48	1.68	4.10	0.80
		θ	0.62	-0.34	0.50	0.66	0.92		0.93	-0.40	0.77	0.62	0.90
		θ M L δ S	4e-4	0.01	3e-4	3e-4	0.94		1e-3	0.02	1e-3	5e-4	0.91
	θϕ S	θML θδ S	0.01 2e-3	-0.05 -0.05	7e-3 2e-4	0.40 1e-3	0.92 0.94		0.02 7e-4 -0.02 -0.10	0.01 3e-4 1e-3 0.92 0.02 0.88
		θ	2e-3	1e-3	2e-3	2e-3	0.93		4e-3 -6e-3 3e-3 0.01	0.90
		θ M L δ S	3e-3	0.02	3e-3	2e-3	0.95		5e-3	0.03	4e-3	3e-3	0.93
	θϕ L	θML θδ S	4e-3 7e-4	0.03 -0.01	4e-3 5e-4	3e-3 3e-3	0.90 0.95		9e-3 2e-3 -0.02 0.05	7e-3 3e-3 2e-3 0.97 4e-3 0.90
		θ	3e-3	-3e-3	3e-3 2e-3	0.91		6e-3 -8e-3 6e-3 7e-3 0.90
		θ M L δ S	7e-4	-0.02	5e-4	3e-4	0.93		2e-3	-0.03	1e-3	1e-3	0.92
	θδ Ab	θML θδ S	2e-3 2e-4	-0.02 0.01	1e-3 1e-4	4e-4 3e-4	0.88 0.95		4e-3 3e-4 2e-3 -0.04	3e-3 3e-4 3e-4 0.96 7e-4 0.88
		θ	4e-4	0.01	3e-4	2e-4	0.90		3e-4 8e-3	3e-4 2e-3 0.89
		θ M L δ S	0.04	-1e-3	0.04	0.07	1	0.15	0.05	0.03	0.05	0.08	1	0.17
	Ψϕ S	θML θδ S	0.11 0.02	0.01 8e-3	0.11 0.02	0.05 0.01	1 0.94	0.17 0.06	0.13 0.02	0.01 2e-3	0.13 0.02	0.25 0.02	1 0.94 0.11 0.21
		θ	0.02	-0.03	0.02	0.02	0.94	0.07	0.02	-0.05	0.02	0.03	0.92 0.08
		θ M L δ S	0.03	0.04	0.02	0.04	1	0.30	0.05	0.03	0.05	0.06	1	0.73
	Ψϕ L	θML θδ S	0.03 0.02	0.05 -0.1	0.02 5e-3	0.04 8e-3	1 0.93	0.60 0.07	0.03 0.02	0.05 -0.10	0.02 0.01	0.07 0.02	1 0.91 0.10 0.74
		θ	0.03	-0.06	0.02	0.01	0.92	0.08	0.03	-0.06	0.02	0.03	0.87 0.12
		θ M L δ S	0.11	0.18	0.08	0.02	1	0.10	0.33	0.41	0.17	0.05	1	0.56
	Ψδ Ab	θML θδ S	0.20 0.10	0.29 -0.30	0.11 0.01	0.02 0.01	1 0.95	0.50 0.03	0.30 0.10	0.34 -0.16	0.19 0.08	0.05 0.06	1 0.91 0.04 0.69
		θ	0.11	-0.27	0.04	0.04	0.95	0.04	0.15	-0.29	0.06	0.10	0.88 0.06
									δ S

  ). Both methods Estimation presented in[START_REF] Pasin | Dynamics of the humoral immune response to a prime-boost ebola vaccine: quantification and sources of variation[END_REF] and via our approach.

		Estimations from Pasin et al (2019) Optimal Control appraoch
	Parameter	Mean IC95%	Mean IC95%
	θ δ S	-0.57 [-1.02, -0.02]	-0.18 [-0.58, 0.22]
	θ ϕ S	7.92 [7.52, 8.30]	7.45 [6.85, 7.96]
	θ ϕ L	2.78 [2.62, 3.01]	2.58 [2.15, 3.01]
	θ δ Ab	-3.54 [-3.62, -3.45]	-3.48 [-3.95, -3.01]
	Ψ ϕ S	0.92 [0.83, 1.01]	0.64 [0.60, 0.70]
	Ψ ϕ L	0.85 [0.78, 0.92]	0.70 [0.55, 0.90]
	Ψ δ Ab	0.30 [0.24, 0.36]	0.25 [0.19, 0.31]
	Table 5

Parameter estimation in NLME-ODEs: an optimal control approach

Acknowledgments. Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/). This manuscript was developed under WP4 of EBOVAC3.

Funding: This work has received funding from the Innovative Medicines Initiative 2 Joint Undertaking under projects EBOVAC1 and EBOVAC3 (respectively grant agreement No 115854 and No 800176). The IMI2 Joint Undertaking receives support from the European Union's Horizon 2020 research and innovation programme and the European Federation of Pharmaceutical Industries and Association.

Code availability: Our estimation method is implemented in R and a code reproducing the examples of Section 3 is available on a GitHub repository located here.