Agostinho Agra
email: aagra@ua.pt

Marielle Christiansen
email: marielle.christiansen@iot.ntnu.no

Rosa Figueiredo
email: rosa.figueiredo@ua.pt

Lars Magnus Hvattum
email: lars.m.hvattum@iot.ntnu.no

Michael Poss
email: mjposs@gmail.com

Cristina Requejo
email: crequejo@ua.pt

Layered Formulation for the Robust Vehicle Routing Problem with Time Windows

Keywords: vehicle routing problem, robust programming, time windows, maritime transportation, layered formulation

This paper studies the vehicle routing problem with time windows where travel times are uncertain and belong to a predetermined polytope. The objective of the problem is to find a set of routes that services all nodes of the graph and that are feasible for all values of the travel times in the uncertainty polytope. The problem is motivated by maritime transportation where delays are frequent and must be taken into account. We present an extended formulation for the vehicle routing problem with time windows that allows us to apply the classical (static) robust programming approach to the problem. The formulation is based on a layered representation of the graph, which enables to track the position of each arc in its route. We test our formulation on a test bed composed of maritime transportation instances.

Introduction

In this paper, we study the vehicle routing problem with time windows in the uncertain context. Given a graph with a special node called depot and a set of vehicles, the vehicle routing problem aims at prescribing routes for the vehicles starting at and returning to the depot in such a way that each remaining node of the graph is visited by exactly one vehicle. The problem has numerous applications in transportation, distribution and logistics, see [START_REF] Golden | The vehicle routing problem : latest advances and new challenges[END_REF]. In this work, we are more particularly interested by an application that arises in maritime transportation [8].

Among the many versions of the vehicle routing problem that have been studied in the literature, we consider the problem where time windows are given for each node of the network, yielding the vehicle routing problem with time windows (V RP T W ). Hence, each node must be serviced during specific time intervals (or time windows) and traveling along the arcs consumes time. Most authors consider in addition that each vehicle has also a capacity that can not be exceeded along its route. In this work however, we consider the problem version without the capacity constraint, often called m-T SP T W in the literature. This assumption is motivated by our application in maritime transportation where each ship carries only one cargo at the time, from the loading port to the unloading port. Since it is straightforward to extend our model and solution method to the problem with the capacity constraint, we keep the notation V RP T W in what follows.

Exact solution methods for the V RP T W have been studied extensively and many integer programming formulations have been proposed for the problem, see the reviews [START_REF] Desaulniers | The Vehicle Routing Problem with Time Windows: State-of-the-Art Exact Solution Methods[END_REF][START_REF] Kallehauge | Formulations and exact algorithms for the vehicle routing problem with time windows[END_REF]. Among them, Bard [START_REF] Bard | A branch-and-cut procedure for the vehicle routing problem with time windows[END_REF] studies the resource inequalities formulation for the problem. The formulation from [START_REF] Bard | A branch-and-cut procedure for the vehicle routing problem with time windows[END_REF] contains two sets of variables: arc variables state which arcs belong to the solution and node variables indicate at what time vehicles arrive at each node. Kallehauge et al. [START_REF] Kallehauge | Path inequalities for the vehicle routing problem with time windows[END_REF] extend to the V RP T W the path inequalities formulation proposed in [START_REF] Ascheuer | A polyhedral study of the asymmetric traveling salesman problem with time windows[END_REF] for the asymmetric traveling salesman with time windows. In [START_REF] Kallehauge | Path inequalities for the vehicle routing problem with time windows[END_REF], routes or paths that cannot satisfy the time windows are cut-off by path inequalities. Formulations based on path variables are also very popular for all versions of vehicle routing problems, including the V RP T W , see [START_REF] Desaulniers | The Vehicle Routing Problem with Time Windows: State-of-the-Art Exact Solution Methods[END_REF][START_REF] Kallehauge | Column Generation, chap. Vehicle routing problem with time window[END_REF], among others. These formulations contain a little number of constraints but a very large number of variables so that efficient branch-and-price algorithms are required.

In this work, we study the V RP T W in the uncertain context where travel times are not known with precision and belong to an uncertainty polytope T . Hence, our approach falls into the framework of robust programming. Robust programming stems from the original work of [START_REF] Soyster | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF] and has witnessed a continuous attention it the last ten years. We refer the interested reader to the book [START_REF] Ben-Tal | Robust optimization[END_REF].

The robust vehicle routing problem with time windows and uncertain travel times (T -V RP T W ) has been mentioned already in [START_REF] Sungur | A robust optimization approach for the capacitated vehicle routing[END_REF]. However, their modeling assumption leads to all travel times taking their maximum values, which is an over-conservative model. In fact, [START_REF] Sungur | A robust optimization approach for the capacitated vehicle routing[END_REF] mainly focus on the robust capacitated vehicle routing problem, see also [START_REF] Ordónez | Robust vehicle routing[END_REF]. The literature on the stochastic version of the V RP T W is also scant, among which [START_REF] Chardy | Handling uncertainties in vehicle routing problems through data preprocessing[END_REF][START_REF] Russell | Vehicle routing with soft time windows and erlang travel times[END_REF] which considers stochastic travel times. Hence, the present work is the first general approach to the robust vehicle routing problem with time windows and uncertain travel times. Travel times belong here to a demand uncertainty polytope, which makes the problem harder to solve, yet tractable, than its deterministic counterpart. The retribution of the addition in complexity is that our model is more flexible than the one from [START_REF] Sungur | A robust optimization approach for the capacitated vehicle routing[END_REF] and leads to less conservative robust solutions.

Our objective in this paper is to make use of the classical framework of (static) robust programming. In that framework, a vector is feasible for the problem if and only if it is feasible for all values of the travel time in the uncertainty polytope T . Hence, the formulation based on resource inequalities cannot be used in that context because it makes no sense to choose arrival times that are independent of the travel times (see Example 1 in the next section). In fact, as explained below, none of the existing formulations for the V RP T W can be used with the classical approach for static robust programming.

The classical approach for static robust programming under polyhedral uncertainty relies on dualizing the constraints that contain uncertain parameters [?]. Hence, the approach requires that the uncertainty parameters appear explicitly in the constraints of the problem. This is not the case of the formulations based, respectively, on path inequalities and path variables. In each of these formulations, the uncertain parameter appear implicitly in the paths that define the inequalities and the variables, respectively. Moreover, other formulations proposed in the literature present one (or both) of the two problems: some of their variables are somehow related to the values taken by travel times (nonstatic robust programming) or the travel times do not appear explicitly in the constraints of the formulations.

For this reason, we propose an extended formulation for the V RP T W that is suitable for robust programming. This formulation is based on two ideas: i) considering a layered representation of the routes followed by the vehicles (see [START_REF] Gouveia | Using variable redefinition for computing lower bounds for minimum spanning and steiner trees with hop constraints[END_REF], among others), and ii) rewriting the time windows in an extended form [START_REF] Chardy | Handling uncertainties in vehicle routing problems through data preprocessing[END_REF]. We apply the dualization technique to this formulation to provide a formulation for the T -V RP T W . We assess our formulation numerically on instances modeling a problem that arises in maritime transportation.

This paper is structured as follows. Next section presents our extended formulation for the V RP T W , denoted by (LF ). The dualization technique is recalled in Section 3.1 and applied to formulation (LF ) in Section 3.2. Our numerical experiments are described in Section 4 while some concluding remarks are given in Section 5.

Extended Formulation for the V RP T W

We are given a directed graph G = (N, A), a set of vehicles K, a cost function c : A × K → R + , and a time function t : A × K → R + for traveling along the arcs of G. The graph contains special depot nodes o and d connected to all other nodes of G, and we denote by N * the set of nodes that are not depots, N * := N \{o, d}. We are given time windows [a i , b i ] with a i , b i ∈ R, for each i ∈ N * . The V RP T W consists of defining routes for the vehicles in K such that the union of all routes passes exactly once by each i ∈ N * . When |K| = 1, the problem contains a unique vehicle and reduces to the Asymmetric Traveling Salesman with Time Windows, see [START_REF] Ascheuer | A polyhedral study of the asymmetric traveling salesman problem with time windows[END_REF].

We first recall the classical resource inequality formulation for the problem and show through an example why it cannot be extended to the (static) robust context. The formulation uses a set of binary flow variables x k ij which indicates whether vehicle k travels from node i ∈ N to node j ∈ N , and a set of continuous variables y k i indicating the arrival time of vehicle k at node i ∈ N . Then, the satisfaction of the time windows is expressed by the following set of constraints

x k ij (y k i + t k ij -y k j ) ≤ 0, (i, j) ∈ A, k ∈ K, (1)
a i ≤ y i ≤ b i , i ∈ N * , (2)
where ( 1) is linearized using classical "big-M" techniques. Extending the formulation to the robust context would require that inequalities (1) and ( 2) be satisfied for all values of t in the uncertainty polytope T . We show in Example 1 that this does not work.

Example 1. Consider N * = {1, 2, 3, 4} and time windows [START_REF] Ascheuer | A polyhedral study of the asymmetric traveling salesman problem with time windows[END_REF][START_REF] Bard | A branch-and-cut procedure for the vehicle routing problem with time windows[END_REF], [START_REF] Ben-Tal | Robust optimization[END_REF][START_REF] Ben-Tal | Robust solutions of uncertain linear programs[END_REF] and [START_REF] Bertsimas | Theory and applications of robust optimization[END_REF][START_REF] Bertsimas | The price of robustness[END_REF] for, respectively, nodes 1,2 and 3. Suppose that the uncertainty polytope T is defined as follows:

T := {(t 12 , t 23 ) = (1 -λ)(3, 2) + λ(2, 3), 0 ≤ λ ≤ 1}. Now, consider a path p := o → 1 → 2 → 3 → d.
It is easy to see that the time windows are feasible along p for all t ∈ T . Consider now the binary vector x such that x o1 = x 12 = x 23 = x 3d = 1. Constraints (1) become y 2 ≥ y 1 + t 12 and y 3 ≥ y 2 + t 23 . Because each of these constraints must be satisfied for all t ∈ T , they become y 2 ≥ y 1 + 3 and y 3 ≥ y 2 + 3. Because the smallest feasible value for y 1 is 1, the smallest feasible value for y 2 is 4 and it is impossible to find a value for y 3 that satisfies (2) for node 3.

The aim of this section is to provide a formulation for the V RP T W that is easily adaptable to (static) robust programming. Hence, the formulation satisfies two properties: all variables are related to the routes taken by the vehicles (to avoid situations as in Example 1), and travel times appear explicitly in the constraints. The formulation is based on the rewriting of the time windows constraints (1) and (2) as performed in [START_REF] Chardy | Handling uncertainties in vehicle routing problems through data preprocessing[END_REF]. Consider path p = i 0 → . . . → i n for vehicle k and a binary vector x ∈ {0, 1} |A| that describes p, that is, x ij = 1 for each (i, j) ∈ p and x ij = 0 otherwise. The authors of [START_REF] Chardy | Handling uncertainties in vehicle routing problems through data preprocessing[END_REF] show that (1) and ( 2) are satisfied along p if and only if the constraints

a i l 1 + l=l1,...,l2-1 t k i l i l+1 ≤ b i l 2 , 0 ≤ l 1 < l 2 ≤ n (3)
are satisfied. To take advantage of constraints (3), we construct a layered graph which keeps track explicitly of the position of each arc along its path. Layered graphs have been used for many network design problems, starting from Gouveia [START_REF] Gouveia | Using variable redefinition for computing lower bounds for minimum spanning and steiner trees with hop constraints[END_REF], and have already been applied to the V RP (see [START_REF] Godinho | Combined route capacity and route length models for unit demand vehicle routing problems[END_REF], among others). However, to the best of our knowledge, layered graphs in the sense proposed by [START_REF] Gouveia | Using variable redefinition for computing lower bounds for minimum spanning and steiner trees with hop constraints[END_REF] have not yet been applied to the V RP T W . The main idea of the formulation below is to model the flow problem associated to each vehicle with a directed graph composed of L = |N | layers as illustrated in Figure 1. Namely, from the original graph G = (N, A), we create a directed layered graph G = (N , A) for each vehicle, where N := N 1 ∪ . . . ∪ N L with N 1 := {o}, N L := {d} and N l := N \{o}, l = 2, . . . , L -1. Let i l be the copy of i ∈ N in the l-th layer of graph G. Then, the arc sets are defined by {1, . . . , L-1}}∪{(d, d, l), l ∈ {2, . . . , L-1}}, see Figure 1. Hence, (i, j, l) denotes an arc between i l ∈ N l and j l+1 ∈ N l+1 . Note that each path between o and d in the layered graph G is composed of exactly L -1 arcs, that corresponds to a path of less than or equal to L -1 arcs in G.

A := {(i, j, l) | (i, j) ∈ A, i l ∈ N l , j l+1 ∈ N l+1 , l ∈ 2 3 o 1 d o 1 2 3 d 1 2 3 d 1 2 3 d d N 1 N 2 N 3 N 4 N 5 (a) (b)
Then, we introduce the set of additional binary flow-position variables z kl ij for each k ∈ K, (i, j, l) ∈ A defined as follows: z kl ij = 1 whenever vehicle k ∈ K services i ∈ N exactly in position l of its path from the artificial origin node o to the artificial destination node d and just before servicing node j ∈ N . The extended arc-flow model for the deterministic case, denoted by (LF ), follows. min

k∈K (i,j)∈A c k ij x k ij (4) s.t.
k∈K j∈N :(i,j)∈A

x k ij = 1, i ∈ N * , (5)
j:(j,i,l-1)∈A z kl-1 ji - j:(i,j,l)∈A z kl ij =    -1 if (i = o) 1 if (i = d and l = L) 0 else , 1 ≤ l ≤ L, i ∈ N l , k ∈ K, (6)
l:(i,j,l)∈A

z kl ij = x k ij , (i, j) ∈ A, k ∈ K, (7)
(i,j)∈A:(i,j,l1)∈A

a i z kl1 ij + l=l1,...,l2-1 (i,j)∈A:(i,j,l)∈A t k ij z kl ij ≤ (i,j)∈A:(i,j,l2)∈A b j z kl2 ij , 1 ≤ l 1 < l 2 < L, k ∈ K, (8)
x k ij ∈ {0, 1}, (i, j) ∈ A, k ∈ K, (9)
z kl ij ∈ {0, 1}, (i, j, l) ∈ A, k ∈ K. ( 10
)
The objective function (4) minimizes the cost of operating the set of vehicles. Constraints [START_REF] Bertsimas | Theory and applications of robust optimization[END_REF] ensure that all i ∈ N * are serviced exactly once. Equations ( 6) are flow balance constraints in the directed layered graph. Variables x and z are linked by constraints [START_REF] Chardy | Handling uncertainties in vehicle routing problems through data preprocessing[END_REF]. Finally, constraints (8) adapt constraints (3) when path p is not yet decided and depends on variables z.

3 Robust Formulation

Dualization Approach

In this work, we consider that the travel times are uncertain and belong to a polytope T . This makes the problem a robust program, a class of optimization problems that has witnessed a tremendous attention in the recent years. Conducting an exhaustive literature review of robust programming is beyond the scope of this paper and we redirect the interested reader to [START_REF] Ben-Tal | Robust optimization[END_REF] and [START_REF] Bertsimas | Theory and applications of robust optimization[END_REF]. We recall below the well-known dualization technique for linear robust programs under polyhedral uncertainty, introduced by [START_REF] Ben-Tal | Robust solutions of uncertain linear programs[END_REF]. Consider the following linear program in {0, 1}-variables

min c T x (P ) s.t. Bx ≤ b, (11)
T x ≤ d, (12)
x ∈ {0, 1} n ,

with c ∈ R n , b ∈ R r , d ∈ R s , B ∈ R rn , and T ∈ R sn .
Suppose that the problem is subject to uncertainty in the sense that matrix T belongs to a polytope T ⊂ R sn . The robust counterpart of (P ) is

min c T x (T -P ) s.t. Bx ≤ b, T x ≤ d, T ∈ T , (13)
x ∈ {0, 1} n , where the s linear constraints in [START_REF] Gouveia | Using variable redefinition for computing lower bounds for minimum spanning and steiner trees with hop constraints[END_REF] must now be satisfied for each value of T ∈ T . Hence, the finite set of constraints [START_REF] Gouveia | Using variable redefinition for computing lower bounds for minimum spanning and steiner trees with hop constraints[END_REF] has been replaced by the infinite set of constraints [START_REF] Kallehauge | Formulations and exact algorithms for the vehicle routing problem with time windows[END_REF].

The method explained next works in two steps. The first step amounts to realize [3, Section 1.2.1] that x satisfies constrains [START_REF] Kallehauge | Formulations and exact algorithms for the vehicle routing problem with time windows[END_REF] if and only if it satisfies the following constraints

T i x ≤ d i , T i ∈ T i , i = 1, . . . , s, (14)
where T i ⊂ R n is the projection of T into the space corresponding to the coefficients of the i-th row of ( 13), for each i = 1, . . . , s. Said differently, the whole set of constraints T x ≤ d is satisfied for each value of the uncertain matrix T ∈ T if and only if each constraint T i x ≤ d i is satisfied for each value of the uncertain vector T i ∈ T i .

For the second step, we need to describe more precisely the uncertainty polytope and its projections, which we suppose non-empty. Let the projections of T be defined as T i = {T i ∈ R n : A i T i ≤ a i , T i ≥ 0} where matrices A i , i = 1, . . . , s and vectors a i , i = 1, . . . , s have appropriate dimensions and, of course, depend on the definition of the uncertainty polytope T . Then, the (infinite) constraint set associated to each row i of ( 14) can be rewritten as max Ti∈Ti T i x ≤ d i . The optimization problem of the left-hand side is equivalent to

max T i x s.t. A i T i ≤ a i T i ≥ 0. (15)
Because ( 15) is always bounded and feasible, linear programming duality ensures us that its optimal solution is equal to the optimal solution of its dual:

min a i u i s.t. (A i ) T u i ≥ x u i ≥ 0.
Thus, each constraint i of inequality set ( 14) is equivalent to

a i u i ≤ d i (A i ) T u i ≥ x u i ≥ 0, (16)
that is, the infinite number of constraints ( 13) is replaced by a finite number of constraints and variables. Moreover, the numbers of new constraints and variables are equal to the dimensions of A i . The dualization technique described above has been applied to numerous robust linear programs subject to polyhedral uncertainty in the literature, see [START_REF] Bertsimas | The price of robustness[END_REF], among others.

Formulation for T -V RP T W

In this section, we apply the methodology recalled in Section 3.1 to the V RP T W . We consider the budget uncertainty polytope studied by Bertsimas and Sim [START_REF] Bertsimas | The price of robustness[END_REF]: we suppose that each component t k ij of t lies between its mean value t k ij and its peak value t k ij + tk ij and that, for each k ∈ K, at most Γ of them can reach their peak values simultaneously. Formally, this is written as

T Γ = × k∈K T k
Γ where all vectors in T k Γ are of the form t k := t k + δ k tk and δ k satisfies the following constraints:

(i,j)∈A δ k ij ≤ Γ, ( 17
) 0 ≤ δ k ij ≤ 1, (i, j) ∈ A. ( 18
)
The robust version of the problem is obtained by replacing ( 8) with (i,j)∈A:(i,j,l1)∈A

a i z kl1 ij + l=l1,...,l2-1 (i,j)∈A t k ij z kl ij ≤ (i,j)∈A:(i,j,l2)∈A b j z kl2 ij , 1 ≤ l 1 < l 2 < L, k ∈ K, t k ∈ T k Γ . ( 19
)
All variables in (LF ) are first-stage variables since they describe the paths taken by the vehicles. Moreover, travel times only appear in [START_REF] Soyster | Convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF]. Since ( 19) must be satisfied for all t k ∈ T k Γ , it is convenient to rewrite these constraints as (i,j)∈A:(i,j,l1)∈A

a i z kl1 ij + max t k ∈T k Γ (i,j)∈A t k ij l=l1,...,l2-1 z kl ij ≤ (i,j)∈A:(i,j,l2)∈A b j z kl2 ij , 1 ≤ l 1 < l 2 < L, k ∈ K. ( 20
)
Let us introduce a dual variable v kl1l2 and u kl1l2 ij for each constraint in ( 17) and ( 18), respectively, associated to each constraint of [START_REF] Suite | Xpress-Optimizer reference manual[END_REF]. Dualizing the maximization problem in [START_REF] Suite | Xpress-Optimizer reference manual[END_REF] as in ( 15)-( 16) yield the following set of constraints:

(i,j)∈A:(i,j,l1)∈A a i z kl1 ij + (i,j)∈A t k ij l=l1,...,l2-1 z kl ij + Γ v kl1l2 + (i,j)∈A u kl1l2 ij ≤ (i,j)∈A:(i,j,l2)∈A b j z kl2 ij , 1 ≤ l 1 < l 2 < L, k ∈ K, ( 21
)
v kl1l2 + u kl1l2 ij ≥ tk ij l=l1,...,l2-1 z kl ij , (i, j) ∈ A, 1 ≤ l 1 < l 2 < L, k ∈ K, ( 22
)
v kl1l2 ≥ 0, 1 ≤ l 1 < l 2 < L, k ∈ K (23) u kl1l2 ij ≥ 0, 1 ≤ l 1 < l 2 < L, k ∈ K, (i, j) ∈ A, (24)
so that a robust version of (LF ) can be formulated as follows: 5) -( 7), ( 9), ( 10), ( 21) -(24), 4 Computational Experiments

min k∈K (i,j)∈A c ij x k ij (T -LF ) s.t. (

Application to the ship routing and scheduling problem

In this section, we apply our model to the Ship Routing and Scheduling problem with Time Windows and uncertain travel times. The deterministic version of this problem is described in [8] where an integer linear programming formulation is presented. Maritime transportation is the major component in international trade and a key part of many economic sectors. Freight transport volumes in maritime systems have been growing for many years and continues to show an upward trend. A great variety of optimization problems is involved in the improvement of maritime transport systems, which makes maritime transportation a challenging research area. Indeed, there has been an increasing research interest in maritime transportation problems over the last years. We refer to [8] for a discussion of practical and theoretical aspects of maritime transportation operations including the description of prescriptive mathematical models and solution approaches.

Time is a crucial factor in maritime transportation and deliveries must occur during predetermined intervals. This makes the maritime transportation problem a special case of the vehicle routing problem with time windows. However, while it can be acceptable to use estimations for travel and service time for some instances of V RP T W , this is not the case in maritime transportation. Delays are important and schedules must account for them.

We consider a heterogeneous fleet of ships with specific ship characteristics including different cost structures and load capacities. We assume that a ship is loaded to its capacity in a loading port and the cargo is transported directly to its unloading port. Only one cargo is transported at a time and the cargo size is less than or equal to the capacity of the ship. The fleet has sufficient capacity to serve all committed cargoes during the planning horizon. The corresponding loading and unloading ports are known. Time windows are imposed for loading cargoes. Herein, the service time of a cargo is the time from the arrival at its loading port until the time of departure from its unloading port. Ships are charged port and channel tolls when visiting ports and passing channels, and these costs depend on the size of the ship. The remaining variable sailing costs consist mainly of fuel and oil costs, and depend usually on the ship size.

Instance Description

The instances used in this paper have been created with a random instance generator made as realistic as possible. The instance generator is based on a real distance matrix that contains 56 ports from around the world, with actual sailing distances between each pair of ports.

Two non-overlapping subsets of ports are selected as pickup ports and delivery ports respectively, to represent the structure of a company operating within deep sea industrial shipping. Cargoes requests are generated between two ports based on a simple inventory model for the delivery port. Time windows are associated with each cargo based on when the request would be generated and an acceptable time before the delivery should be made.

The instance generator also specifies the possible delay in sailing time for each arc in the network. This delay is calculated based on the time normally required to perform the transportation represented by the arc. Since the planning horizon is long, there is a significant risk of a ship being delayed at some point during its route, but the probability of experiencing longer travel times for all legs would be small. Hence it makes sense to make routes that can handle some delays, with Γ equal to some small number.

In the computational testing we generate five instances for each combination of values for the number of cargoes and number of ships.

Reducing the number of layers

Some vehicle routing problems [START_REF] Godinho | Combined route capacity and route length models for unit demand vehicle routing problems[END_REF] consider that there exists a constraint on the number of nodes that any route can follow. These problems can be modeled through layered formulations such as the one used in this paper with the difference that the number of layers is a small integer L (part of the problem input), instead of being the total number of nodes of the graph.

The V RP T W studied in this paper does not present this additional constraint: any vehicle can visit an arbitrary number of nodes, as long as time windows are satisfied. Hence, we use a layered formulation in this paper in a different purpose: our aim herein is to present a formulation that is suitable for the robust programming dualization approach. The drawback of this approach is that the number of layers used in the formulation is equal to the number of nodes in the graph, which yields very large number of variables and constraints.

We cannot simply reduce the number of layers because this may cut-off the optimal solution to the problem. However, due to the presence of time windows, the vehicles are usually not able to visit all nodes of the graph. With this idea in mind, we apply a pre-processing to the problem which, for each vehicle, computes the longest path that satisfies the time windows. The pre-processing step was solved by an integer programming formulation based on the MTZ-inequalities [START_REF] Miller | Integer programming formulation of traveling salesman problems[END_REF].

Numerical Results

We present in this section computational results for formulations (LF ) and (T -LF ) on instances with 10 and 20 cargoes and a number of ships varying between 1 and 5. The formulations have been coded using the modeling language Xpress Mosel 3.2.3 and solved by Xpress Optimizer 22.01.09 [START_REF] Suite | Xpress-Optimizer reference manual[END_REF]. The experiments were run on computer equipped with a processor Intel Core i5 at 2.53 GHz and 4 GB of RAM memory.

For 10 cargoes, the average longest path for 2, 3, 4, and 5 ships is equal to 6.6, 6.8, 6.1, and 4.96, respectively. For 20 cargoes, these values are 14.8, 12.7, 10.7, and 9.96, respectively. Hence, the reduction increases with the number of ships. This was expected because the instances are generated in such a way that all ships are necessary. Hence, more ships lead to time windows harder to satisfy and smaller feasible paths. For problems with one ship (|K| = 1), the pre-processing has no effect since the ship must visit all nodes of the graph.

Tables 1 and2 present the average solution times for instances with 10 and 20 cargoes, respectively. Column Reduction provides the average solution times necessary to compute the longest paths for each vehicle, while columns below With reduction and Without reduction provide the average solution times necessary to solve the problem to optimality, respectively with and without reducing the number of layers using the pre-processing. A time limit of 3600 seconds has been set. Then, unsolved instances within this limit are written in parentheses and the average solution times are computed without them. We see from Tables 1 and2 that the pre-processing allows to reduce significantly the solution times, especially for instances with 20 nodes. In particular, the 7 instances that could not be solved within the time limit can be solved after the reduction.

Conclusion

In this paper, we present the first robust formulation for the vehicle routing problem with time windows and uncertain travel times. To this aim, we introduce a new layered formulation for the vehicle routing problem with time windows that enables us to apply the classical static robust programming approach. We test this approach on instances that describe a maritime transportation problem. We see that this methodology can solve instances with up to 20 nodes. In order to solve larger instances, future work will address different approaches for the problem that rely on the use of alternative robust programming techniques, such as adjustable robust programming.

Fig. 1 .

 1 Fig. 1. Basic Network for commodity k (a) and its Layered Representation (b).

Table 1 .

 1 Results for instances with 10 cargoes.

			reduction	With reduction	Without reduction
		|K|\Γ			0	1		2	0	1	2
		1	-		-	-		-	0.0746 0.0798 0.0846
		2	0.168		0.319 0.162 0.169 0.299 0.25	0.192
		3	0.615		0.19	0.847 0.769 0.388 0.976 0.907
		4	1.2		0.182 0.426 0.557 0.38	0.8	0.914
		5	1.34		0.104 0.312 0.326 0.466 0.658 0.681
		reduction		With reduction			Without reduction
	|K|\Γ		0	1	2	4	6	0	1	2	4	6
	1	-	-	-	-	-	-	7.16	5.13	8.56	7.58	9.95
	2	6.09	18.8 54.4 140 134 192 22.3	174	269	528	536
	3	8.59	44.4 332 908 1990 592 73.2	265(1) 307(1) 609(1) 741
	4	81.9	13.2 89.1 160 237 253 101	760	828(1) 1161(1) 857(1)
	5	334	13.6 52.8 67 188 153 57.9	218	323	193(1) 1487

Table 2 .

 2 Results for instances with 20 cargoes.

Acknowledgements

This research was carried out with financial support from the DOMinant II project, partly funded by the Research Council of Norway. Agostinho Agra, Rosa Figueiredo, and Cristina Requejo are supported by FEDER founds through COMPETE-Operational Programme Factors of Competitiveness and by Portuguese founds through the CIDMA (University of Aveiro) and FCT, within project PEst-C/MAT/UI4106/2011 with COMPETE number FCOMP-01-0124-FEDER-022690. Michael Poss is supported by FCT under the postdoctoral scholarship SFRH/BPD/76331/2011 and the grant PEst-C/MAT/UI0324/2011.


Sophie Donnet

Adeline Samson

A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models

Keywords: Stochastic differential equations, Pharmacokinetic, Pharmacodynamic, population approach, maximum likelihood estimation, Kalman Filter, EM algorithm, Hermite expansion, Gauss quadrature, Bayesian estimation

This paper is a survey of existing estimation methods for pharmacokinetic/pharmacodynamic (PK/PD) models based on stochastic differential equations (SDEs). Most parametric estimation methods proposed for SDEs require high frequency data and are often poorly suited for PK/PD data which are usually sparse. Moreover, PK/PD experiments generally include not a single individual but a group of subjects, leading to a population estimation approach. This review concentrates on estimation methods which have been applied to PK/PD data, for SDEs observed with and without measurement noise, with a standard or a population approach. Besides, the adopted methodologies highly differ depending on the existence or not of an explicit transition density of the SDE.

Introduction

Pharmacokinetics (PK) aims at describing the relationship between the dose administered and the exposure to the drug, i.e. the total concentration of drug in the body. Pharmacodynamics (PD) quantifies the relationship between the drug exposure and the response to this exposure. PK/PD models are often described by differential systems derived from physiology. In general, the proposed models are deterministic, that is, the observed kinetic/dynamic is driven exclusively by internal deterministic mechanisms. However, real pharmacological processes are always exposed to influences that are not completely understood or not feasible to model explicitly. Ignoring these phenomena in the modeling may affect the estimation of PK/PD parameters and the derived conclusions. Therefore there is an increasing need to extend the deterministic models to models including a stochastic component. A natural extension of deterministic differential equations model is a system of stochastic differential equations (SDEs), where relevant parameters have been modeled as suitable stochastic processes, or stochastic processes have been added to the driving system equations [START_REF] Ditlevsen | Introduction to stochastic models in biology[END_REF].

The first papers encouraging the introduction of random fluctuations in PK/PD were published by D' Argenio and Park (1997) and Ramanathan (1999a,b). The authors underline that PK/PD have contributions from both deterministic and stochastic components: drug concentrations follow determinable trends but the exact concentration at any given time is not completely determined.

For example, [START_REF] Ramanathan | A method for estimating pharmacokinetic risks of concentrationdependent drug interactions from preclinical data[END_REF] proposes a stochastic one-compartment PK model with a variable elimination rate. More sophisticated PK/PD models have then been proposed with multiple compartments, non-linear or time-inhomogeneous absorption or elimination (see for example [START_REF] Ferrante | A stochastic compartmental model with long lasting infusion[END_REF]Tornøe et al., 2004a;Ditlevsen and De Gaetano, 2005b;Ditlevsen et al., 2005;[START_REF] Picchini | Modeling the euglycemic hyperinsulinemic clamp by stochastic differential equations[END_REF].

Parameter estimation for SDE has been highly tackled in the statistical literature, often motivated by financial applications (see [START_REF] Sørensen | Parametric inference for diffusion processes observed at discrete points in time: A survey[END_REF], for a review). However, many suggested solutions require high frequency data and are not suited for PK/PD data where designs are usually sparse. Especially, estimation methods based on approximations of the continuous-time observation likelihood (namely the Girsanov formula), which require a high number of data and a small time step between two successive observations, are not adapted. Moreover, PK/PD data are more and more analyzed through a population approach when data from several subjects are considered simultaneously. This yields to PK/PD models with random parameters. Combining SDE with a population approach is quite appealing but raises inference challenges.

In this paper, we concentrate our review on estimation methods adapted to the particular characteristics of PK/PD data. After a short presentation of some examples of PK/PD SDEs in Section 2, we introduce some preliminary comments on the likelihood functions depending on the considered observation model (with and without measurement noise) in Section 3. Section 4 is about estimation methods for standard PK/PD SDE: when SDEs are directly observed, the reader is introduced to techniques based on (i) exact maximum likelihood estimator when explicit solution is available or (ii) Hermite expansion of the transition density, (iii) approximation of the spectral density if the SDE has no explicit solution. When the SDE is observed with an additive measurement error, methods are based on (iv) Kalman filter and its extended version or (v) Monte Carlo approximation of the likelihood. Section 5 discusses estimation methods within a population approach: we detail methods based on (i) exact maximum likelihood estimator when linear SDE with random effect and no measurement noise are considered, (ii) Gauss-Hermite quadrature to approximate the likelihood, possibly coupled with (iii) Hermite expansion of the transition density, (iv) Bayesian approach, (v) Kalman filter and linearization of the likelihood, (vi) Expectation-Maximisation algorithm. The paper finishes with some discussion (Section 6).

Stochastic PK/PD models

In this section, we present some stochastic compartmental PK/PD models that have been proposed in the literature. This list is far from being exhaustive but aims at presenting typical situations, each of them involving a different level of statistical inference difficulty.

where V is the volume of the compartment. This equation has an explicit solution: C t = D V e -ket . Now, assume that k e is not constant in time but randomly fluctuates around a mean value as k e + ξ t , where ξ t is a Gaussian white noise process. Then ξ t dt can be written as γdB t , where B t is a Brownian motion and γ is a constant parameter. Incorporating this noise into the deterministic model, C t becomes a stochastic process, solution of the following SDE:

dC t = -k e C t dt + γC t dB t , C 0 = D V . (1)
That process -known as geometric Brownian motion-has an explicit expression

C t = D V e -ket exp - γ 2 2 t + γB t .
This stochastic process, which is log-normal, only takes positive values, which is noticeable when modeling concentration. Parameters to be estimated are θ = (k e , V, γ).

A stochastic one compartment PK model with first-order absorption has also been considered

by [START_REF] Ferrante | A stochastic compartmental model with long lasting infusion[END_REF]:

dC t = (k a /V -k e C t )dt + γdB t , C 0 = D V , (2)
where k a is the absorption rate. This process -known as an Ornstein-Uhlenbeck process-has an explicit expression

C t = D V e -ket + k a V k e (1 -e -ket ) + γ t 0 e -ke(t-s) dB s ,
and is Gaussian with explicit mean and variance. Parameters to be estimated are θ = (k e , k a , V, γ).

Multi compartments PK models have also been extended to stochastic versions. For example, a two compartments model is described by the deterministic system

dC1t dt = -k 1 C 1t + k 2 C 2t -k e C 1t dC2t dt = k 1 C 1t -k 2 C 2t (3)
where C 1t , C 2t are both compartments concentrations, and k 1 , k 2 are transfer constants. [START_REF] Cuenod | Parameter estimation and change-point detection from dynamic contrast enhanced mri data using stochastic differential equations[END_REF] propose a stochastic version of (3) with independent Brownian motions B 1t , B 2t on each equation and diffusion coefficients γ 1 , γ 2 :

dC 1t = (-k 1 C 1t + k 2 C 2t -k e C 1t ) dt + γ 1 dB 1t dC 2t = (k 1 C 1t -k 2 C 2t ) dt + γ 2 dB 2t . (4)
The yielding stochastic process has an explicit Gaussian solution. Parameters to be estimated are Ditlevsen and De Gaetano (2005b) add a Brownian motion only on the second compartment C 2t of (3), with a diffusion coefficient γ 2 :

θ = (k 1 , k 2 , k e , γ 1 , γ 2 ).
dC 1t = (-k 1 C 1t + k 2 C 2t -k e C 1t ) dt dC 2t = (k 1 C 1t -k 2 C 2t ) dt + γ 2 dB 2t . (5)
This leads to a hypoelliptic SDE, that is, a stochastic differential system in which only a few equations include a volatility term, the other equations being of ODE type. No explicit solution exists for this SDE. Parameters to be estimated are then θ = (k 1 , k 2 , k e , γ 2 ).

SDEs have also been proposed for glucose kinetics using data from intravenous glucose tolerance test. For example, [START_REF] Tornøe | Grey-box modelling of pharmacokinetic/pharmacodynamic systems[END_REF] develop a stochastic non-linear multi-compartments glucose minimal model. No explicit solution exists for this complex SDE.

Stochastic coupled PK/PD models

Several papers propose a stochastic model for simultaneous analysis of PK/PD data. In that case, either only the PK model, or only the PD model or both models are considered stochastic.

For example, [START_REF] Ferrante | A stochastic formulation of the gompertzian growth model for in vitro bactericidal kinetics: parameter estimation and extinction probability[END_REF] consider a stochastic PD model (Gompertz) of the bacterial count N t under antibiotic effect coupled to a two dimensional deterministic PK model for the antibiotic concentration C t :

dN t = (r -bN t log N t -kC t )N t dt + γdB t C t = Dka V (ka-ke) (e -ket -e -kat ), (6)
with r the intrinsic growth rate, b the growth deceleration rate, k the bacterial effect of the drug, D the dose of antibiotic, V the volume of distribution, k a and k e the absorption and elimination constants, respectively. Log-transformation of this non-linear time-inhomogenous SDE has an explicit solution (see [START_REF] Ferrante | A stochastic formulation of the gompertzian growth model for in vitro bactericidal kinetics: parameter estimation and extinction probability[END_REF]. Parameters to be estimated are θ = (k e , k a , V, r, b, k, γ).

Another example is given by [START_REF] Tornøe | Grey-box modelling of pharmacokinetic/pharmacodynamic systems[END_REF] with a stochastic first-order elimination one-compartment PK model and an indirect response PD model (Emax):

dC t = -k e C t + γ C dB 1t dR t = k in -k out 1 + Ct EC50+Ct R t dt + γ R dB 2t , (7)
where R t is the state variable for the PD response, EC 50 is the drug concentration causing 50% of maximal stimulation, γ C and γ R are two diffusion coefficients. This system has no explicit solution.

Parameters to be estimated are θ = (k e , k in , k out , EC 50 , γ C , γ R ).

More complex PK/PD systems have also been introduced in the literature. [START_REF] Picchini | Modeling the euglycemic hyperinsulinemic clamp by stochastic differential equations[END_REF] propose for instance a two dimensional SDE for the oscillations of glycemia occurring in response to hyperinsulinisation and to constant glucose infusion.

Preliminary comments on estimation

We now introduce a general model and the notations used through out the whole paper. Two observations models are considered, depending if the PK/PD SDE is directly observed or observed with an additive measurement noise. The likelihood functions are presented in these two contexts.

Corresponding estimation methods are detailed in Section 4. Observation models with a PK/PD population approach and corresponding estimation methods are gathered in Section 5.

SDE model

The studied biological process is assumed to be modeled by (X t ) t≥0 a stochastic process in X ∈ R p defined as the solution of the multidimensional SDE:

dX t = µ(X t , t, θ)dt + Γ(X t , θ)dB t X 0 = x 0 (θ) (8)
where θ ∈ Θ ⊂ R d is the parameter vector and Θ is the set of parameters, µ(X t , t, θ) : R p ×R×R d → R p is the drift function, Γ(X t , θ) : R p × R d → M(R p ), a p × p matrix, is the volatility function, (B t ) t≥0 denotes a p-dimensional Brownian motion and x 0 (θ) is the initial condition. If the drift and the volatility functions are independent of time t, then the system is said to be homogenous, otherwise it is inhomogenous. The initial condition may be random but for sake of simplicity, we assume it deterministic in this paper. The aim is to estimate the parameter θ from (maybe noisy and partial) observations of X t .

We assume the solution of SDE (8) to exist, and the law of X t given θ and X s = x s for any s < t, to have a strictly positive density, the transition density, w.r.t. the Lebesgue measure on X :

x → p X (x, t -s|x s , θ) > 0, x ∈ X .
This existence is ensured by regularity properties on the drift and volatility functions (see Lipster and Shiryaev, 2001, for instance). Note that this assumption does not imply the transition density to be explicit, it happens only in few cases.

Observation models and likelihood functions

In the first papers published on the estimation of PK/PD SDEs, processes were supposed to be exactly observed, no observation noise was considered. This assumption, rather unrealistic for PK/PD data, has the great advantage to simplify the likelihood, as presented in section 3.2.1. In section 3.2.2, the more realistic assumption of discrete and noisy observations of PK/PD SDE is introduced, as well as the corresponding likelihood.

Discretely and directly observed diffusions

Assume that (X t ) t≥0 is directly observed at discrete times t 1 , . . . , t n . We denote indistinctly by X j or X tj the observation at time t j and

X 1:n = (X t1 , • • • , X tn ) = (X 1 , • • • , X n ). Set ∆ j = t j -t j-1 .
Let p X denote the density of X 1:n given θ. By Markov property, p X is the product of the n -1 transition densities. Thus the likelihood function L(θ;

X 1:n ) is L(θ; X 1:n ) = p X (X 1:n |θ) = n j=1 p X (X tj , ∆ j |X tj-1 , θ). (9)
When p X is explicit, the likelihood (9) is explicit and exact maximum likelihood estimators (MLE) can be computed (see Section 4.1.1). Difficulties in the parameter estimation arise when the transition density p X is not explicit. In that case, approximate estimators have been proposed

in the literature (see Section 4.1.2).

Discretely observed diffusions with observation noise

Let y 1:n = (y 1 , . . . , y n ) be the noisy observations of the diffusion process (X t ) t≥0 at time t 1 , . . . , t n .

We assume that y j ∈ R q , with q ≤ p, is derived from the following statistical model

y j = h(X tj , θ) + g(X tj , θ) ε j , ε j ∼ i.i.d. N (0, I q ) dX t = µ(X t , t, θ)dt + Γ(X t , θ)dB t X 0 = x 0 (θ) (10)
where h(x, θ) and g(x, θ) are two functions from R p ×R d in R q and ε j is a measurement error random q-vector. For the sake of simplicity, we consider the usual assumption of a Gaussian distribution on ε j with identity variance matrix I q . Standard choices for g are g(x, θ) = σ resulting in the homoscedastic model and g(x, θ) = (a+σh(x, θ)) resulting in the heteroscedastic model. Function h is generally the identity function when X t in one-dimensional and a linear application modeling that only a restricted number or a linear combination of components of (X t ) t≥0 is observed, otherwise.

The likelihood function L(θ; y 1:n ) computed from data y 1:n is the marginal distribution p Y :

L(θ; y 1:n ) = p Y (y 1:n |θ) = X p Y |X (y 1:n |X 1:n , θ)p X (X 1:n |θ)dX 1:n , (11)
where p Y |X is the conditional density of y 1:n given X 1:n . If the transition density is Gaussian and h(X, θ) is linear, L(θ; y 1:n ) has an explicit form (see Section 4.2.1). Otherwise, either the transition density, or the integration over X are not explicit. Then the intractable likelihood L(θ; y 1:n ) has to be approximated. These methods are presented in Section 4.2.2.

Estimation for PK/PD SDE

We first present estimation methods adapted to directly observed PK/PD SDE (section 4.1). Then estimation for SDE observed with noise is detailed in section 4.2.

Discretely but directly observed diffusions

Likelihood ( 9) is explicit for few SDEs and analytic estimators have then been derived, as presented in Section 4.1.1. Approximate estimators suited for general SDEs are reviewed in Section 4.1.2.

Explicit or exact estimators

(C t ) t∈[0,T ] .
As data belong to an infinite dimensional space, the continuous-time observation likelihood has to be computed with the Girsanov formula [START_REF] Lipster | Statistics of random processes I : general theory[END_REF]:

L(θ; (C t ) t∈[0,T ] ) = exp - 1 γ 2 T 0 (k a /V -k e C t )dC t - 1 2 T 0 (k a /V -k e C t ) 2 dt where θ = (k e , k a , V, γ). The MLE based on continuous data (C t ) t∈[0,T ] is defined as θ[0,T ] = arg max θ∈Θ {L(θ; (C t ) t∈[0,T ] )}.
It can be proved that θ[0,T ] is asymptotically unbiased, normally distributed and efficient [START_REF] Rao | Statistical Inference for Diffusion Type Processes[END_REF]. However, PK/PD data are measured at discrete times rather than continuously. Discretisation of θ[0,T ] can be considered but it is generally difficult to study the properties of such estimators analytically, especially when the number of data is small and fixed. As an alternative, [START_REF] Ferrante | A stochastic compartmental model with long lasting infusion[END_REF] consider the likelihood (9) based on discrete-time observations C 1:n . They can compute it exactly because PK model ( 2) has an explicit transition density. More precisely, the process (C t ) is Gaussian and the transition density is

p C (y, t -s|z, θ) = P (C t = y|C s = z, θ) = 1 2πv θ (t -s) exp - (y -m θ (z, t -s)) 2 2v θ (t -s) ,
where

m θ (z, t -s) = ka V ke + (z -ka V ke ) exp(-k e (t -s)) and v θ (t -s) = γ 2 2ke 1 -exp(-2k e (t -s)). Then the exact likelihood function is L(θ; C 1:n ) = n j=1 p C (C j , ∆ j |C j-1 , θ). The MLE based on discrete-time observations is defined as θn = arg max θ∈Θ L(θ; C 1:n ).
He is not explicit but can be obtained from a numerical optimization of the likelihood. The MLE θn is asymptotic normal and efficient [START_REF] Lipster | Statistics of random processes I : general theory[END_REF].

In [START_REF] Ferrante | Parameter estimation in a gompertzian stochastic model for tumor growth[END_REF][START_REF] Ferrante | A stochastic formulation of the gompertzian growth model for in vitro bactericidal kinetics: parameter estimation and extinction probability[END_REF], the same estimation method is used for the PK/PD SDE (6) to model bacterial growth and bacteria-drug interactions. The log transformation X t = log N t of SDE ( 6) has a Gaussian explicit solution. The transition density p X (y, t -s|z) is thus explicit and the exact MLE is derived from the likelihood function ( 9) of the discrete data (X 1:n , C 1:n ).

Estimators based on approximation

Transition density of more complex PK/PD SDEs is unknown and approximation has to be used to estimate parameters. Several approaches have been proposed. The simplest one is the Gaussian Euler-Maruyama approximation [START_REF] Oksendal | Stochastic differential equations: an introduction with applications[END_REF] but this scheme is efficient only if ∆ j = t j -t j-1 is small, which is not usual for PK/PD data. We present here two other strategies suited for PK/PD data, one based on an Hermite expansion of the transition density and the second based on the spectral density. This second approach has been introduced for an hypoelliptic SDE. To our knowledge, this is the only estimation method adapted to hypoelliptic SDE.

Hermite expansion of the transition density. This estimation method has been first used in a PK/PD SDE context by [START_REF] Picchini | Maximum likelihood estimation of a timeinhomogeneous stochastic differential model of glucose dynamics[END_REF]. This approach has been originally proposed by [START_REF] Aït-Sahalia | Maximum likelihood estimation of discretely sampled diffusions: a closedform approximation approach[END_REF] for unidimensional and time-homogeneous equations and then extended to inhomoegeneous equations by [START_REF] Egorov | Maximum likelihood estimation of time-inhomogeneous diffusions[END_REF] and to multi-dimensional equations by [START_REF] Aït-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF]. [START_REF] Picchini | Maximum likelihood estimation of a timeinhomogeneous stochastic differential model of glucose dynamics[END_REF] consider a one-dimensional (q = 1) inhomogeneous SDE of glucose dynamics after the steady state of insulin concentration has been reached. The transition density is not explicit, neither the likelihood. Their estimation method works as follows.

They assume the volatility function Γ(x, θ) to be bounded below by a strictly positive function and introduce the Lamberti transform of X t : (X t ) = Xt du Γ(u,θ) . The Hermite expansion p S X of order S of the unknown transition density p X is defined by

p S X (X j , ∆ j |X j-1 , θ) = 1 ∆ j Γ(X j , θ) f 0 (Z j ) S s=0 m s ( (X j-1 , ∆ j ), θ) H s (Z j ), (12)
where H s is the Hermite polynom of order s for s = 1 . . . S, f 0 is the standard normal density function,

Z j = (Xj )-(Xj-1) √ ∆j and m s ( (X j-1 ), ∆ j , θ) = 1 s! E [H s (Z j ) | (X j-1 )]
is the conditional moment of the Hermite polynom H s . Then the approximated likelihood of order S is

L S (θ, X 1:n ) = n j=1 p S X (X j , ∆ j |X j-1 , θ) (13)
and the estimator of θ is defined as

θS n = arg max θ∈Θ L S (θ, X 1:n ).
Aït-Sahalia ( 2002) and [START_REF] Egorov | Maximum likelihood estimation of time-inhomogeneous diffusions[END_REF] prove that under regularity conditions on the drift and volatility functions, θS n is consistent. This method is computationally efficient.

Monte Carlo approximation of the spectral density. Ditlevsen et al. (2005[START_REF] Ditlevsen | Parameter estimation of feedback gain in a stochastic model of renal hemodynamics: differences between spontaneously hypertensive and sprague-dawley rats[END_REF] propose this method for a hypoelliptic PD SDE modeling tubuloglomerular feedback mechanism in a rat nephron. The SDE system consists in 6 equations with Brownian motion on only one coordinate.

There is no analytic expression for the transition density and so for the likelihood function.

Instead of estimating the parameters by maximizing an approximation of the likelihood function,

pK (u, θ) = 1 K K k=1 pk (u, θ).
The strong law of large numbers ensures convergence of pK (u, θ) to p(u, θ). Then the parameters are estimated by nonlinear least squares as

θK n = arg min θ∈Θ u∈[0,n/(2T )] log p(u) -log pK (u, θ) 2 .
Approximate confidence intervals on θK n can also be obtained by Monte Carlo simulations. No theoretical properties have been proved for this estimator.

Discretely observed diffusions with observation noise

Noisy observations of PK/PD SDEs as defined by model ( 10) have been considered in the literature the restriction of a constant diffusion coefficient Γ(X t , θ) = Γ θ . We first present an exact estimation method proposed for linear SDEs (Section 4.2.1), then two approximate estimation methods for more general SDEs (Section 4.2.2).

Explicit or exact estimators

Linear PK/PD SDEs with homoscedastic measurement noise have been studied by D'Argenio and Park (1997); Tornøe et al. (2004a); [START_REF] Favetto | Parameter estimation for a bidimensional partially observed Ornstein-Uhlenbeck process with biological application[END_REF]; [START_REF] Cuenod | Parameter estimation and change-point detection from dynamic contrast enhanced mri data using stochastic differential equations[END_REF]. Model (10) reduces to:

y j = H θ X j + σε j , ε j ∼ iid N (0, I q ), j = 1, . . . , n dX t = (A θ X t + C θ (t))dt + Γ θ dB t (14) setting h(X t , θ) = H θ X t , with H θ a q × p matrix, g(X t , θ) = σ, µ(X t , t, θ) = A θ X t + C θ (t). The p -vector C θ (t) may depend on time t.
By recursive conditioning, likelihood (11) can be decomposed into

L(θ; y 1:n ) = p(y 1 ; θ) n j=2 p(y j |y 1:j-1 ; θ), (15)
where p(y j |y 1:j-1 ; θ) is the conditional law of y j given y 1:j-1 . In the particular case of model ( 14), y j |y 1:j-1 is Gaussian. Indeed, the solution of the linear SDE ( 14) is Gaussian, as well as y j by linearity of h(x, θ). Let us denote m j|1:j-1 and R j|1:j-1 the conditional expectation and variance of y j given y 1:j-1 respectively. We have:

p(y j |y 1:j-1 ; θ) = exp{-1 2 (y j -t m j|1:j-1 )(R j|1:j-1 ) -1 (y j -m j|1:j-1 )} |2πR j|1:j-1 | .
Then m j|1:j-1 and R j|1:j-1 can be computed using the Kalman filtering procedure. Indeed, by linearity of h(x, θ), we have:

m j|1:j-1 = E(y j |y 1:j-1 ) = H θ E(X j |y 1:j-1 ) R j|1:j-1 = V ar(y j |y 1:j-1 ) = H θ V ar(X j |y 1:j-1 ) t H θ + σ 2 I q (16)
where t H θ is the transposed matrix of H θ . Consequently, m j|1:j-1 and R j|1:j-1 rely on the predictive expectations and variances:

Xj|1:j-1 = E(X j |y 1:j-1 ) and P j|1:j-1 = E((X j -Xj|1:j-1 ) t (X j -Xj|1:j-1 ))

whose computations are details below.

Predictive moments. Noting that (X t ) t≥tj verifies X t = t tj-1 µ(X s , s, θ)ds + Γ θ t tj dBs, for any t ≥ t j-1 yields to the following integral equation for the predictive conditional expectation:

E(X t |y 1:j-1 ) = t tj-1 E(µ(X s , s, θ)|y 1:j-1 )ds. ( 17
)
Consequently, drift µ being linear for model ( 14), x t := E(X t |y 1:j-1 ) verifies x t = t tj-1 µ(x s , s, θ)ds, which can be written in a differential form as

dx t /dt = µ(x t , t, θ) x tj-1 = E(X tj-1 |y 1:j-1 ) = Xj-1|1:j-1 . ( 18
)
Using similar arguments, we also obtain an ODE verified by P t := V ar(X t |y 1:j-1 ) for t ≥ t j-1 .

Let J µ denote the Jacobian of µ w.r.t. x, we get

dP t /dt = J µ (x t , t, θ)P t + P t t J µ (x t , t, θ) + Γ θ t Γ θ P tj-1 = V ar(X tj-1 |y 1:j-1 ) = P j-1|1:j-1 . (19)
Because of the linearity of µ, equations (18-19) have explicit solutions. Moreover, noting that Xj|1:j-1 = x tj and P j|1:j-1 = P tj , we obtain Xj|1:j-1 = e A θ ∆j Xj-1|1:j-1 + Cθ (t j-1 , ∆ j )

P j|1:j-1 = e A θ ∆j P j-1|1:j-1 e t A θ ∆j + Γθ (t j-1 , ∆ j ) (20)
where Xj-1|1:j-1 = E(X j-1 |y 1:j-1 ) and P j-1|1:j-1 = V ar(X j-1 |y 1:j-1 ) are the filtering moments and Cθ (t j-1 , ∆ j ), Γθ (t j-1 , ∆ j ) are given in Tornøe et al. (2004a) and [START_REF] Favetto | Parameter estimation for a bidimensional partially observed Ornstein-Uhlenbeck process with biological application[END_REF].

Filtering moments. A similar reasoning yields the update (filtering) step:

Xj|1:j = Xj|1:j-1 + K j (y j -H θ Xj|1:j-1 ),

P j|1:j = (I -K j H θ )P j|1:j-1 ,
where

K j = P j|1:j-1 t H θ (H θ P j|1:j-1 t H θ + σ 2 ) -1 is the Kalman gain.
Finally, iterating the computations of predictive and filtering moments (the so-called Kalman filtering) leads to the exact likelihood ( 15) and the MLE is then defined as

θn = arg max θ∈Θ L(θ; y 1:n ).
Favetto and Samson (2010) prove that θn , which is the exact MLE, is consistent and asymptotically Gaussian. This estimation method is implemented in the Danish Technical University project CTSM [START_REF] Kristensen | CTSM 2.0, UserÕs Guide[END_REF].

This Remark 1 Model (10) with explicit transition density for X t is a state space model or a hidden Markov model [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Corresponding estimation methods have not yet been all applied to PK/PD data.

Approximate estimators

When either the transition density or likelihood L(θ; y 1:n ) is not explicit, an approximation has to be proposed. A first solution relies on an extension of the Kalman filter and has been introduced

by [START_REF] Tornøe | Grey-box modelling of pharmacokinetic/pharmacodynamic systems[END_REF]. Another approach was proposed by Ditlevsen and De Gaetano (2005b) with an Monte Carlo approximation of the likelihood to estimate the volatility parameters, the drift parameters being estimated using the deterministic version of the model.

Extended Kalman Filter. This approach, based on a generalization of the Kalman filter, has been proposed by [START_REF] Tornøe | Grey-box modelling of pharmacokinetic/pharmacodynamic systems[END_REF] for the general model (10) with the restriction of timehomogeneous drift µ(X t , t, θ) = µ(X t , θ):

y j = h(X tj , θ) + σ ε j , ε j ∼ i.i.d. N (0, I q ) dX t = µ(X t , θ)dt + Γ θ dB t X 0 = x 0 (θ) (21)
Once again, the likelihood can be decomposed into L(θ; y 1:n ) = p(y 1 ; θ) n j=2 p(y j |y 1:j-1 ; θ). However, whereas each conditional distribution p(y j |y 1:j-1 ; θ) is Gaussian for linear model ( 14), it is not true for model ( 21). The idea of the extended Kalman Filter (EKF) is thus to approximate this conditional distribution by a Gaussian distribution:

p(y j |y 1:j-1 , θ) exp{-1 2 (y j -t m j|1:j-1 )(R j|1:j-1 ) -1 (y j -m j|1:j-1 )} |2πR j|1:j-1 |
where, as for linear model, m j|1:j-1 = E(y j |y 1:j-1 ) = E(h(X j , θ)|y 1:j-1 ) and R j|1:j-1 = V ar(y j |y 1:j-1 ) =

V ar(h(X j , θ)|y 1:j-1 ) + σ 2 . But in that case, h being not linear, m j|1:j-1 and R j|1:j-1 are not explicit. Thus, EKF considers a local linearisation of h around E(X j |y 1:j-1 ). More precisely, let J h be the Jacobian function of h with respect to x, we can write h(X j , θ) h(E(X j |y 1:j-1 ), θ) + J h (E(X j |y 1:j-1 ), θ)(X j -E(X j |y 1:j-1 )) and approximate the conditional moments by:

m j|1:j-1 h( Xj|1:j-1 , θ) (22)
R j|1:j-1 J h ( Xj|1:j-1 , θ)P j|1:j-1 t J h ( Xj|1:j-1 , θ) + σ 2
where, as before, Xj|1:j-1 and P j|1:j-1 are the conditional predictive expectation and variance.

The difficulty for the nonlinear model ( 21) arises in the computation of these predictive quantities. Using integral equation ( 17), one would like to deduce a solution for Xj|1:j-1 , as in the linear case. However, µ being not linear, another approximation is used, writting:

E(µ(X s , θ)|y 1:j-1 ) µ(E(X s |y 1:j-1 ), θ).
Then, one can say that an approximation of x t := E(X t |y 1:j-1 ) is a solution of differential equation (18). Similarly, using an approximation of the variance, P t := V ar[X t |y 1:j-1 ] can be approximate by a solution of differential equation ( 19). However, except for linear drift µ, equations (18)(19) have no close solution. A new approximation is performed in EKF, consisting in a linearization of the ODEs using Taylor expansion of µ. In [START_REF] Tornøe | Grey-box modelling of pharmacokinetic/pharmacodynamic systems[END_REF], a first order Taylor expansion is considered, leading to the following approximate predictive quantities Xj|1:j-1 = Xj-1|1:j-1 + µ( Xj-1|1:j-1 , θ)∆ j P j|1:j-1 = P j-1|1:j-1 + ∆ j J µ ( Xj-1|1:j-1 , θ)P j-1|1:j-1 + P j-1|1:j-1 J µ ( Xj-1|1:j-1 , θ) + Γ θ t Γ θ

Remark 2 Delattre and Lavielle (submitted) propose a higher order Taylor approximation of equations equations (18-19) (see Section 5.3.3).

The update step is performed using similar Gaussian approximation and linearization:

Xj|1:j = Xj|1:j-1 + K j (y j -h( Xj|1:j-1 , θ))

P j|1:j = (I -K j J h ( Xj|1:j-1 , θ))P j|1:j-1
where K j is the approximated Kalman gain

K j = P j|1:j-1 t J h ( Xj|1:j-1 , θ)(J h ( Xj|1:j-1 , θ)P j|1:j-1 t J h ( Xj|1:j-1 , θ) + σ 2 ) -1 .
Finally, the EKF supplies an iterative procedure to approximate the likelihood, and maximize it yields an estimator of θ. This method is implemented in the CTSM project. However, to our knowledge, no theoretical results have been proved rigorously for this estimator, which can behave badly when both SDE and observation model are highly nonlinear.

This estimation method has been used for PK/PD models. For example, Tornøe et al. (2004b) use PK/PD SDE (7) on simulated data to illustrate the properties of the estimator. Then, for the modeling of intravenous glucose tolerance test with a glucose minimal model (GMM), they use a nonlinear two-dimensional SDE with two compartments (glucose and insulin). The observation is partial and consists in discrete noisy measure of the glucose rate. Data are analyzed both with the SDE and the deterministic GMM. The diffusion term is estimated significantly different from zero. This indicates that the deterministic GMM does not fully capture the glucose/insulin dynamics.

Another example is given in [START_REF] Kristensen | Using stochastic differential equations for PK/PD model development[END_REF] for PK/PD model development.

Monte Carlo approximation of the likelihood. Ditlevsen and De Gaetano (2005b) consider this method with the hypoelliptic two-dimensional PK SDE (5) for the modeling of the uptake of dodecanedoic acid. A diffusion coefficient appears only on one of the two equations. Only the concentration C 1 in the first compartment is observed, with a measurement noise. Ditlevsen and De Gaetano (2005b) consider the heteroscedastic error model

y j = C 1,tj ε j
where ε j are iid log-normal variables with variance var(log ε j ) = σ 2 . The hypoellipticity of PK SDE (5) makes estimation more complex. Indeed, contrary to the twodimensional PK SDE considered by [START_REF] Favetto | Parameter estimation for a bidimensional partially observed Ornstein-Uhlenbeck process with biological application[END_REF], this hypoelliptic linear SDE has no explicit solution or transition densities. Therefore, the likelihood is not available. Ditlevsen and De Gaetano (2005b) thus propose a two-step estimation method. As the deterministic corresponding PK model has an explicit solution, they estimate PK parameters k 1 , k 2 from the deterministic PK model by usual nonlinear least squares estimators. Then, with these parameters fixed, they proceed as a second step to the estimation of γ and σ. However, even when k 1 , k 2 are fixed, the likelihood L(γ, σ; y 1:n ) = n j=1 p(y j |C 1,tj , σ)p(C 1,t1 , . . . , C 1,tn ; γ)dC 1,t1 . . . dC 1,tn has no close form due to the hypoellipticity of (5). Therefore, they approximate the likelihood by Monte-Carlo simulations, as first proposed by [START_REF] Pedersen | A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations[END_REF]. A large number K of trajectories

(C (k)
1,1:n ) k=1,...,K are simulated given γ with the Milstein scheme from SDE (5). Then the likelihood is approximated by

L K (γ, σ; y 1:n ) = 1 K K k=1 n j=1 p(y j |C (k) 1,tj , σ)
and numerically maximized. Estimators of γ, σ are asymptotically equivalent to the MLE. The same approach is used in [START_REF] Picchini | Modeling the euglycemic hyperinsulinemic clamp by stochastic differential equations[END_REF] for modeling the euglycemic hyperinsulinemic clamp. This approach is computationally intensive: K = 500.000 was required to obtain a good approximation of the likelihood. It is also more appropriate to estimate all parameters in a single optimization step.

SDE for population PK/PD models

PK-PD studies often include not a single individual but a group of subjects. In the first papers considering SDE models on a population of subjects, individuals were treated independently without proposing any global model (see for instance [START_REF] Picchini | Modeling the euglycemic hyperinsulinemic clamp by stochastic differential equations[END_REF]Ditlevsen and De Gaetano, 2005b). However, the success of mixed effects models -which allow to discriminate between the inter and intra subjects variability-has encouraged their use for PK/PD SDEs.

We first introduce the population PK/PD SDE in Section 5.1. The estimation methods for population PK/PD SDEs observed without measurement noise are presented in Section 5.2 and for population PK/PD SDEs observed with measurement noise in Section 5.3.

Population PK/PD SDE

Assume that PK/PD data are available from I subjects. Let (X (i) ) t≥0 ∈ R p be the process of interest for subject i, for i = 1, . . . I. We assume that the processes {(X (i) t ) t≥0 , i = 1 . . . I} verify the same SDE but with individual parameters φ i :

dX (i) t = µ(X (i) t , t, φ i , η)dt + Γ(X (i) t , η)dB (i) t X (i) 0 = X 0 (φ i ) φ i ∼ i.i.d. p Φ (•; θ pop ),
where η ∈ R d0 is a parameter common to all the individuals, φ i ∈ R d are the individual random parameters, also called random effects, distributed with a density p Φ depending on a population parameter θ pop ∈ Θ and (B (i) t ) are I independent Brownian motions. A standard choice for p Φ is a Gaussian distribution

φ i ∼ i.i.d. N (µ, Ω)
with mean µ ∈ R d and variance Ω a definite positive matrix of size d. In that case, θ pop = (µ, Ω). However, other distributions can be considered, such as Gamma distributions to ensure the positivity of the parameter, or discrete distributions. The parameters to be estimated are θ = (θ pop , η) from direct discrete observations of X (i) t or noisy discrete observations. Statistical inference for such models is a challenging issue. Some of the proposed methodologies derive from the non-linear mixed effects models literature whereas [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF] extend their method developed for the observation of one diffusion (Section 4.1.2) to the observation of several individual diffusions. Obviously, the difficulty is not the same if the transition density is explicit or not and if the process is observed directly or with observation noise.

Population PK/PD SDE without observation noise

We first assume that the processes (X (i) t ) t≥0 , for i = 1, . . . , I are directly observed on the I individuals at discrete times t ij , for j = 1, . . .

n i . Set ∆ ij = t ij -t ij-1 . Let X (i) 1:ni = (X (i) ti1 , . . . X (i) tin i
) denote the observations of individual i and X = (X (i) 1:ni ) i=1...I , the vector of whole data. By independence of the I individuals, the likelihood is:

L(θ; X) = I i=1 L i (θ; X (i) 1:ni ) = I i=1 Φ p X (X (i) 1:ni |φ i , η)p Φ (φ i |θ pop )dφ i . ( 23
)
where p X (X

(i) 1:ni |φ i , η) = ni j=1 p X (X (i) j , ∆ ij |X (i) j-1 , φ i , η
) similarly as equation ( 9). In a few cases, the likelihood ( 23) can be expressed in a close form, leading to the exact MLE of θ (section 5.2.1). Otherwise, (i) if the transition density p X (x, t -s|x s , φ, γ) has an explicit solution but integrating out the φ i 's in ( 23) is unfeasible; or (ii) if the transition density can not be expressed explicitly, likelihood (23) has no close form and approximate solutions have to be proposed. In case (i), the integral has to be approximated (see Section 5.2.2) whereas in case (ii), the transition density must be approximated and maybe also the integral (see Section 5.2.3).

Exact estimator

Ditlevsen and De Gaetano (2005a) consider the one-compartment linear elimination PK SDE (1) with one random effect:

dC (i) t = -k e,i C (i) t dt + γC (i) t dB (i) t (24) k e,i ∼ N (k e , ω 2 ).
In that case, φ i = k e,i , θ pop = (k e , ω 2 ), η = γ 2 and θ = (k e , ω 2 , γ 2 ). Data are discrete observations

C (i)
tij for i = 1, . . . , I and j = 1, . . . , n i . Using a log-transformation of the data, the likelihood ( 23) can be computed exactly. An estimation method adapted from linear mixed model [START_REF] Pinheiro | Mixed-effect models in S and Splus[END_REF] is applied. Namely, the likelihood is splitted into two parts, that are independent of and dependent of the random effect φ i , respectively. The integral of the dependent part is the integral of a Gaussian density because SDE ( 24) is linear and has a log-normal solution. The explicit expression of the likelihood and of the exact MLE of the parameter θ follows. The MLE is Gaussian, and Ditlevsen and De Gaetano (2005a) give explicit variance of the MLE, leading to explicit confidence intervals.

Estimation for explicit transition density SDE

Excepted for particular SDEs such as (24), even if the transition density p X (x, t -s|x s , φ, γ) has an explicit expression, the difficulty to estimate θ comes from the integration of the individual random effects density. 1 strategy is to approximate the likelihood, that is the integral with respect to the random effects density, and then use a standard maximization procedure. This has been proposed by [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF] with a Gaussian quadrature of the likelihood and by [START_REF] Picchini | Practicle estimation of high dimensional stochastic differential mixed-effects models[END_REF] with a Laplace approximation of the likelihood. The first method has been applied to the one-compartment linear elimination PK SDE with one random effect ( 24) and the second one to the stochastic two-compartment PK model (4) with random effects on k 1 , k 2 , k e . We now detail the two techniques.

Gaussian quadrature. We restrict to the case of a one-dimensional Gaussian individual parameter φ (φ i ∼ i.i.d. N (µ, ω 2 )), we refer the reader to [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF] for the case of any continuous distribution. Assume that φ → p X (X

(i) 1:ni |φ, η) is C 2R (R). Each individual likelihood L i (θ; X (i) 1:ni ) = Φ p X (X (i) 1:ni |φ i , η)p Φ (φ i |θ pop )dφ i is approximated by the Gauss-Hermite quadrature of order R: L R i (θ; X (i) 1:ni ) = R r=1 π r p X (X (i) 1:ni |ω √ 2z r + µ, η)
where z r , r = 1, . . . , R are the zeros of the Hermite polynomial

H R (•) of degree R and π r = 2 R-1 R! R 2 (H R-1 (zr)
) 2 are adequate weights. Neither the zeros z r nor the weights π r depend on the individual. The approximate likelihood is then defined as

L R (θ; X) = I i=1 L R i (θ; X (i) 1:ni ) = I i=1 R r=1 π r p X (X (i) 1:ni |ω √ 2z r + µ, η) (25)
The convergence of L R (θ; X) towards L(θ; X) is ensured when R tends to infinity and the domain of integration is compact. The parameter θ can be estimated by maximizing this approximate

likelihood θ R n = ( η R n , ω R n ) = arg min θ∈Θ - I i=1 log R r=1 π r p X (X (i) 1:ni |ω √ 2z r + µ, η)
The estimator θ R n is obtained by a standard optimization algorithm. No theoretical convergence has been proved for this estimator.

Estimation for unknown transition density SDE

When the transition density has no explicit expression, [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF] propose to approximate it with a Hermite expansion, as detailed in Section 4.1.2. Then this approximate transition density is injected in the approximate marginal distribution (25).

More precisely, the transition likelihood is approximated by p S X (X

(i) j , ∆ ij |X (i) j-1 , φ i , η
) defined by equation ( 12) where S is the order of the Hermite expansion. The conditional likelihood p X (X (i) 1:ni |φ i , η) of subject i can be approximated by:

p S X (X (i) 1:ni |φ i , η) = ni j=2 p S X (X (i) j , ∆ ij |X (i) j-1 , φ i , η) (26)
exactly as equation ( 13). Then, the integration over the individual parameters φ i is obtained using the Gauss-Hermite quadrature (25). Finally, the approximate likelihood is:

L R,S (θ; X) = I i=1 R r=1 π r p S X (X (i) 1:ni |ω √ 2z r + µ, η).
As before, the MLE estimator of θ is obtained applying an optimization routine on the approximate likelihood:

θ R,S n = arg min θ∈Θ -log L R,S (θ; X)
No theoretical result has been proved but [START_REF] Picchini | Practicle estimation of high dimensional stochastic differential mixed-effects models[END_REF] illustrate the efficiency of θ R,S n on several uni-and multi-dimensional examples (note that the multidimensionality comes from both the dimension of φ and the dimension of the diffusion X t ), especially for the stochastic two-compartment PK model (4) with random effects on k 1 , k 2 , k e .

Population SDE models with measurement error

We now consider that the processes X (i)

t , for i = 1, . . . , I, are observed with a measurement error. Let y (i) j denote the j-th observation of subject i at time t ij . The model is the following:

y (i) j = h(X (i) tij ) + g(X (i) tij , σ)ε ij , ε ij ∼ N (0, I q ) dX (i) t = µ(X (i) t , t, φ i , η)dt + Γ(X (i) t , η)dB (i) t , X (i) 0 = X 0 (φ i ) φ i ∼ p Φ (•; θ pop ) (27)
The parameters to be estimated are θ = (θ pop , η, σ). Let us denote y = (y

(i) 1:ni ) i=1.
..I , the whole vector of data. The likelihood is:

L(θ; y) = I i=1 Φ X p Y (y (i) 1:ni |X (i) 1:ni , σ)p X (X (i) 1:ni |φ i , η)dX (i) 1:ni p Φ (φ i |θ pop )dφ i (28)
where p Y (y

(i) 1:ni |X (i) 1:ni , σ) = ni j=1 p(y (i) j |X (i) j , σ)
is the conditional density of the observations y given the diffusion X. We also denote φ = (φ 1 , . . . , φ I ).

Difficulties pointed out with direct observations of X (i) t still hold but are exacerbated by the necessity to integrate over the hidden trajectories X Donnet and Samson (2008, submitted); Delattre and Lavielle (submitted), are gathered in Section 5.3.3.

Bayesian inference

In a Bayesian inference, a prior distribution is set on the parameters of interest θ ∼ π(•). The aim is to estimate the posterior distribution of θ given the observation y, which is given by the Bayes formula p(θ|y) ∝ f (y; θ)π(θ). Estimators of θ are then obtained with the posterior mean E(θ|y)

or the median of the posterior distribution, for instance. Except in very few cases, the posterior distribution p(θ|y) has no explicit expression. A standard strategy is to generate a sample from that distribution of interest, using a Monte Carlo Markov Chain (MCMC) algorithm [START_REF] Robert | Monte Carlo statistical methods[END_REF]. This approach has been studied for modeling growth curves with a population Gompertzian SDE by [START_REF] Donnet | Bayesian analysis of growth curves using mixed models defined by stochastic differential equations[END_REF].

The principle of the MCMC algorithm is to generate a Markov Chain whose marginal stationary distribution is the distribution of interest p(θ|y). In model ( 27), the simplest MCMC algorithm is to alternatively generate, at iteration k:

[1]. φ (k) |X (k-1) , y, θ (k-1) [2]. X (k) |φ (k) , y, θ (k-1) [3]. θ (k-1) |X (k) , φ (k) , y,
with stationary distribution p(φ, X, θ|y). One can prove that the stationary distribution of the Markov Chain produced by this algorithm is p(φ, X, θ|y) [START_REF] Robert | Monte Carlo statistical methods[END_REF]. As a consequence, after a high number of iterations, the sampled parameters θ (k) are assumed to be distributed under the marginal posterior distribution p(θ|y). However, in population PK-PD SDE, the conditional distributions p(φ|X, θ, y), p(X|φ, y, θ) and p(θ|X, φ, y) are not explicit. Thus, the direct simulation is replaced by a Metropolis-Hastings algorithm, that is, a candidate is simulated with a proposal distribution q and accepted with a probability whose expression ensures the stationarity of the distribution of interest p(φ, X, θ|y). For instance, the Metropolis-Hastings algorithm for step [1] which generates φ (k) |X (k-1) , y, θ (k-1) at iteration k of the MCMC algorithm is

• Propose a candidate φ c with a proposal distribution φ c ∼ q(•|φ (k-1) )

• Compute the probability of acceptation:

ρ(φ c |φ (k-1) ) = min 1, p(φ c |θ (k) , X (k-1) , y) p(φ (k-1) |θ (k) , X (k-1) , y) q(φ (k-1) |φ c ) q(φ c |φ (k-1) ) • Set φ (k) = φ c with probability ρ(φ c |φ (k-1) ) φ (k-1) with probability 1 -ρ(φ c |φ (k-1) )
Similar algorithms are proposed to simulate steps This algorithm can be used as soon as the transition density p(X (k-1) |φ (k-1) , θ c ) of the SDE is explicit. If not, [START_REF] Donnet | Bayesian analysis of growth curves using mixed models defined by stochastic differential equations[END_REF] propose to perform the Bayesian inference on an approximate model issued from the Euler-Maruyama approximation of the solution of the SDE. They study the distance between the approximate posterior distribution p h (θ|y), where h is the step size of the Euler-Maruyama approximation, and the exact one p(θ|y). This distance is bounded by h. Following the idea of recursive conditioning used for standard SDEs, the likelihood L(θ; y) (28) can be re-written as:

Linearisation and extended Kalman filter

L(θ; y) = I i=1 L i (θ; y (i) 1:ni ) = I i=1 Φ   ni j=1 p(y (i) j |y (i) 1:j-1 , φ i , σ, η)   p Φ (φ i |θ pop )dφ i where y (i) 1:j = (y (i) 1 , . . . , y (i) j ) and y (i) 0 = x 0 (φ i ) by convention. As in Section 4.2.2, the EKF approx- imates p(y (i) j |y (i) 1:j-1 , φ i , η) by a Gaussian distribution, whose mean m (i) j|1:j-1 and variance R (i) j|1:j-1
are computed recursively (see equation ( 22)) and depend on φ i , σ and η. The approximate likelihood for subject i is: φi) dφ i where i (φ i ) is the approximate conditional log-density of the random effects given the observation (recall that m (i) j|1:j-1 and R (i) j|1:j-1 depend on φ i ). Then i is approximated by a second-order Taylor expansion, where the expansion is made around the value φi that maximizes i (FOCE algorithm).

L EKF i (θ; y (i) 1:ni ) = Φ ni j=1 exp{-1 2 (y (i) j -m (i) j|1:j-1 ) T (R (i) j|1:j-1 ) -1 (y (i) j -m (i) j|1:j-1 )} |2πR (i) j|1:j-1 | p Φ (φ i |θ pop )dφ i = Φ e i(
At the value φi , the first derivative of i is zero. So finally, the likelihood approximated by EKF and by the Taylor expansion reduces to

L EKF,F OCE (θ, y) = I i=1 1 √ 2π e i( φi) |-H( i )| -1/2
where H( i ) is the Hessian of i . The Hessian is not computed exactly, but with an approximation similar to FOCE (see [START_REF] Pinheiro | Mixed-effect models in S and Splus[END_REF], for more details). Once that approximation has been done, the maximization of L EKF,F OCE (θ, y) is performed by FOCE. This methodology suffers from the drawbacks inherent to EKF and FOCE, and no theoretical convergence has been proved.

This algorithm is implemented in NONMEM [START_REF] Tornøe | Stochastic differential equations in NONMEM: implementation, application, and comparison with ordinary differential equations[END_REF], in a Matlab package [START_REF] Mortensen | A matlab framework for estimation of NLME models using stochastic differential equations. applications for estimation of insulin secretion rates[END_REF] and in the R package PSM [START_REF] Klim | Population stochastic modelling (PSM) -An R package for mixed-effects models based on stochastic differential equations[END_REF].

EM algorithm for mixed effects models defined by SDE

Another approach to estimate parameters of mixed effects models is based on the Expectation-Maximization (EM) algorithm, originally proposed by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Stochastic versions of EM have been adapted to nonlinear mixed models [START_REF] Walker | An EM algorithm for non-linear random effects models[END_REF][START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]. Donnet and Samson (2008, submitted) and Delattre and Lavielle (submitted) propose to extend it to population SDEs. This approach has been applied to a one-compartment PK model (2) with three random effects and measurement noise. We now detail their method.

The EM algorithm is a way to circumvent the problem of integrating the likelihood (28) over the non-observed data X and φ, in situations where the maximization of the conditional expectation of the likelihood of the complete data (y, X, φ) is possible:

Q(θ|θ ) = E [log p(y, X, φ; θ)|y; θ ] .
EM is an iterative procedure: at the -th iteration, given the current value θ -1 of the parameters, the E-step evaluates

Q (θ) = Q(θ | θ -1 ) while the M-step updates θ -1 by maximizing Q (θ).
This algorithm is not explicit for population SDEs: the E-step is not explicit and the conditional distribution p(X, φ|y; θ) is not explicit. A stochastic version of the EM algorithm has been proposed by [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF], the Stochastic Approximation EM algorithm (SAEM), for nonlinear mixed model with deterministic regression function by [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] and for mixed model based on SDEs by Donnet and Samson (2008, submitted) and Delattre and Lavielle (submitted).

The SAEM algorithm for population SDE (27) proceeds as follows. The E-step is divided into a simulation step (S-step) of the non-observed data (X ( ) , φ ( ) ) with the conditional distribution p(X, φ |y; θ -1 ) and a stochastic approximation step (SA-step):

Q (θ) = Q -1 (θ) + α log p(y, X ( ) , φ ( ) ; θ) -Q -1 (θ) ,
where (α ) ∈N is a sequence of positive numbers decreasing to zero. The simulation step can not be performed exactly, for the reasons exposed in Section 5.3.1. When the transition density p X (x, t-s|x s , φ, η) is explicit, [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] propose to use a MCMC algorithm similar to the one used in the Bayesian framework. An improvement has also been proposed more recently When the transition density is not explicit, [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] propose to approximate the solution of the SDE with the Euler-Maruyama scheme. More precisely, if the time interval between two observation points is small, the transition density p X (x, t -s|x s , φ, η) can be approximated by a Gaussian distribution. If the time intervals between two observation instants are too large to obtain a good approximation of the transition density, a natural approach is to introduce a set of auxiliary latent data points between every pair of observations, as first proposed by [START_REF] Pedersen | A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations[END_REF]. Under general assumptions, [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] prove that the SAEM-MCMC algorithm theoretically converges to a maximum of the likelihood of the Euler-Maruyama approximate model. The error between the estimator issued from this approximate model and the "true" model ( 27) can be bound by a function of the step size of the Euler-Maruyama approximation. However, in practice the use of a Euler-Maruyama scheme implies the introduction of a large amount non-observed data, which can slow the convergence rate.

An alternative has been proposed by Delattre and Lavielle (submitted) with an approximation of the E-step of EM by the EKF algorithm. Their SAEM algorithm is based on the simulation of only the individual parameters φ i under the conditional distribution p(φ i |y (i) ; θ), the diffusion trajectories are not simulated and are directly integrated out in the conditional distribution p(φ i |y (i) ; θ).

Since direct simulation of this conditional distribution is generally impossible, a MCMC algorithm is used with a Metropolis-Hastings approach. However, the acceptance probability computed in Metropolis-Hastings requires the knowledge of the expression of p(y (i) |φ i , σ) which has no close form. Delattre and Lavielle (submitted) propose a Gaussian approximation of this conditional likelihood, based on the EKF. The implementation of the EKF is similar to (22) but with a higher order Taylor expansion for the prediction equations. This algorithm is much more less computationally intensive than those proposed by Donnet and Samson (2008, submitted) but the theoretical convergence of their estimator is not proved.

Conclusion and discussion

Due to the specific nature of PK/PD data, standard estimation methods for SDE models proposed in the financial literature can not be directly applied. In this paper, we restrict our review to PK/PD oriented estimation methods, eliminating those relying on continuous observation of the SDE. Moreover, dynamical systems used in biology are quickly of high complexity as soon as a precise description of the biological processes of interest is considered. As a consequence, ad hoc estimation methods have been proposed in the literature, sometimes suffering from a lack of theoretical justifications but always with practical objectives.

In most papers, estimators are obtained by maximizing the likelihood. If the transition density of the SDE is explicit, exact MLE can be computed, using Kalman filter if SDE is observed with measurement noise. Otherwise, if the transition density is not explicit, the likelihood has to be approximated. We can find in the literature easy to implement methods with theoretical validity (Hermite expansion) when the SDE is directly observed, easy to implement methods but without theoretical validity (extended Kalman filter) and methods with theoretical validity but computationally intensive (Monte Carlo approximation).

In a population context, there exist four estimation methods: Hermite expansion of the marginal likelihood, the combination of FOCE with the extended Kalman filtering and the use of a stochastic version of the EM algorithm and Bayesian approach.

Note that Hermite expansion can not be directly extended when the diffusion is observed with noise whereas Kalman filter, Extended Kalman filter and FOCE methods are strictly designed for models with observation noise. Methods relying on the EM algorithm can be applied to models with or without noise, provided the complete likelihood belongs to the exponential family. Bayesian approach is finally the most adaptable method since it can be applied to any model (one or several trajectories, with or without noise observation) as soon as the transition density is explicit. If not, a Euler approximation has to be performed and the introduction of set of auxiliary latent data points between every pair of observations can dramatically jeopardize the convergence rate of the algorithm of simulation. We propose in Table 6 a resume of existing methods for population SDE models.

[Table 1 about here.]

Ferrante

  et al. (2003) consider the stochastic one compartment PK model (2). They introduce the ideal maximum likelihood estimator (MLE) based on continuous-time observed trajectory

  Ditlevsen et al. (2005) propose an estimator based on the spectral density p(u, θ), where u is a frequency in the Fourier domain, the spectral density characterizing uniquely the stochastic process. An estimation of θ is obtained by minimizing the distance between the theoretical spectral density p(u, θ) of the SDE and the estimated spectral density p(u) obtained from the data by the periodogram at the Fourier frequency u. The exact spectral density p(u, θ) of the PD model having no explicit form, the authors propose to approximate it by simulation. If K independent trajectories are simulated with a Euler-Maruyama approximation of the SDE for a fixed parameter value θ, a periodogram pk (u, θ) can be estimated for each trajectory k and the spectral density p(u, θ) is then approximated by Monte Carlo:

  estimator has been used in a PK context byTornøe et al. (2004a) with a stochastic PK/PD model describing the dynamic of insulin and glucose from an euglycaemic clamp study. The PK model of insulin kinetic is a single compartment with linear rates. The insulin effect is modeled by adding an effect compartment. The PK part contains three equations describing the amount of insulin remaining to be absorbed from the tissue, the amount of insulin in the central compartment and the amount of insulin in the effect compartment. The PD of the amount of infused glucose is linked to the insulin concentration with a sigmoidal E max model. The estimation results show that the diffusion coefficients are all estimated to zero, which indicates that the measured data can eventually be modeled by a deterministic PK/PD model.Another example is given by[START_REF] Cuenod | Parameter estimation and change-point detection from dynamic contrast enhanced mri data using stochastic differential equations[END_REF]. They consider the stochastic two-compartment PK model (4) to describe tissue micro-vascularization and angiogenesis via Dynamic Contrast Enhanced Imaging (DCE-imaging) techniques. Results illustrate that the stochastic PK model is more stable than its deterministic version.

  . Suggested solutions are derived from the ones presented in Section 4.2 and combined with estimation methods developed for nonlinear mixed effects models. We first present, in Section 5.3.1, a Bayesian approach, proposed by Donnet et al. (2010). Then we detail how the extended Kalman Filter (EKF) can be coupled with the First Order Conditional Estimate (FOCE) algorithm to estimate θ, as first proposed by Overgaard et al. (2005); Tornøe et al. (2005) (Section 5.3.2). Finally, estimation methods based on the Expectation-Mazimization (EM) algorithm, proposed by

  [2] and[3]. Metropolis-Hastings does not require the exact expression of the targeted distribution. Thanks to the Bayes formula, only the expression of p Y (y(i) 1:ni |X (i) 1:ni , σ), p X (X (i)1:ni |φ i , η) and p Φ (φ i |θ pop ) are used. We refer the reader to[START_REF] Donnet | Bayesian analysis of growth curves using mixed models defined by stochastic differential equations[END_REF] for more details.

Overgaard

  et al. (2005); Tornøe et al. (2005) propose to adapt the well-known First Order Conditional Estimate (FOCE) algorithm to population SDEs. FOCE has been introduced by Lindstrom and Bates (1990) for nonlinear mixed models with deterministic regression function. The algorithm is based on a local linearization of the regression function, leading to a linear mixed model, for which estimation of parameters is more easy. Overgaard et al. (2005); Tornøe et al. (2005) use this idea, coupled with the extended Kalman filter (EKF) already mentioned in Section 4.2.2. Their methodology has been used in a collection of papers, see for instance Mortensen et al. (2007); Overgaard et al. (2007); Klim et al. (2009);[START_REF] Berglund | Investigations of a compartmental model for leucine kinetics using non-linear mixed effects models with ordinary and stochastic differential equations[END_REF]. We detail their approach below.

by

  Donnet and Samson (submitted) based on the Particle MCMC (PMCMC) algorithm proposed by[START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF], which combined the MCMC algorithm with particle filter techniques.The theoretical convergence of both algorithms SAEM-MCMC and SAEM-PMCMC have been established if the complete log-likelihood belongs to the exponential family (with respect to the parameter θ).

Laplace approximation. When the individual parameters are of dimension d greater than 1, [START_REF] Picchini | Practicle estimation of high dimensional stochastic differential mixed-effects models[END_REF] use a Laplace approximation of the likelihood. More precisely, the log-likelihood log L i (θ;

where φi = arg max φi log p X (X

denotes the determinant of a matrix and H(φ i |θ) is the Hessian matrix:

Then the approximate likelihood is numerically optimized. This methods requires the computation of the I Hessian matrices H(φ i |θ) at each step of the optimization, which can be burdensome. The authors give in their paper various numerical and practical solutions to tackle that point.

Obs