Raphaël Kuate
email: raphael.kuate@inria.fr

Marc Lavielle
email: marc.lavielle@inria.fr

Eric Blaudez
email: eric.blaudez@lixoft.net

Kaelig Chatel
email: kaelig.chatel@lixoft.net

Jerome Marquet
email: jerome.marquet@lixoft.net

Jean-François Si Abdallah
email: jean-francois.si_abdallah@lixoft.net

A delay differential equation solver for Monolix & MlxPlore

Keywords: delay differential equations, solver, algorithm, Runge-Kutta formulas, modeling, pharmacokinetics, pharmacodynamics, Monolix, MlxPlore équations différentielles avec retard, solveur, algorithmique, formules de Runge-Kutta, modélisation, pharmacocinétique, pharmacodynamique, Monolix, MlxPlore

We describe the implementation of a Delay differential Equation (DDE) solver in Monolix, a platform for population modeling of longitudinal data, and MlxPlore, a tool for the exploration of complex models. We use explicit Runge-Kutta schemes. Several examples for MlxPlore and Monolix are proposed.

Introduction

Monolix (MOdèles NOn LInéaires à effets miXtes) [2, [START_REF] Lavielle | Estimation of population pharmacokinetic parameters of saquinavir in HIV patients with the MONOLIX software[END_REF][START_REF] Phylindal | The use of the SAEM algorithm in MONOLIX software for estimation of population pharmacokinetic-pharmacodynamic-viral dynamics parameters of maraviroc in asymptomatic HIV subjects[END_REF][START_REF] Ruiz | Ceftazidime dosage regimen recommendations in burn patients based on a monolix population pharmacokinetic study[END_REF] is a platform of reference for model-based drug development. It combines the most advanced algorithms with unique ease of use. Pharmacometricians of preclinical and clinical groups can rely on Monolix for population analysis and to model pharmacokinetic/pharmacodynamic (PK/PD) and other complex biochemical and physiological processes [START_REF] Makowski | Using SAEM to estimate parameters of models of response to applied fertilizer[END_REF][START_REF] Lavielle | An improved SAEM algorithm for maximum likelihood estimation in mixtures of non linear mixed effects models[END_REF][START_REF] Radojka | Implementation and evaluation of the SAEM algorithm for longitudinal ordered categorical data with an illustration in pharmacokinetics-pharmacodynamics[END_REF]. Monolix is an easy, fast and powerful tool for parameter estimation in non-linear mixed effect models, model diagnosis and assessment, and advanced graphical representation.

MlxPlore [1] is a graphical and interactive software for the exploration and visualization of complex pharmacometric models. MlxPlore also includes the ability to study the statistical variability of the models, and to model and study complex administrations designs. MlxPlore is based on the Mlxtran model description language popularized by the Monolix software, therefore benefiting from the remarkable flexibility and power of Mlxtran to easily encode complex pharmacometric, PKPD and statistical models. MlxPlore does not require Monolix, although they make for a powerful combination, enabling to use the same, human-readable model description, to finely explore the properties of the model on the one hand, and on the other hand use the same model for advanced parameter estimation in the context of population analysis and mixed effect statistics.

Modeling complex biological phenomenons with ordinary differential equations (ODEs) sometimes involves a large number of variables and parameters, which makes the analysis of these models cumbersome. Modeling such systems with delay differential equations (DDEs) helps to describe the important dynamics with fewer variables and parameters [START_REF] Srividhya | The effects of time delays in a phosphorylation-dephosphorylation pathway[END_REF]. Another advantage of using a DDE model over an ODE model is that the parameters in the DDE model usually have a direct biological interpretation [START_REF] Baker | Mathematical modelling of the interleukin-2 T-cell system: A comparative study of approaches based on ordinary and delay differential equations[END_REF]. Thus, providing Monolix and MlxPlore with a DDE solver is of a major importance.

Observing that the Taylor expansion of higher order of the solution of a DDE system needs derivatives of the delayed terms, we augment the DDE system being solved with the derivatives of the delayed terms. Then we apply classical explicit Runge-Kutta schemes on the new augmented system. We have implemented two variants of timestep: an adaptive timestep and a fixed timestep of which we propose its computation with respect to the CFL condition. We also propose an algorithm for the generation of the fixed time grid such that no interpolation of delayed terms be computed.

The document is organized as follows. In the second section we describe the implementations of the algorithms. Two variants of timestep are used: a fixed timestep and an adaptive timestep. In the fixed timestep case, a version uses the delayed solution, computed simultaneously during the simulation and a second version uses interpolated delayed solution. The numerical schemes implemented in the fixed timestep case are explicit Runge-Kutta formulas of order 2, 3 and 4. In the adaptive timestep version, the numerical schemes coded are Runge-Kutta embedded pairs [START_REF] Prince | High order embedded Runge-Kutta formulae[END_REF][START_REF] Dormand | Runge-Kutta triples[END_REF][START_REF] Bogacki | A 3(2) pair of Runge-Kutta formulas[END_REF][START_REF] Sofroniou | Construction of explicit Runge-Kutta pairs with stiffness detection[END_REF] of type 2(1), 3(2) and 4(3), using of course interpolated delayed solution. In the third section, we give some examples of numerical tests on DDEs, solved using MlxPlore. In the fourth section, we present examples of estimation of parameters of DDE models, computed with Monolix.

Algorithms

The algorithms implemented use the well-known Runge-Kutta numerical schemes [START_REF] Christopher | Computing stability regions-Runge-Kutta methods for delay differential equations[END_REF][START_REF] Butcher | A history of Runge-Kutta methods[END_REF][START_REF] Butcher | Runge-Kutta methods: Some historical notes[END_REF][START_REF] Calvo | Stepsize selection for tolerance proportionality in explicit Runge-Kutta codes[END_REF][START_REF] Vanden Berghe | Exponentially-fitted explicit Runge-Kutta methods[END_REF][START_REF] Sofroniou | Construction of explicit Runge-Kutta pairs with stiffness detection[END_REF][START_REF] Ch | Runge-Kutta interpolants for high precision computations[END_REF][START_REF] Ch | On modified Runge-Kutta trees and methods[END_REF].

Runge-Kutta schemes

ODEs Let Y be a time-dependent vector of real unknowns and let Ẏ be the time derivative of Y . Let the vector F (t, Y ) be a function of Y and of the time t, and let Y 0 be the initial condition at the starting time t 0 . The usual form of a differential equation of order 1 which can be found in PK/PD or biological [START_REF] Gabrielsson | Methodology for pharmacokinetic/pharmacodynamic data analysis[END_REF][START_REF] Pérez-Urizar | Pharmacokinetic-pharmacodynamic modeling: Why?[END_REF][START_REF] Rosenbaum | Basic pharmacokinetics and pharmacodynamics: An integrated textbook and computer simulations[END_REF][START_REF] Smith | Mathematical model of a three-stage innate immune response to a pneumococcal lung infection[END_REF][START_REF] Krzyzanski | Lifespan based indirect response models[END_REF] models writes:

Ẏ (t) = F (t, Y ), t ∈ [t 0 , T ], Y (t 0 ) = Y 0 .
(2.1.1)

The numerical approximation of the solution of equation (2.1.1) using standard explicit Runge-Kutta numerical scheme is done as follows. Let t n ∈ [t 0 , T ] be a value of the time grid, and denote Y n the value of Y at time t n . Let s, p ≥ 1 be two integers and let h n = t n+1 -t n be the timestep. An explicit Runge-Kutta scheme of s stages and of order p writes:

           Y n+1 = Y n + h n s i=1 b i K i , K i = F t n + c i h n , Y n + i-1 j=1 a i,j K j , i = 2, ..., s, K 1 = F (t n , Y n ), (2.1.2)
where the coefficients a ij , b i , c i ∈ R are obtained by solving a system of equations [START_REF] John | Coefficients for the study of Runge-Kutta integration processes[END_REF][START_REF] John | On the attainable order of Runge-Kutta methods[END_REF] which results of an identification of (2.1.2) with the Taylor expansion of order p of Y , with respect to (2.1.1).

DDEs Many developments have been made [START_REF] Christopher | Developing a delay differential equation solver[END_REF][START_REF] Baker | Issues in the numerical solution of evolutionary delay differential equations[END_REF][START_REF] Luzyanina | Computation, continuation and bifurcation analysis of periodic solutions of delay differential equations[END_REF][START_REF] Shampine | Solving DDEs in matlab[END_REF][START_REF] Shampine | Solving ODEs with MATLAB[END_REF][START_REF] Bellen | Numerical methods for delay differential equations[END_REF][START_REF] Castro | Quantization-based integration methods for delay-differential equations[END_REF][START_REF] Bauer | Solving delay differential equations in S-ADAPT by method of steps[END_REF] on numerical methods for delay differential equations. We use standard explicit Runge-Kutta formulas, adapted to DDEs. For sake of simplicity, we describe the derivation of an explicit Runge-Kutta numerical scheme for DDEs, in the case s = p = 2.

Let N be a positive integer and τ ∈ R N + be the vector of delays and denote Y τi = Y (t -τ i ). Let Y 0 be a function of the time t, and of τ i , 1 ≤ i ≤ N . The usual form of a delayed differential equation of order 1 which can be found in PK/PD or biological [START_REF] Baker | Mathematical modelling of the interleukin-2 T-cell system: A comparative study of approaches based on ordinary and delay differential equations[END_REF][START_REF] Alexander Y Mitrophanov | Control of streptococcus pyogenes virulence: Modeling of the CovR/S signal transduction system[END_REF][START_REF] Srividhya | The effects of time delays in a phosphorylation-dephosphorylation pathway[END_REF][START_REF] Koch | Multi-response model for rheumatoid arthritis based on delay differential equations in collagen-induced arthritic mice treated with an anti-GM-CSF antibody[END_REF] models writes:

Ẏ (t) = F (t, Y, Y τ1 , ..., Y τ N ), t ∈ [t 0 , T ], Y (t) = Y 0 (t, τ ), t ≤ t 0 . (2.1.3)
The Taylor expansion of order 2 of Y gives:

Y n+1 = Y n + h n Ẏn + 1 2 h 2 n Ÿn + O(h 3 n ). (2.1.4)
Denote F Y the matrix of the derivatives of F with respect to Y . The time derivative of Y of order 2, expressed in terms of F with respect to (2.1.3), writes:

Ÿ = Ḟ + F Y F + N i=1 F Yτ i Ẏτi .
(2.1.5)

If one does not neglect Ẏτi , it appears from (2.1.5) that one needs Ẏτi in order to identify the coefficients of an explicit Runge-Kutta formula for the numerical approximation of (2.1.3). Thus, system (2.1.3) is augmented as follows:

   Ẏ (t) = F (t, Y, Y τ1 , ..., Y τ N ), t ∈ [t 0 , T ], Ẏτi (t) = F (t -τ i , Y τi , Y τ1+τi , ..., Y τ N +τi ), t ∈ [t 0 , T ], 1 ≤ i ≤ N. Y (t) = Y 0 (t, τ ), t ≤ t 0 .
(2.1.6)

Inria Define Z := {Y, Y τ1 , ..., Y τ N }, set Z 0 (t) = Y 0 (t, τ ), t ≤ t 0 , and define
F (t, Z) :=              F (t, Y, Y τ1 , ..., Y τ N ) F (t -τ 1 , Y τ1 , Y τ1+τ1 , ..., Y τ N +τ1 ) F (t -τ 2 , Y τ2 , Y τ1+τ2 , ..., Y τ N +τ2 )
. . .

F (t -τ N , Y τ N , Y τ1+τ N , ..., Y τ N +τ N )             
.

One has:

Ż(t) = F (t, Z), t ∈ [t 0 , T ], Z(t) = Z 0 (t), t ≤ t 0 . (2.1.7) System (2.1.7
) is used instead of (2.1.3) for the computations with standard Runge-Kutta formulas. Thus, an explicit Runge-Kutta scheme of s stages and of order p for the numerical approximation of (2.1.3) writes:

           Z n+1 = Z n + h n s i=1 b i K i , K i = F t n + c i h n , Z n + i-1 j=1 a i,j K j , i = 2, ..., s, K 1 = F (t n , Z n ), (2.1.8)
where the coefficients a ij , b i , c i ∈ R are obtained by solving a system of equations which results of an identification of (2.1.8) with the Taylor expansion of order p of Z, with respect to (2.1.7).

Solver features

Fixed timestep

The solver has a feature which uses a fixed timestep, with two sub-features. The first sub-feature uses cubic Hermite interpolation for the computation of delayed terms, since the solution and its time first order derivative are computed simultaneously. The second sub-features does not use interpolation. Let Ω be the discretization of the time interval.

Ω = {t k , t k ∈ [t 0 , T ]} is computed such that: ∀t n ∈ Ω, t n -τ i ≥ t 0 =⇒ t n -τ i ∈ Ω, 1 ≤ i ≤ N. (2.2.1)
The table T Ω containing elements of Ω is returned with a table of numbers of the entries of T Ω which contains t n -τ i , ∀t n ∈ Ω, 1 ≤ i ≤ N. Thus for any Y n = Y (t n ) computed, one knows the delayed terms Y (t n -τ i ) already computed, since statement (2.2.1) is always satisfied.

In the second sub-feature where the time grid contains delayed terms, the computation of T Ω is done using Algorithm 1, for a given a timestep ∆ t .

Computation of the timestep In numerical simulation of time-dependent equations, numerical instabilities can occur, particularly for explicit numerical schemes [START_REF] Rogers | Stability and convergence of approximation schemes[END_REF][START_REF] Gopal | A review of recent developments in solving ODEs[END_REF][START_REF] Sanz-Serna | Regions of stability, equivalence theorems and the Courant-Friedrichs-Lewy condition[END_REF][START_REF] Van Der Houwen | Explicit Runge-Kutta (-Nyström) methods with reduced phase errors for computing oscillating solutions[END_REF][START_REF] Christopher | Computing stability regions-Runge-Kutta methods for delay differential equations[END_REF][START_REF] Petzold | Numerical solution of highly oscillatory ordinary differential equations[END_REF][START_REF] Kalogiratou | Runge-Kutta type methods with special properties for the numerical integration of ordinary differential equations[END_REF]. Since the timestep ∆ t is fixed, one must compute ∆ t with respect to the Courant-Friedrichs-Lewy [START_REF] Courant | Über die partiellen differenzengleichungen der mathematischen physik[END_REF][START_REF] Courant | On the partial difference equations of mathematical physics[END_REF] condition. Let #z be the size of Z and let I #z be the identity matrix of the size of Z. The first order time derivative of Z is approached from (2.1.7) with the following:

Ż ≈ F Z Z,
of which an explicit numerical discretization writes:

Z n+1 ≈ I #z + ∆ tn {F Z } n Z n .

RR n°8489

Algorithm 1 Time grid computation.

1: T Ω inserts t 0 , t 1 ← T , inserts t 1 into T Ω in descending order. 2: while the biggest stepsize between t 1 and t 0 is greater than ∆ t do 3:

t j ← t 1 , t k ← t 1 4:
while t j ≥ t 0 do 5:

for all i, 1 ≤ i ≤ N do 6:

t i ← t j -τ i 7: if t i ≥ t 0 then 8:
put a flag on t i 9:

inserts t i into T Ω in descending order. in T Ω , remove the flag on t j 14:

t k ← t j 15:
end while 16:

t 1 ← t 1 -∆ t 17:
if t 1 ≥ t 0 then 18:

inserts t 1 into T Ω in descending order.

19:

end if 20: end while 21: rearrange T Ω in ascending order

The timestep ∆ t is then computed such that (2.2.2) be satisfied.

I #z + ∆ t F Z Z = CF L < 1.
(2.2.2)

Let z i , 1 ≤ i ≤ #z, be the i-th entry of Z. The timestep ∆ t is obtained from (2.2.2) with

z i = 1 √ #z , 1 ≤ i ≤ #z, i.e., Z = 1.
Let Z i , 1 ≤ i ≤ #z, be the vector of the size #z containing z i at the i-th entry and zeros elsewhere. Since F Z may be time dependent, let F k Z be its value for t = t k . The i-th column of F k Z is obtained by evaluating the following expression

F i,k Z = F (t k , Z + αZ i ) -F (t k , Z) α , 1 ≤ i ≤ #z, O < α ≪ 1.
For each t k chosen, one obtains a value ∆ t k of ∆ t from (2.2.2). The time interval [t 0 , T ] is uniformly meshed and

t k = t 0 + k #nk (T -t 0 ), 0 ≤ k ≤ #nk, #nk < 20.
The timestep ∆ t is taken as the minimum of the ∆ t k 's. However, one checks if F Z is not time dependent by comparing F (t 0 , Z) with F (t 1 , Z), and the computations are stopped if F (t 0 , Z) ≡ F (t 1 , Z).

Inria

Coefficients of the Runge-Kutta formulas In the fixed timestep feature of the solver, the numerical schemes implemented are explicit Runge-Kutta formulas of order 2, 3 and 4:

           Z n+1 = Z n + h n s i=1 b i K i , K i = F t n + c i h n , Z n + i-1 j=1 a i,j K j , i = 2, ..., s, K 1 = F (t n , Z n ).
The coefficients a i,j , b i and c i can be stored in a table called Butcher tableau defined in table 1.

c 2 a 2,1 . . . . . . . . . c s-1 a s-1,1 a s-1,2 . . . c s a s,1 a s,2 . . . a s,s-1 b 1 b 2 . . . b s-1 b s
Table 1: The Butcher tableau of the Runge-Kutta coefficients.

The coefficients used are given in the Butcher tableaux, table 2, and can be found in [START_REF] John | Coefficients for the study of Runge-Kutta integration processes[END_REF][START_REF] Butcher | Runge-Kutta methods: Some historical notes[END_REF][START_REF] Butcher | A history of Runge-Kutta methods[END_REF].

Order 2:

1 2 1 2 1 2 1 2
Order 3:

1 2 1 2 1 2 -1 2 1 6 2 3 1 6
Order 4:

1 2 1 2 1 2 0 1 2 1 0 0 1 1 6 1 3 1 3 1 6
Table 2: Some Butcher tableaux of Runge-Kutta explicit scheme of order 2, 3 and 4.

Adaptive timestep

The model considered is the augmented system (2.1.7). In an adaptive timestep method, the timestep is computed such that the method of order p be more accurate than an embedded method of order p ′ < p. Generally, one choses p ′ = p -1. In a step n, the local error

E n = Z n -Z ′
n between the approximation Z of order p and the approximation Z ′ of order p ′ is used to validate or not the step and to compute the timestep for the next step.

Z n+1 = Z n + h n s i=1 b i K i , Z ′ n+1 = Z ′ n + h n s i=1 b ′ i K i .
The coefficients of the Runge-Kutta formula for these two approximations are computed such that the coefficients a i,j , and c i (2.1.8) of the two approximations be the same. Some embedded pairs also use the FSAL (first same as the last) strategy i.e., when a step is accepted, the last coefficient K s becomes the first coefficient K 1 of the next step, which reduces the method to a s-1 stage method. Many strategies [START_REF] Calvo | Stepsize selection for tolerance proportionality in explicit Runge-Kutta codes[END_REF][START_REF] Shampine | Solving DDEs in matlab[END_REF][START_REF] Sofroniou | Construction of explicit Runge-Kutta pairs with stiffness detection[END_REF][START_REF] Hindmarsh | Sundials: Suite of nonlinear and differential/algebraic equation solvers[END_REF][START_REF] Ch | Runge-Kutta interpolants for high precision computations[END_REF][START_REF] Ch | On modified Runge-Kutta trees and methods[END_REF] have been developed for the computation of adaptive timesteps. The strategy used is a standard formula described in [START_REF] Calvo | Stepsize selection for tolerance proportionality in explicit Runge-Kutta codes[END_REF][START_REF] Ch | Runge-Kutta interpolants for high precision computations[END_REF]. Given a tolerance T OL the algorithm computes the next timestep h n+1 with the following formula:

h n+1 = θ h n T OL E n 1 p , θ ∈ (0, 1).
The next stepsize h n+1 is accepted if T OL ≥ E n . If T OL < E n the current step n is rejected and recomputed again with h n ← h n+1 . The first timestep h 1 is computed [START_REF] Calvo | Stepsize selection for tolerance proportionality in explicit Runge-Kutta codes[END_REF] as

h 1 = T OL max F (t 0 , Z 0 ) , 10 -p 1 p .
The Butcher tableau for Runge-Kutta embedded pairs is defined in table 3.

c 2 a 2,1 . . . . . . . . . c s-1 a s-1,1 a s-1,2 . . . c s a s,1 a s,2 . . . a s,s-1 b 1 b 2 . . . b s-1 0 b ′ 1 b ′ 2 . . . b ′ s-1 b ′ s Table 3:
The Butcher tableau of coefficients of the Runge-Kutta embedded pairs.

Runge-Kutta coefficients In the adaptive timestep feature of the solver, the numerical schemes implemented are Runge-Kutta embedded FSAL pairs of M. Sofroniou and G. Spaletta [START_REF] Sofroniou | Construction of explicit Runge-Kutta pairs with stiffness detection[END_REF], of type 2(1), 3(2) and 4(3). The Butcher tableaux are given in tables 4 and 5.

Inria

Type 2(1):

1 1 1 1 2 1 2 1 2 1 2 0 1 - 1 6 1 6
Type 3(2): 3 Simulation of DDEs using MlxPlore

1 2 1 2 1 -1 2 1 1 6 2 3

Example 3.2

The model solved is the Interleukin-2: model studied by Baker et al. in [START_REF] Baker | Mathematical modelling of the interleukin-2 T-cell system: A comparative study of approaches based on ordinary and delay differential equations[END_REF] and defined by:

                                   İ2 (t) = -α I2 I 2 (t) -n I2 b T I2 I 2 (t) I 2 (t)/I * 2 + 1 T A (t), ṪA (t) = ρb T I2 I 2 (t -τ D ) I 2 (t -τ D )/I * 2 + 1 T A (t -τ D ) -b T I2 I 2 (t -τ S ) I 2 (t -τ S )/I * 2 + 1 T A (t -τ S ) -α AR T A (t), ṪD (t) = b T I2 I 2 (t -τ S ) I 2 (t -τ S )/I * 2 + 1 T A (t -τ S ) -b T I2 I 2 (t -τ D ) I 2 (t -τ D )/I * 2 + 1 T A (t -τ D ), ṪR (t) = α AR T A (t) -α R T R (t), t ∈ [0, 100]. (3.2.1)
The parameters of model (3.2.1) are

α I2 = 0, α R = 1.5 • 10 -2 , α AR = 6.6 • 10 -2 , n I2 = 4755, b T I2 = 1.8 • 10 -11 , I * 2 = 6 • 10 10 , τ D = 6.6, τ S = 10, ρ = 2. (3.2.2)
The model (3.2.1) is compared in [START_REF] Baker | Mathematical modelling of the interleukin-2 T-cell system: A comparative study of approaches based on ordinary and delay differential equations[END_REF] with an ODE model of the same phenomenon, defined as follows:

                   İ2 (t) = -α I2 I 2 (t) -n I2 b T I2 I 2 (t) I 2 (t)/I * 2 + 1 T A (t), ṪA (t) = ρb D T D (t) -b T I2 I 2 (t) I 2 (t)/I * 2 + 1 T A (t) -α AR T A (t), ṪD (t) = b T I2 I 2 (t) I 2 (t)/I * 2 + 1 T A (t) -b D T D (t), ṪR (t) = α AR T A (t) -α R T R (t), t ∈ [0, 100].
(3.2.3)

RR n°8489

The parameters of model (3.2.3) are

α I2 = 0, α R = 3.5 • 10 -2 , α AR = 0.02, n I2 = 4755, b T I2 = 1.0 • 10 -11 , I * 2 = 6 • 10 10 , b D = 1/8, ρ = 2.
The history of models (3.2.1) and (3.2.3) is given by

I 2 (t) = 2 • 10 10 , T A (t) = 3.8 • 10 5 , T D (t) = 0, T R (t) = 1.2 • 10 5 , t ≤ 0.
The MlxPlore script of the dde model of this example is given in the section below Example 3.3

The model solved is the PKPD model for unperturbed and perturbed arthritis development of Koch et al. [START_REF] Koch | Multi-response model for rheumatoid arthritis based on delay differential equations in collagen-induced arthritic mice treated with an anti-GM-CSF antibody[END_REF] defined by:

       Ġ(t) = k 3 -σ 1 exp (-σ 2 c(t)) + σ 3 c(t)G(t) - k 1 k 2 1 -exp (-k 2 t) G(t), İ(t) = k 4 G(t) -k 4 G(t -τ ), Ḋ(t) = k 4 G(t -τ ) -k 5 D(t), t ∈ [0, 35].
The history is given by

G(t) = a exp(bs), I(t) = I 0 , D(t) = 0, t ≤ 0.
The parameters are a = 1, b = 0.5, k 1 = 0.183, k 2 = 0.092, k 3 = 5, k 4 = 0.064, k 5 = 0.016, τ = 11.2.

The parameter c(t) is given by a PK data dose of a standard linear 2-compartment i.v. model. The equation used in [START_REF] Koch | Multi-response model for rheumatoid arthritis based on delay differential equations in collagen-induced arthritic mice treated with an anti-GM-CSF antibody[END_REF] is c(t) = dose A iv exp(-αt) + B iv exp(-βt) .

However, c(t) can be fitted with MlxPlore, knowing the volumes V 1 , V 2 of the compartments, the parameters α, β, CL, the dose (amount in Mlxtran) and the administration times of the PK model. The MlxPlore script of this example is given in section the below MlxPlore script

<MODEL> [PREDICTION]

Table 5 :

 5 The Butcher tableau of Runge-Kutta FSAL pairs of M. Sofroniou and G. Spaletta[START_REF] Sofroniou | Construction of explicit Runge-Kutta pairs with stiffness detection[END_REF], of type 4(3).

Example 3. 1 Figure 1 :

 11 Figure 1: The SEIR epidemic model of Genik & van den Driessche, example 4.4.1 of [40].

Figure 2 :

 2 Figure 2: The Interleukin-2: model studied by Baker et al. in [4]

Figure 3 :

 3 Figure 3: The PKPD model for unperturbed and perturbed arthritis development of Koch et al. [24], Experiment A with the dosing schedules: 10 mg/kg on day 1, 8, 15.

  j ← the next entry of T Ω with a flag, after t k in descending order

	10:	end if
	11:	end for
	12:
	13:

t

Table 4 :

 4 The Butcher tableaux of Runge-Kutta FSAL pairs of M. Sofroniou and G. Spaletta[START_REF] Sofroniou | Construction of explicit Runge-Kutta pairs with stiffness detection[END_REF], of types 2(1) and 3(2).

	1
	6

RR n°8489

rue Honoré d'Estienne d'Orves Bâtiment Alan Turing Campus de l'École Polytechnique 91120 Palaiseau Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

 [START_REF] Shampine | Solving ODEs with MATLAB[END_REF]defined by:

The history is given by [2] Monolix A software for the analysis of nonlinear mixed effects models.

http://www.lixoft.eu/products/monolix/documentation/.

Inria