Keywords: ARMA process, EM algorithm, Hidden Markov Model, Kalman filter, Maximum likelihood estimation, Ornstein-Uhlenbeck process, Partial observations

We consider a bidimensional Ornstein-Uhlenbeck process to describe the tissue microvascularisation in anti-cancer therapy. Data are discrete, partial and noisy observations of this stochastic differential equation (SDE).

Our aim is the estimation of the SDE parameters. We use the main advantage of a one-dimensional observation to obtain an easy way to compute the exact likelihood using the Kalman filter recursion. We also propose a recursive computation of the gradient and hessian of the log-likelihood based on Kalman filtering, which allows to implement an easy numerical maximisation of the likelihood and the exact maximum likelihood estimator (MLE). Furthermore, we establish the link between the observations and an ARMA process and we deduce the asymptotic properties of the MLE. We show that this ARMA property can be generalised to a higher dimensional underlying Ornstein-Uhlenbeck diffusion. We compare this estimator with the one obtained by the well-known EM algorithm on simulated data.

Introduction

Stochastic continuous-time models are a useful tool to describe biological or physiological systems based on continuous evolution (see e.g. [START_REF] Ditlevsen | Stochastic vs. deterministic uptake of dodecanedioic acid by isolated rat livers[END_REF], Ditlevsen et al. (2005), [START_REF] Picchini | Modeling the euglycemic hyperinsulinemic clamp by stochastic differential equations[END_REF]). The biological context of this work is the modeling of tissue microvascularisation in anticancer therapy. This microcirculation is usually modeled by a bidimensional deterministic differential system which describes the circulation of a contrast agent between two compartments (see [START_REF] Brochot | Evaluation of antiangiogenic treatment effects on tumors' microcirculation by bayesian physiological pharmacokinetic modeling and magnetic resonance imaging[END_REF] and appendix A).

However, this deterministic model is unable to capture the random fluctuations observed along time. In this paper, we consider a stochastic version of this system to take into account random variations around the deterministic solution by adding a Brownian motion on each compartment. This leads to a bidimensional stochastic differential equation (SDE) defined as:

dP (t) = (αa(t) -(λ + β)P (t) + (k -λ)I(t))dt + σ 1 dW 1 (t) dI(t) = (λP (t) -(k -λ)I(t))dt + σ 2 dW 2 (t) (1)
where P (t) and I(t) represent contrast agent concentrations in each compartment, a(t) is an input function assumed to be known, α, β, λ and k are unknown positive parameters, W 1 and W 2 are two independent Brownian motions on R, and σ 1 , σ 2 are the constant diffusion terms. We assume that (P (0), I(0)) is a random variable independent of (W 1 , W 2 ). In our biological context, only the sum S(t) = P (t) + I(t) can be measured. So ( 1) is changed into: dS(t) = (αa(t) -βS(t) + βI(t)) dt + σ 1 dW 1 (t) + σ 2 dW 2 (t)

dI(t) = (λS(t) -kI(t)) dt + σ 2 dW 2 (t) (2)
Noisy and discrete measurements (y i , i = 0, . . . , n) of S(t) are performed at times t 0 = 0 < t 1 < . . . < t n = T . The observation model is thus:

y i = S(t i ) + σε i , ε i ∼ N (0, 1)
where (ε i ) i=0,...,n are assumed to be independent and σ is the unknown standard deviation of the Gaussian noise. To evaluate the effect of the treatment on a patient, it is of importance to have a proper estimation of all unknown parameters from this data set. The aim of this paper is to investigate this problem both theoretically and numerically on simulated data.

Parametric inference for discretely observed general SDEs has been widely investigated. [START_REF] Genon-Catalot | On the estimation of the diffusion coefficient for multi-dimensional diffusion processes[END_REF] and [START_REF] Kessler | Estimation of an ergodic diffusion from discrete observations[END_REF] propose estimators based on minimization of suitable contrasts and study the asymptotic distribution of these estimators when the sampling interval tends to zero as the number of observations tends to infinity. For fixed sampling interval, [START_REF] Bibby | Martingale estimation functions for discretely observed diffusion processes[END_REF] propose martingale estimating functions. In a biological context, Ditlevsen et al. (2005) propose an estimation method based on simulation. [START_REF] Picchini | Maximum likelihood estimation of a time-inhomogeneous stochastic differential model of glucose dynamics[END_REF] propose estimators based on the Hermite expansion of the transition densities. When combining the case of discrete, partial and noisy observations, parameter estimation is a more delicate statistical problem.

In this context, it is classical to estimate the unobserved signal (filtering) (see e.g. [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]). However, our aim is the estimation of SDE parameters. In this paper, we use the main advantage of a one-dimensional observation y and the Gaussian framework of all distributions to obtain an easy way to compute the exact likelihood. For this, we solve and discretize the SDE (2). Then we use the Kalman filter recursion to compute the exact likelihood as proposed by [START_REF] Pedersen | Statistical analysis of gaussian diffusion processes based on incomplete discrete observations[END_REF] and implemented in the Danish Technical University project CTSM. We also obtain a recursive computation of the exact gradient and hessian of the log-likelihood based on Kalman filtering, which allows us to implement an easy numerical maximisation of the likelihood using a gradient method and to compute the exact maximum likelihood estimator. The exact observed Fisher information matrix is also directly obtained. As our model is a hidden Markov model, we develop a second approach based on the EM algorithm, which is widely used in this context since the so-called complete likelihood (observed, unobserved) is generally explicit whereas the exact likelihood (observed) is generally not explicit. This method has been first proposed by [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF] and [START_REF] Segal | A new method for evaluating the loglikelihood gradient, the hessian, and the fisher information matrix for linear dynamic systems[END_REF]. [START_REF] Segal | A new method for evaluating the loglikelihood gradient, the hessian, and the fisher information matrix for linear dynamic systems[END_REF] claim that the EM algorithm is computationally more efficient that the Kalman filters. Thus we compare the EM algorithm and the Kalman filter approach in our context.

In Section 2, we study the SDE. We detail in Section 3 the computation of the exact likelihood, the score and hessian functions. We present the EM method in Section 4. In Section 5, we establish the link between the observations and an ARMA process. This allows to deduce the asymptotic properties of the maximum likelihood estimator. Section 6 contains numerical results based on simulated data. This allows to compare the two estimation methods. Appendix A describes briefly the biological background. Appendix B, C, D, E and F contain some proofs and auxiliary results. In particular, the ARMA property can be generalised to a higher dimensional underlying Ornstein-Uhlenbeck diffusion.

Study of the stochastic differential equation

Introducing U (t) = (S(t), I(t)) where denotes the transposed matrix, (2) can be written in a matrix form:

dU (t) = (F (t) + G U (t))dt + ΣdW (t), U (0) = U 0 y i = J U (t i ) + σε i
where J = (1 0) and

F (t) = αa(t) 0 , G = -β β λ -k , dW (t) = dW 1 (t) dW 2 (t) , Σ = σ 1 σ 2 0 σ 2
The process (U (t)) is a bidimensional Ornstein-Uhlenbeck diffusion, which can be explicitly solved. From the biological context (see Appendix A), the parame-ters satisfy β, k, λ > 0 and λ < k. This implies that G is diagonalizable with two distinct negative eigenvalues. Setting d = (β -k) 2 + 4βλ > 0, the eigenvalues of G are distinct and equal to:

µ 1 = -(β + k) - √ d 2 and µ 2 = -(β + k) + √ d 2
The diagonal matrix D of eigenvalues and the matrix P of eigenvectors are:

D = µ 1 0 0 µ 2 , P = 1 1 β-k- √ d 2β β-k+ √ d 2β with D = P -1 GP.
Proposition 1 Let X(t) = P -1 U (t) and Γ = (Γ kj ) 1≤k,j≤2 = P -1 Σ. Then, for t, h ≥ 0, we have:

X(t + h) = e Dh X(t) + B(t, t + h) + Z(t, t + h) (3)
where t+h) t+h t e -Ds ΓdW s .

B(t, t + h) = e D(t+h) t+h t e -Ds P -1 F (s)ds (4)
Z(t, t + h) = e D(
(5)

Therefore, the conditional distribution of X(t + h) given X(s), s ≤ t is

N 2 e Dh X(t) + B (t, t + h) , R (t, t + h)
where

R(t, t + h) = e (µ k +µ l )h -1 µ k + µ l (ΓΓ ) kl 1≤k,l≤2 (6)
If a(t) ≡ c ≥ 0 with c a constant, (X(t)) has a Gaussian stationary distribution with mean equal to

M = -D -1 P -1 F
and covariance matrix equal to

V = 1 -(µ k + µ l ) (ΓΓ ) kl 1≤k,l≤2
Proof. See Appendix B.

Parameter estimation by maximum likelihood

Our aim is to estimate the unknown parameters α, β, λ, k, σ 1 , σ 2 and σ from observations y 0:n = (y 0 , . . . , y n ). As the law of ((X(t)), ε i , i = 0, . . . , n) is Gaussian, the likelihood of y 0:n can be explicitly evaluated. However, the direct maximization of this likelihood requires the inversion of a matrix of dimension 2(n + 1) × 2(n + 1) (the covariance matrix of (X(t i ))). This inversion can be numerically instable. In this section, we present an alternative method for the computation of the exact likelihood based on Kalman filtering, which does not require any matrix inversion. This is due to the fact that data are onedimensional. Moreover, it is worth stressing that we need not come back to the initial process (U (t)) for computing the likelihood. Indeed, as (U (t)) is not observed, we can use either (U (t)) or any other transformation of (U (t)) even involving unknown parameters. As (X(t)) is simpler, we consider the following transformed model:

dX(t) = (DX(t) + P -1 F (t))dt + ΓdW t , X(0) = P -1 U 0 = X 0 y i = J P X(t i ) + σε i (7)
Given the particular form of our vector J = (1 0) and the fact that the eigenvectors can be chosen up to a proportionnality constant, we have

H = JP = (1 1).
It is especially interesting for further computations of the gradient and hessian of the likelihood that H does not depend of any unknown parameter. From model ( 7) and ( 3)-( 6), we deduce the following discrete-time evolution system where X i = X(t i ):

X i = A i X i-1 + B i + η i , η i ∼ N (0, R i ) y i = HX i + σε i ( 8
)
where

A i = exp(D(t i -t i-1 )), B i = B(t i-1 , t i ), R i = R(t i-1 , t i ).

Computation of the exact likelihood

This discrete model is a hidden Markov model (HMM) [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]:

(X i
) is a hidden Markov chain on R 2 and, conditionally on (X i ), observations (y i ) are independent. Genon-Catalot and Laredo ( 2006) study the maximum likelihood estimator for general HMM. They specialize the exact likelihood in the case where the unobserved Markov chain is a Gaussian one-dimensional AR(1)

process. We generalize this computation to the case where the unobserved Markov chain is a bidimensional AR(1) process. Let φ denote the vector of unknown parameters and y 0:i = (y 0 , . . . , y i ) the vector of observations until time t i . By recursive conditioning, it is sufficient to compute the distribution of y i given y 0:i-1 :

L(φ, y 0:n ) = p(y 0 ; φ) n i=1 p(y i |y 0:i-1 ; φ).
But the conditional law of y i given y 0:i-1 can be evaluated by

p(y i |y 0:i-1 ; φ) = p(y i |X i ; φ)p(X i |y 0:i-1 ; φ)dX i
Then, as the innovation noise η i of the hidden Markov chain, and the observation noise ε i are Gaussian variables, by elementary computations on Gaussian laws, we are able to get the law of y i given y 0:i-1 if we know the mean and covariance of the Gaussian conditional law of X i given y 0:i-1 . This conditional distribution can be exactly computed using Kalman recursions as proposed by [START_REF] Pedersen | Statistical analysis of gaussian diffusion processes based on incomplete discrete observations[END_REF] and implemented in the Danish Technical University project CTSM.

This computation is described below.

Kalman filter

To ease the reading, the parameter φ is omitted. The Kalman filter is an iterative procedure which computes recursively the following conditional distributions

L(X i |y 0:i-1 ) = N 2 ( X- i , P - i ) (prediction) L(X i |y 0:i ) = N 2 ( Xi , P i ) (f ilter)
where

X- i = E(X i |y 0:i-1 ) and P - i = E((X i -X- i )(X i -X- i ) ) Xi = E(X i |y 0:i ) and P i = E((X i -Xi )(X i -Xi ) )
Let us assume that the law of X 0 is Gaussian. Initial values for the algorithm are:

X 0 ∼ N ( X- 0 , P - 0 )
with X-0 = 0, P - 0 = 0 (from theoretical point of view, one can choose the stationary distribution X-0 = M , P - 0 = V without changes). Next we have the recursive formulae obtained using (8):

X- i = A i Xi-1 + B i , P - i = A i P i-1 A i + R i , i ≥ 1 (9) Xi = X- i + K i (y i -H X- i ), P i = (I -K i H)P - i , i ≥ 0
where [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]).

K i = P - i H (HP - i H + σ 2 ) -1 (see e.g.

Computation of the exact likelihood of the observations

The conditional distribution of y i given y 0:i-1 is Gaussian and one-dimensional.

Let m i (φ) = E φ (y i |y 0:i-1 ) and V i (φ) = V ar φ (y i |y 0:i-1
) denote its mean and variance which are given using ( 8) by

m i (φ) = H X- i , V i (φ) = HP - i H + σ 2
where Xi and P - i depend on φ. The exact likelihood of y 0:n is thus equal to

L(φ, y 0:n ) = n i=0 1 2πV i (φ) exp - 1 2 (y i -m i (φ)) 2 V i (φ) . ( 10
)
Relations (9) imply that there exist two functions F φ and G φ such that

m i (φ) = F φ (m i-1 (φ)), V i (φ) = G φ (V i-1 (φ)) (11)
These iterative relations are used to compute the derivatives of the log-likelihood.

Computation of the maximum likelihood estimator

Pedersen (1994) and the Danish Technical University project CTSM propose to approximate the maximum likelihood estimator (MLE) using a quasi-Newton maximisation method based on the approximation of the gradient and the hessian of the log-likelihood. In our model, we show that it is possible to compute the exact MLE. We use a conjugate gradient method, which relies on the explicit knowledge of the gradient and hessian of the log-likelihood. Both can be exactly computed using formula (10) and observing that the derivatives of m i (φ) and

V i (φ) can be explicitly and recursively computed by derivating formulae (11).

New parametrization

In order to simplify the derivatives of ( 11), from now on, we assume that observation times are equally spaced and set

∆ = t i -t i-1 , ∀i = 1, . . . , n.
Hence we have

A i = A, R i = R.
For the sake of simplicity we set a(t)

≡ c ≥ 0,
where c is a known constant, corresponding to an intravenous injection in our biological framework. Hence we have

B i = B = -(I -A)D -1 P -1 F = (I -A)M . Set Z i = X i -M and m = HM . Therefore, model (8) becomes Z i = AZ i-1 + η i , η i ∼ N (0, R) y i = HZ i + m + σε i (12)
The exact likelihood (10) of y 0:n is thus equal to

L(φ, y 0:n ) = n i=0 1 2πV i (φ) exp - 1 2 (y i -m -H Ẑ- i (φ)) 2 V i (φ) . ( 13
)
Moreover, instead of φ = (α, β, λ, k, σ 1 , σ 2 , σ 2 ), we propose a new parametrization fitted with the discretization. We consider θ = (θ 1 , θ 2 , θ 3 , θ 4 , θ 5 , θ 6 ) where θ i = e µi∆ , i = 1, 2 and θ 3 , θ 4 and θ 5 are explicit functions of µ 1 , µ 2 , σ 1 , σ 2 and ∆ such that

A = A(θ) = θ 1 0 0 θ 2 and R = R(θ) = θ 3 θ 5 θ 5 θ 4
and θ 6 = m. We set ϑ = (θ, σ 2 ). Our aim is to maximize the likelihood L(ϑ, y 0:n )

with respect to ϑ. Given an estimation θ, φ can be obtained by solving numerically the equation f ( φ) = θ where f is the mapping φ → f (φ) = ϑ (see Appendix C for details). Note that, later on, we will see that only six out of the seven parameters can be consistently estimated.

Computation of the exact gradient and hessian of the loglikelihood

Let W i (ϑ) = y i -ϑ 6 -H Ẑ- i (
ϑ) and l 0:i (ϑ) = log L(ϑ, y 0:i ). Using (10), it comes:

l 0:i (ϑ) = l 0:i-1 (ϑ) - 1 2 log(2πV i (ϑ)) - 1 2 W i (ϑ) 2 V i (ϑ) . ( 14
)
Thus for i = 1, . . . , n, q = 1, . . . , 7:

∂l 0:i ∂ϑ q (ϑ) = ∂l 0:i-1 ∂ϑ q (ϑ) - 1 2 1 V i (ϑ) ∂V i ∂ϑ q (ϑ) - W i (ϑ) V i (ϑ) ∂W i (ϑ) ∂ϑ q + 1 2 W i (ϑ) 2 V i (ϑ) 2 ∂V i (ϑ) ∂ϑ q (15)
where

∂V i (ϑ) ∂ϑ q = H ∂P - i (ϑ) ∂ϑ q H , 1 ≤ q ≤ 6, ∂V i (ϑ) ∂σ 2 = H ∂P - i (ϑ) ∂σ 2 H + 1 ∂W i (ϑ) ∂ϑ q = - ∂m ∂ϑ q -H ∂ X- i (ϑ) ∂ϑ q , 1 ≤ q ≤ 7
Furthermore, the derivatives of Xi (ϑ) and P - i (ϑ) can be obtained using Kalman recursions (see appendix D). With a more cumbersome computation, second order derivatives of l 0:n (ϑ) can be analogously deduced from (15). For i = 1, . . . , n, q, r = 1, . . . , 7,

∂ 2 l 0:i ∂ϑ r ∂ϑ q (ϑ) = ∂ 2 l 0:i-1 ∂ϑ r ∂ϑ q (ϑ) - 1 2 1 V i (ϑ) ∂ 2 V i ∂ϑ r ∂ϑ q (ϑ) + 1 2 1 V 2 i (ϑ) ∂V i ∂ϑ r (ϑ) ∂V i ∂ϑ q (ϑ) - 1 2 2 W i (ϑ) V i (ϑ) ∂ 2 W i (ϑ) ∂ϑ r ∂ϑ q - W i (ϑ) 2 V i (ϑ) 2 ∂ 2 V i (ϑ) ∂ϑ r ∂ϑ q (16) - ∂W i (ϑ) ∂ϑ r ∂W j (ϑ) ∂ϑ q 1 V i (ϑ) - W i (ϑ) V i (ϑ) 2 ∂W i (ϑ) ∂ϑ r ∂V i (ϑ) ∂ϑ q + ∂W i (ϑ) ∂ϑ q ∂V i (ϑ) ∂ϑ r W i (ϑ) V i (ϑ) 2 - W i (ϑ) 2 V i (ϑ) 4 ∂V i (ϑ) ∂ϑ r V i (ϑ) ∂V i (ϑ) ∂ϑ q
where (see appendix D for details)

∂ 2 V i (ϑ) ∂ϑ r ∂ϑ q = H ∂ 2 P - i (ϑ) ∂ϑ r ∂ϑ q H and ∂ 2 W i (ϑ) ∂ϑ r ∂ϑ q = -H ∂ 2 X- i (ϑ) ∂ϑ r ∂ϑ q .
Hence, we obtain an explicit expression of -∂ 2 l0:n ∂ϑr∂ϑq (ϑ) 1≤q,r≤7

.

Maximisation of the exact likelihood

To compute the maximum likelihood estimator, the conjugate gradient algorithm is applied to minimize l - 0:n (ϑ) = -l 0:n (ϑ) (see Stoer and Bulirsch (1993)). Let ∇l -denote the gradient of l - 0:n and Hess l -its hessian evaluated by ( 15)-( 16). Starting with an arbitrary initial vector ϑ 0 , we set as descent direction u 0 = ϑ 0 . At iteration k, given ϑ k and u k , the parameter and descent direction are updated by

ϑ k+1 = ϑ k - u k , ∇l -(ϑ k ) u k , Hess l -(ϑ k )u k u k , u k+1 = -∇l -(ϑ k+1 ) + ∇l -(ϑ k+1 ) ∇l -(ϑ k ) u k .
Classical stopping conditions are used. The sequence (ϑ k ) k converges towards the maximum of the likelihood l 0:n (ϑ).

Parameter estimation by Expectation Maximization algorithm

An alternative method to estimate ϑ = (θ, σ 2 ) is the Expectation Maximization (EM) algorithm, proposed by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], see also [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF] and [START_REF] Segal | A new method for evaluating the loglikelihood gradient, the hessian, and the fisher information matrix for linear dynamic systems[END_REF]. The EM algorithm is a classical approach to estimate parameters of models with non-observed or incomplete data, especially it is widely used for HMMs. In our case, the non-observed data are the (Z i )'s, the complete data are the (y i , Z i )'s. The principle is to maximize

ϑ → Q(ϑ|ϑ * ) = E(log p(y 0:n , Z 0:n ; ϑ)|y 0:n ; ϑ * )
with Z 0:n = (Z 0 , . . . , Z n ). This is often easier than the maximization of the observed data log-likelihood since the log-likelihood of the complete data is generally simpler. Moreover, according to Wu (1983), as our model is an exponential family, the EM estimate sequence (ϑ k ) k converges towards a (local) maximum of the data likelihood. The EM algorithm uses two steps: the Expectation step (E-step) and the Maximization step (M-step). Starting with an initial value

(ϑ 0 ), the k-th iteration is • E-step: evaluation of Q k (ϑ) = Q(ϑ | ϑ k ) • M-step: update of ϑ k by ϑ k+1 = arg max Q k (ϑ).
In our model, function Q has an explicit expression. Recall that we have assumed that observations times are equally spaced and that a(t) ≡ c. For simplicity, we also set for the initial variable Z 0 = 0 ( Ẑ-0 = 0 and P - 0 = 0). The complete data log-likelihood is thus equal to:

log p(y 0:n , Z 0:n ; ϑ) = - n + 1 2 log(2πσ 2 ) - 1 2σ 2 n i=0 (y i -ϑ 6 -HZ i ) 2 - 1 2 n i=1 log(2π|R(ϑ)|) - 1 2 n i=1 (Z i -A(ϑ)Z i-1 ) R(ϑ) -1 (Z i -A(ϑ)Z i-1 ).
Function Q consists in taking the conditional expectation given y 0:n under P ϑ * .

This conditional distribution is the so-called smoothing distribution at ϑ * . In our model, it is Gaussian and characterized by M i|0:n (ϑ * ) = E ϑ * (Z i |y 0:n ) and

Σ i|0:n (ϑ * ) = V ar ϑ * (Z i |y 0:n ), Σ i-1,i|0:n (ϑ * ) = Cov ϑ * (Z i-1 , Z i |y 0:n )
These can be obtained through a forward-backward algorithm (see Appendix E). Thus function Q is equal to:

Q(ϑ|ϑ * ) = - n + 1 2 log(2πσ 2 ) - 1 2σ 2 n i=0 (y i -ϑ 6 -HM i|0:n (ϑ * )) 2 + HΣ i|0:n (ϑ * )H - n 2 log(2π|R(ϑ)|) - 1 2 T r R(ϑ) -1 [C(ϑ * ) -T (ϑ * )A (ϑ) -A(ϑ)T (ϑ * ) + A(ϑ)S(ϑ * )A (ϑ)]
where

T (ϑ * ) = n i=1 Σ i-1,i|0:n (ϑ * ) + M i|0:n (ϑ * )M i-1|0:n (ϑ * ) S(ϑ * ) = n i=1 Σ i-1|0:n (ϑ * ) + M i-1|0:n (ϑ * )M i-1|0:n (ϑ * ) C(ϑ * ) = n i=1 Σ i|0:n (ϑ * ) + M i|0:n (ϑ * )M i|0:n (ϑ * ) .
The matrices A, R, θ 6 and σ 2 are updated as

A(ϑ k ) = diag(T (ϑ k-1 )S -1 (ϑ k-1 )) R(ϑ k ) = 1 n (C(ϑ k-1 ) -T (ϑ k-1 )S -1 (ϑ k-1 )T (ϑ k-1 )) ϑ 6k = 1 n + 1 n i=0 (y i -HM i|0:n (ϑ k-1 )) σ 2 k = 1 n + 1 n i=0 (y i -ϑ 6k-1 -HM i|0:n (ϑ k-1 )) 2 + HΣ i|0:n (ϑ k-1 )H
5 Properties of the exact maximum likelihood estimator in the stationary case Recall that we have assumed that a(t) ≡ c. Moreover in this paragraph we assume that the initial variable X 0 has the stationary distribution N 2 (M, V )

given in Proposition 1. This implies that the joint process (X i , y i ) is strictly stationary. Let (y i ) i∈Z be its extension to a process indexed by Z.

Link with an ARMA model

We generalize the result of [START_REF] Genon-Catalot | Conditional likelihood estimators for hidden Markov models and stochastic volatility models[END_REF] to the bidimensional case and also to the multidimensional case (see Appendix F).

Proposition 2 Let (ỹ i ) = (y i -θ 6 ) define the centered process. The process

(ỹ i ) i∈Z is centered Gaussian and ARMA(2,2). Proof. Evidently (ỹ i ) is centered Gaussian. We easily check that ỹi -(θ 1 + θ 2 )ỹ i-1 + θ 1 θ 2 ỹi-2 = ξ i (17)
where ξ i is defined by

ξ i = HA(θ)η i-1 + Hη i + σε i -(θ 1 + θ 2 )Hη i-1 -(θ 1 + θ 2 )σε i-1 + θ 1 θ 2 σε i-2
As the (η i ) i and (ε i ) i are mutually independent, we get that:

Cov(ξ i , ξ i+k ) = 0, ∀k ≥ 3.
This implies that (ξ i ) is MA(2). Hence the result.

Proposition 3 The spectral density f (u, ϑ) of (ỹ i ) has the explicit form:

f (u, ϑ) = σ 2 + H(ARA + (1 + (θ 1 + θ 2 ) 2 )R)H + 2 cos(u)(HA -(θ 1 + θ 2 )H)RH 1 + (θ 1 + θ 2 ) 2 + θ 2 1 θ 2 2 -2(θ 1 + θ 2 )(1 + θ 1 θ 2 ) cos(u) + 2 cos(2u)θ 1 θ 2 (18) with A = A(θ), R = R(θ). Proof. Let γ(k) = Cov(ξ i , ξ i+k ). Elementary computations show that γ(0) = HARA H + HRH (1 + (θ 1 + θ 2 ) 2 ) + σ 2 (1 + (θ 1 + θ 2 ) 2 + θ 2 1 θ 2 2 ) γ(1) = (HA -(θ 1 + θ 2 )H)RH -σ 2 (θ 1 + θ 2 )(1 + θ 1 θ 2 ) γ(2) = σ 2 θ 1 θ 2 γ(k) = 0 ∀ k ≥ 3
The spectral density (with respect to du 2π ) h(u, ϑ) of

(ξ i ) is h(u, ϑ) = n∈Z γ(n) exp(-inu) = γ(0) + γ(1)2 cos(u) + γ(2)2 cos(2u).
For the AR(2) part, we set:

p(x) = 1-(θ 1 +θ 2 )x+θ 1 θ 2 x 2 (recall that θ 1 , θ 2 < 1). Then f (u, ϑ) = h(u, ϑ) |p(exp(-iu))| 2 = γ(0) + γ(1)2 cos(u) + γ(2)2 cos(2u) 1 + (θ 1 + θ 2 ) 2 + θ 2 1 θ 2 2 -2(θ 1 + θ 2 )(1 + θ 1 θ 2 ) cos(u) + 2 cos(2u)θ 1 θ 2
See [START_REF] Brockwell | Time Series: Theory and Methods[END_REF] for technical details.

The number of parameters which are identifiable on the spectral density is precised by the following proposition Proposition 4 The identifiable quantities are σ 2 , θ 1 +θ 2 and θ 1 θ 2 , and at most two out of three parameters among θ 3 , θ 4 and θ 5 .

When ∆ is small, exactly two out of the three paramaters θ 3 , θ 4 and θ 5 are identifiable.

Proof. See Appendix G.

Asymptotic properties of maximum likelihood estimator

We now deduce the asymptotic properties of the maximum likelihood estimators when (y i ) is stationary.

We denote by ϑ 0 the true value of parameters vector (θ 1 , . . . , θ 6 , σ 2 ). We assume that the parameter set Θ is an open subset of R 7 . We denote by

ϑ -= (θ 1 , . . . , θ 5 , σ 2 ) ∈ Θ -the projection of ϑ on R 6 .
The following result is classical to estimate the parameter ϑ 0 6 (see e.g. in Brockwell and Davis (1991)).

Proposition 5 (Mean estimator) Let ȳn = 1 n n i=1 y i be the empirical mean. Under the assumption of stationarity of (y i ), ȳn → θ 0 6 a.s. as n → ∞. Moreover, √ n(ȳ n -θ 0 6 ) converges in distribution:

√ n(ȳ n -θ 0 6 ) -→ n→∞ N (0, J(ϑ 0 -))
where J(ϑ 0 -) = f (0, ϑ 0 -) = γ(0)+2γ(1)+2γ( 2)

((1-θ1)(1-θ2)) 2
We now consider the centered process (ỹ i ). Consider the two assumptions (which can be checked up to some technicities)

A1 (u, ϑ -) -→ f (u, ϑ -) is a C 3 -function on a neighborhood of [-π, π] × Θ - A2 ϑ --→ f (., ϑ -) is one to one As (ỹ i ) is a ARMA(2,2) process, its spectral density is positive for every (u, ϑ -) ∈ [-π, π] × Θ -.
Proposition 6 (Information matrix) Let l0:n (ϑ -) = log L(ϑ -, ỹ0:n ) be the loglikelihood of the centered process (ỹ 0:n ). Under the assumption A1, we have

P ϑ 0 --a.s. lim n→∞ - 1 n ∂ 2 ∂θ i ∂θ j l0:n ϑ 0 - 1≤i,j≤6 = I(ϑ 0 -) (19)
Proposition 7 (Consistency and asymptotic normality of the MLE)

Let θn be a maximum likelihood estimator of ϑ 0 -based on ỹ0:n . Under the assumptions A1 and A2, θn → ϑ 0 -a.s. as n → ∞. Moreover, if

I(ϑ 0 -) is invertible, √ n( θn -ϑ 0 -) converges in distribution: √ n( θn -ϑ 0 -) -→ n→∞ N (0, I -1 (ϑ 0 -))
Proof. This result may be found e.g. in [START_REF] Brockwell | Time Series: Theory and Methods[END_REF].

Remark 1 As it is done usually, for further numerical considerations, the empirical mean estimator ȳn is plugged in the likelihood for parameter θ 6 in the Kalman-recursion approach. The EM algorithm estimates this parameter during the algorithm iterations.

Simulation study

We compare the performances of the two estimation methods on simulated data sets. The exact maximum likelihood estimators and the EM estimators are com- Identifiability results of Section 5 show that parameters θ 6 and σ 2 are identifiable, and that at most 4 among the five parameters (θ 1 , . . . , θ 5 ) are identifiable.

In our biological context, it is reasonnable to assume that σ 2 is known. Therefore σ 2 is fixed to its true value. One parameter among (θ 1 , . . . , θ 5 ) has to be fixed: we choose to fix θ 5 to its true value and to estimate θ 1 , θ 2 , θ 3 , θ 4 . For the MLE estimation, parameter θ 6 is previously estimated by the empirical mean. Then the other parameters are estimated based on the empirically centered observations.

Results are given in Table 1 for ∆ = 0.2 and Table 2 for ∆ = 0.2 and ∆ = 0.04.

Results obtained on one simulated data set (n = 200, ∆ = 0.2, σ 2 = 3) are plotted in Figure 1 and show the ability of the method to estimate the trajectories (S(t), P (t), I(t)). The EM algorithm is around three times quicker than the MLE algorithm. The mean computation time to estimate parameters of a data set with n = 200 observations is around 4 seconds (CPU time) for the EM algorithm and around 13 seconds (CPU time) for the MLE for the same precision of convergence, on an Apple MacPro 2 × 3 Ghz with 5 Go of RAM.

The Matlab code are given at http://www.mi.parisdescartes.fr/~favetto.

[Figure 1 about here.]

The EM estimates are often less biased for parameters θ 1 and θ 2 than the MLE estimates. Variance parameters (θ 3 , θ 4 ) are estimated with less bias by the MLE method, θ 4 is always estimated with a large bias by the EM algorithm. The standard errors of the EM estimates are lower than those of the MLE estimates.

The standard errors reduce with the increase of the number of observations n.

The bias and standard errors of all parameters decrease with the decrease of the observation noise σ. When ∆ decreases, bias and standard errors decrease for a small observation noise σ = 1. MLE estimates are very satisfactory in this case. The exact Fisher information matrix provides standard errors of the estimates which are close to the empirical ones, especially those obtained with the MLE approach. Note that the theoretical study of the exact MLE allows to deduce the identifiable parameters. On the contrary, the EM algorithm misses completely the problem.

[Table 1 about here.]

[Table 2 about here.]

Conclusion

The Kalman filter is classical in the field of noisy, discretely and partially observed stochastic differential equations. In this paper, we have shown that it can be used for the estimation of the parameters by maximum likelihood. In the particular case of an Ornstein-Uhlenbeck process, this method computes the exact likelihood, its gradient and hessian. We have also shown that the EM algorithm, which is classical in the field of hidden Markov models, combined with a smoother algorithm can be used for the parameter estimation.

We study some theoretical properties of the model. We show that only six out of the seven parameters are identifiable and we deduce the asymptotic properties of the maximum likelihood estimate. We illustrate the two methods on simulated data. The identifiability problem is confirmed on the simulation study: the observed Fisher information matrix computed by the Kalman method is not invertible when we estimate the seven parameters.

The next step of this work is its application to real data in anti-cancer therapy. This work could also be extended to the case of non-Gaussian observation errors. For a unidimensional Markov chain (X i ) observed with non Gaussian errors, [START_REF] Ruiz | Quasi-maximum likelihood estimation of stochastic volatility models[END_REF] proposes a quasi-maximum likelihood estimator based on the Kalman filter and shows the normality asymptotic distribution of this estimator.

This approach can be extended to our bidimensional model.

A Physiological model

We focus on the evaluation of anti-angiogenesis treatments in anti-cancer therapy. These treatments take effect on the vascularization of tissue. The in vivo evaluation of their efficacy is based on the estimation of the tissue microvascularization parameters. The experiment consists in the injection of a contrast agent to the patient, followed by the recording of a medical images sequence which measures the evolution of the concentration of contrast agent along time.

The contrast agent pharmacokinetic is modeled by a bidimensional differential system. The contrast agent pulsates in the plasma and interstitium cells. Let a(t), P (t) and I(t) denote respectively the quantity of contrast agent at time t in the artery, the plasma and the interstitium and 1 -h, V P and V I the volume of artery, plasma and interstitium (h is the hematocrit rate). The initial condition at time t 0 = 0 is P (0) = 0, I(0) = 0. The contrast agent is injected in vein at time t 0 , transits in the artery and arrives in plasma, with a tissue perfusion flow equal to F tp . The contrast agent is eliminated from plasma with the perfusion flow F tp , proportionally to the concentration of contrast agent in plasma. The quantity of contrast agent exchanging from plasma through interstitium is equal to K trans times the concentration of contrast agent in plasma, where K trans is and σ 2 2 are solution of the system

               k = -(µ 1 + µ 2 + β), λ = -( µ1µ2 β + k), (µ 2 + β) 2 σ 2 1 + µ 2 2 σ 2 2 = θ3 , (µ 1 + β) 2 σ 2 1 + µ 2 1 σ 2 2 = θ4 , σ 1 β 2 + σ 1 (µ 1 + µ 2 )β + θ5 + (σ 2 1 + σ 2 2 )µ 1 µ 2 = 0.

D Gradient and hessian of the log-likelihood

The gradient and the hessian of the loglikelihood are computed with explicit recursions. We denote ϑ 7 = σ 2 . The first order derivatives of A(ϑ) are equal to:

∂A(ϑ) ∂ϑ 1 = 1 0 0 0 , ∂A(ϑ) ∂ϑ 2 = 0 0 0 1 and ∂A(ϑ) ∂ϑ q = 0 0 0 0 , q = 3, 4, 5, 6, 7.
for R(ϑ) we get: and for m we get: ∂m ∂ϑq = 0, q = 1, . . . , 5, 7 and ∂m ∂ϑ6 = 1. The second order derivatives of A(ϑ) and R(ϑ) are null. The first order derivatives of Xi (ϑ) and P - i (ϑ) with respect to ϑ q , q = 1, . . . , 7 can be deduced:

∂R(ϑ) ∂ϑ 1 = ∂R(ϑ) ∂ϑ 2 = ∂R(ϑ) ∂ϑ 6 = ∂R(ϑ) ∂ϑ 7 = 0 0 0 0 , ∂R(ϑ) ∂ϑ 3 = 1 0 0 0 , ∂R (ϑ
∂ X- i (ϑ) ∂ϑ q = ∂A(ϑ) ∂ϑ q Xi-1 (ϑ) + A(ϑ) ∂ Xi-1 (ϑ) ∂ϑ q ∂P - i (ϑ) ∂ϑ q = ∂A(ϑ) ∂ϑ q P i-1 (ϑ)A(ϑ) + A(ϑ) ∂P i-1 (ϑ) ∂ϑ q A(ϑ) + A(ϑ)P i-1 (ϑ) ∂A(ϑ) ∂ϑ q + ∂R(ϑ) ∂ϑ q .
Then we get the derivatives of the mean and the covariance of the filter:

∂ Xi-1 (ϑ) ∂ϑ q = ∂ X- i-1 (ϑ) ∂ϑ q + ∂P - i-1 (ϑ) ∂ϑ q H W i-1 (ϑ) V i-1 (ϑ) + P - i-1 (ϑ)H V i-1 (ϑ) ∂W i-1 (ϑ) ∂ϑ q - ∂V i-1 (ϑ) ∂ϑ q W i-1 (ϑ) V i-1 (ϑ) ∂P i-1 (ϑ) ∂ϑ q = I - P - i-1 H H V i-1 (ϑ) ∂P - i-1 (ϑ) ∂ϑ q - ∂P - i-1 (ϑ) ∂ϑ q H V i-1 (ϑ) - P - i-1 (ϑ)H V i-1 (ϑ) 2 ∂V i-1 (ϑ) ∂ϑ q HP - i-1 (ϑ)
The computation of second order derivatives can be deduced. Because second order derivatives of A and R are null, we have for q, r = 1, . . . , 7, :

∂ 2 X- i (ϑ) ∂ϑ q ∂ϑ r = A(ϑ) ∂ 2 Xi-1 (ϑ) ∂ϑ q ∂ϑ r + ∂A(ϑ) ∂ϑ r ∂ Xi-1 (ϑ) ∂ϑ q + ∂A(ϑ) ∂ϑ q ∂ Xi-1 (ϑ) ∂ϑ r ∂ 2 P - i (ϑ) ∂ϑ q ∂ϑ r = ∂A(ϑ) ∂ϑ r ∂P i-1 (ϑ) ∂ϑ q + ∂A(ϑ) ∂ϑ q ∂P i-1 (ϑ) ∂ϑ r A(ϑ) + ∂A(ϑ) ∂ϑ r P i-1 (ϑ) ∂A(ϑ) ∂ϑ q +A(ϑ) ∂ 2 P i-1 (ϑ) ∂ϑ q ∂ϑ r A(ϑ) + ∂P i-1 (ϑ) ∂ϑ r ∂A(ϑ) ∂ϑ q + ∂P i-1 (ϑ) ∂ϑ q ∂A(ϑ) ∂ϑ r + ∂A(ϑ) ∂ϑ q P i-1 (ϑ) ∂A(ϑ) ∂ϑ r
The second order derivatives of the mean and the covariance of Kalman filter are:

∂ 2 Xi-1 (ϑ) ∂ϑ q ∂ϑ r = ∂ 2 X- i-1 (ϑ) ∂ϑ q ∂ϑ r + ∂ 2 P - i-1 (ϑ) ∂ϑ q ∂ϑ r H W i-1 (ϑ) V i-1 (ϑ) + P - i-1 (ϑ) H V i-1 (ϑ) ∂ 2 W i-1 (ϑ) ∂ϑ q ∂ϑ r - ∂V i-1 (ϑ) ∂ϑ q 1 V i-1 (ϑ) ∂W i-1 (ϑ) ∂ϑ r - ∂ 2 V i-1 (ϑ) ∂ϑ q ∂ϑ r W i-1 (ϑ) V i-1 (ϑ) + 2 ∂V i-1 (ϑ) ∂ϑ r ∂V i-1 (ϑ) ∂ϑ q W i-1 (ϑ) V 2 i-1 (ϑ) - ∂V i-1 (ϑ) ∂ϑ r 1 V i-1 (ϑ) ∂W i-1 (ϑ) ∂ϑ q + ∂P - i-1 (ϑ) ∂ϑ r H V i-1 (ϑ) ∂W i-1 (ϑ) ∂ϑ q - ∂V i-1 (ϑ) ∂ϑ q W i-1 (ϑ) V i-1 (ϑ) + ∂P - i-1 (ϑ) ∂ϑ q H V i-1 (ϑ) ∂W i-1 (ϑ) ∂ϑ r - ∂V i-1 (ϑ) ∂ϑ r W i-1 (ϑ) V i-1 (ϑ) ∂ 2 P i-1 (ϑ) ∂ϑ q ∂ϑ r = I - P - i-1 (ϑ)H H V i-1 (ϑ) ∂ 2 P - i-1 (ϑ) ∂ϑ q ∂ϑ r - ∂ 2 P - i-1 (ϑ) ∂ϑ q ∂ϑ r - P - i-1 (ϑ) V i-1 (ϑ) ∂ 2 V i-1 (ϑ) ∂ϑ q ∂ϑ r - ∂P - i-1 (ϑ) ∂ϑ r ∂V i-1 (ϑ) ∂ϑ q + ∂P - i-1 (ϑ) ∂ϑ q ∂V i-1 (ϑ) ∂ϑ r 1 V i-1 (ϑ) + 2P - i-1 (ϑ) V 2 i-1 (ϑ) ∂V i-1 (ϑ) ∂ϑ q ∂V i-1 (ϑ) ∂ϑ r H HP - i-1 (ϑ) V i-1 (ϑ) - ∂P - i-1 (ϑ) ∂ϑ r H - P - i-1 (ϑ)H V i-1 (ϑ) ∂V i-1 (ϑ) ∂ϑ r H V i-1 (ϑ) ∂P - i-1 (ϑ) ∂ϑ q - ∂P - i-1 (ϑ) ∂ϑ q H -P - i-1 (ϑ) H V i-1 (ϑ) ∂V i-1 (ϑ) ∂ϑ q H V i-1 (ϑ) ∂P - i-1 (ϑ) ∂ϑ r

E Smoother algorithm

The Kalman smoother is calculated recursively with a forward-backward algorithm (see e.g. [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]). The forward algorithm is the classical Kalman filter which computes

M i|0:i-1 = Ẑ- i = E(Z i |y 0:i-1 ), Σ i|0:i-1 = P - i = V ar(Z i |y 0:i-1 ), M i|0:i = Ẑi = E(Z i |y 0:i ) and Σ i|0:i = P i = V ar(Z i |y 0:i ). Then, in order to calculate M i|0:n = E(Z i |y 0:n ), Σ i|0:n = V ar(Z i |y 0:n ), Σ i-1,i|0:n =
Cov(Z i-1 , Z i |y 0:n ), one performs the set of backward recursions i = n, n -1, . . . , 1:

J i-1 = Σ i-1|0:i-1 A (Σ i|0:i-1 ) -1 M i-1|0:n = M i-1|0:i-1 + J i-1 (M i|0:n -M i|0:i-1 ) Σ i-1|0:n = Σ i-1|0:i-1 + J i-1 (Σ i|0:n -Σ i|0:i-1 )J i-1
To calculate Σ i-1,i|0:n , we have

Σ n-1,n|0:n = (I -K n H)AΣ n-1|0:n-1
and the following backward recursions, for i

= n -1, n -2, . . . , 1 Σ i-1,i|0:n = Σ i|0:i J i-1 + J i (Σ i,i+1|0:n -AΣ i|0:i )J i-1

F ARMA property of multidimensional process

In our model, (y i ) i∈Z is an ARMA(2,2) process and asymptotic properties of the maximum likelihood estimator are derived. This result can be generalised to the case where X i is p-dimensional under weak assumptions. We consider the model:

y i = HX i + σε i , X i = AX i-1 + η i , X 0 ∼ ν where X i is p-dimensional, A is a diagonal matrix with diagonal coefficients (θ k , k = 1, . . . , p) such that θ k = θ l for k = l and θ i ∈ (0, 1) for i = 1 . . . p, (η i ) i≥0 is a sequence of independent N p (0, R) random variables, H is a (1, p)-matrix and
(ε i ) is a sequence of i.i.d N (0, 1) random variables. Up to a transformation of (X i ), H is equal to H = (1 . . . 1 0 . . . 0) with its first d coordinates equal to 1 and its p -d next coordinates equal to 0 (1 ≤ d ≤ p). Consequently, we observe with additive noise the partial sum of the first d coordinates of X i . Since A is diagonal and θ i ∈ (0, 1) for i = 1 . . . p, the process (X i ) i≥0 admits a stationary distribution ν. Then with X 0 ∼ ν, the process (X i ) i≥0 is stationary. Denote (y i ) i∈Z the extension of y to Z by stationarity.

Proposition 8 The process (y i ) i∈Z is ARMA(d,d).

Proof. Denote S j , j = 1, . . . , d, the j-th symmetric function of θ 1 , . . . , θ d :

S j = 1≤i1<•••<ij ≤d θ i1 . . . θ ij
and S 0 = 1. Define the polynomial P such that P (θ 1 ) = . . . = P (θ d ) = 0:

P (x) = d k=0 (-1) k S k x d-k
Set L the one-lag operator :

Ly i = y i-1 , Lε i = ε i-1 . Set ξ i = P (L)(y i ) = d k=0 (-1) k S k y i-k . By recursive computation for 0 ≤ k ≤ d, we have y i-k = HA d-k X i-d + d-k-1 j=0 HA j η i-k-j + σε i-k
We deduce

ξ i = HP (A)X i-d + d k=0 (-1) k S k   d-k-1 j=0 HA j η i-k-j   + σP (L)(ε i )
But HP (A) = (P (θ 1 ) . . . P (θ d ) 0 . . . 0) = (0 . . . . . . 0). Thus ξ i only depends on (η j , ε j ) j≤i . Therefore (y i ) verifies an AR(d) equation. Moreover, as the (η i ) i and (ε i ) i are mutually independent, we get that:

Cov(ξ i , ξ i+k ) = 0, ∀k ≥ d.
This implies that (ξ i ) is MA(d). Hence the result.

G First order identifiability

The spectral density can be rewritten as

f (u, ϑ) = σ 2 + d(ϑ)e iu + c(ϑ) + d(ϑ)e -iu (1 -(θ 1 + θ 2 )e iu + θ 1 θ 2 e 2iu )(1 -(θ 1 + θ 2 )e -iu + θ 1 θ 2 e -2iu ) where c(ϑ) = H(ARA +(1+(θ 1 +θ 2 ) 2 )R)H and d(ϑ) = (HA-(θ 1 +θ 2 )H)RH . The equality f (u, ϑ) = f (u, ϑ ) ∀u ∈ (0, 2π) implies σ 2 = σ 2 , θ 1 + θ 2 = θ 1 + θ 2 , θ 1 θ 2 = θ 1 θ 2 , c(ϑ) = c(ϑ ), d(ϑ) = d(ϑ )
We deduce that σ 2 , θ 1 + θ 2 , θ 1 θ 2 are identifiable and that at most two of three parameters among θ 3 , θ 4 and θ 5 are identifiable from c(ϑ) and d(ϑ).

We can prove that there are exactly two parameters identifiable when ∆ is small.

Set

g(z, ϑ) = d(ϑ) + c(ϑ)z + d(ϑ)z 2 . so that for z = e iu f (z, ϑ) = σ 2 + z g(z, ϑ) (1 -θ 1 z)(θ 1 -z)(1 -θ 2 z)(θ 2 -z)
Note that the product of the roots of g is equal to 1. Thus, if θ 1 and θ 2 are not roots of g, then no simplification is possible in the expression of f and exactly two out of the three parameters θ 3 , θ 4 , θ 5 are identifiable (from c(ϑ) and d(ϑ)).

µ 1 µ 2 µ 2 -µ 1 ) + O(∆) Hence it comes g(θ 1 ) = O(∆ 2 ) + 4C∆ + o(∆), g(θ 2 ) = O(∆ 2 ) + 4C∆ + o(∆).
with

C = θ 3 + θ 4 + 2θ 5 But C is a positive constant because C = (1 1)R(1 1) .
where R is a covariance matrix thus positive. Hence g(θ 1 ) and g(θ 2 ) are positive when ∆ is small. We have thus a non degenerate ARMA(2,2) process and exactly five parameters are identifiable. obtained with the exact MLE and the EM algorithms and exact standard errors obtained from the Fisher information matrix, evaluated on 1000 simulated data with n = 200 and n = 1000 observations and σ 2 = 1 or σ 2 = 3 (σ 2 and θ 5 fixed to their true values).

  puted as described in Section 3 and Section 4, respectively. Data are simulated using equally spaced observation times (∆ = 0.2) and n = 200 or n = 1000 observations. Values of parameter θ are deduced from values of biological parameters (α, β, λ, k) estimated on real data in Thomassin (2008), and σ 2 1 = 0.5, σ 2 2 = 0.125, c = 50:θ 1 = 0.6, θ 2 = 0.9, θ 3 = 0.7, θ 4 = 0.2, θ 5 = 0.1, θ 6 =20 Two levels of observation noise are used: σ 2 = 1 or σ 2 = 3. Thousand replications are performed for each design (n = 200 or n = 1000 observations, σ 2 = 1 or σ 2 = 3). The influence of the time scale ∆ is evaluated on simulated data with ∆ = 0.04 and n = 1000 observations with parameter values equal to θ 1 = 0.91, θ 2 = 0.99, θ 3 = 0.2, θ 4 = 0.03, θ 5 = 0.01, θ 6 = 20 Thousand replications are performed for each observation noise (σ 2 = 1 or σ 2 = 3) in this case. Mean estimates and their empirical standard errors are computed on the 1000 replications of each design. The exact standard errors obtained from the asymptotic information matrix (computed by Kalman-based recursions) are also provided.

  Stoer, J. and Bulirsch, R. (1993). Introduction to numerical analysis, vol. 12 of Texts in Applied Mathematics. Springer-Verlag, New York, second edn. Thomassin, I. (2008). Etude des tumeurs annexielles du pelvis feminin en IRM fonctionnelle, mise au point des techniques et applications cliniques. Physical PhD, Ecole doctorale Sciences et Technologies de l'Information des Télécomunication et des Systèmes . Wu, C. (1983). On the convergence property of the EM algorithm. Ann. Statist. 11, 95-103.The corresponding author for negotiations concerning the manuscript is Adeline Samson (adeline.samson@parisdescartes.fr)

Figure 1 :

 1 Figure 1: Noisy observations of one simulated data set (n = 200, ∆ = 0.2, σ 2 = 3) are plotted with stars. True simulated trajectories (thin solid lines), mean estimated trajectories (thick solid lines) and estimated 95% confidence intervals (dotted lines) obtained with the Kalman algorithm are plotted with dark lines for S(t), light lines for P (t) and very light line for I(t).

Table 1 :

 1 Mean estimated values (with empirical standard errors in bracket) obtained with the exact MLE and the EM algorithms and exact standard errors obtained from the Fisher information matrix, evaluated on 1000 simulated data with n = 200 and n = 1000 observations and σ 2 = 1 or σ 2 = 3 (σ 2 and θ 5 fixed to their true values).

	2

Table 2 :

 2 Mean estimated values (with empirical standard errors in bracket)

Acknowledgments

The authors thank V. Genon-Catalot for her constructive advices and help.

The authors thank C.A. Cuenod, D. Balvay and Y. Rozenholc for their helpful discussions on the biological problem. This project was supported by a grant BQR from University Paris Descartes managed by Y. Rozenholc.

the volume transfer constant. Inversely, the quantity of contrast agent exchanging from interstitium through plasma is equal to K trans times the concentration of contrast agent in interstitium. Lastly, the two-compartment model is:

For statistical accommodations, we use a new parameterization and set:

Model (20) can thus be transformed as follows:

B Proof of Proposition 1

The process X(t) = P -1 U (t) is solution of:

Applying Ito's formula, we obtain From this equation, we deduce:

where B(t, t+h) and Z(t, t+h) are given in Proposition 1. Using that W 1 , W 2 are independent and that X 0 is independent of (W 1 , W 2 ), we obtain the conditional

The stationary distribution can be deduced from equation (3) with a(t) = c. As the two elements of D are negative, we have

If X 0 ∼ N 2 (M, V ), an elementary computation shows that (X(t)) is strictly stationary.

C Link between the two parametrisations

New parameters (θ 1 , . . . , θ 5 ) are given as functions of initial parameters (β, λ, k, σ 1 , σ 2 ) in Section 3.2.1. Now we deduce (β, λ, k, σ 1 , σ 2 ) from (θ 1 , . . . , θ 5 ). From the definition of µ 1 and µ 2 (Section 2), we have

Thus we rewrite the matrix P and its inverse

and finally comes the matrix R R =

Notice that µ 1 = log(θ 1 )/∆ and µ 2 = log(θ 2 )/∆. The parameters β, λ, k, σ 2 1


Yazdan Arya

Guanghui Panah

Didier Fu

Bruno Dormont

Olivier Stankoff

Colliot

SCROLLING 2D U-NET FOR THE 3D SEGMENTATION OF THE INTERNAL CAROTID ARTERY

INTRODUCTION

There are several challenges attached with segmenting brain vasculature from angiographic images, including small vessels, uneven contrast, sparse structures, and unique anatomy. Precise and connected segmentations are necessary for analyzing cerebral vasculature and guiding the treatment of cerebral aneurysms. The main challenge is to cover all vessels branching from the main feeding artery in each image.

METHODS

Dataset & pre-processing

The dataset contained 35 rotational angiographies from subjects with intracranial aneurysms captured in Tours, France using a GE INNOVA 3D machine. The images have a resolution of 0.227x0.227x0.227mm 3 and include the cranium bone without subtraction. Radiopaque contrast dye highlights relevant vessel structures. Two experts in neurovascular anatomy and X-ray imaging delineated cerebral arteries branching from the Internal Carotid Artery (ICA). Images are re-scaled to [-1, 1], padded/cropped to (512,512,510), and reoriented to the canonical (RAS+) voxel orientation.

Model architecture

The segmentation model uses a modified 2D U-Net [1] with group normalization instead of batch-norm, "Leaky ReLu" instead of "ReLu", and up-sampling instead of up-convolutions. The model has 16 filters at the first convolution block output and performs three downsampling steps using average pooling. It takes 10 channels as input.

Training procedure

The dataset was divided into 22, 6, and 7 subjects for training, validation, and testing, respectively. The model was trained for 800 epochs using a batch size of 32, an Adam optimizer, and an initial learning rate of 0.001. The loss used is an average of the Dice loss and the binary cross Entropy (BCE). The training procedure is outlined under Algorithm 1. We have trained the model by sliding it along the inferior-superior direction (IS) and extracting a stack of N channels = 10 slices of the image.

Inference

Images are processed through the model in all N dirs = 6 anatomical directions: anterior-posterior, posterior-anterior, left-right, right-left, superior-inferior, and inferior-superior.

Algorithm 1 Training during 1 epoch, on IS direction X ← pad N (Batch image ) ▷ (B, 1, 512, 512, 510 + 2N ) Y ← pad N (Batch mask ) C ← (DiceLoss() + BCE())/2 for i = [N : 510 + N [ do x ← reshape(X[..., i : i + N ]) ▷ (B, N, 512, 512) y gth ← reshape(Y [..., i : i + N ]) y pred ← model(x) L ← C(y gth , y pred )/510
▷ gradient accumulation Backpropagate L end for

Step the optimizer All model outputs are summed and normalized by a factor N channels × N dirs . After visual inspection, we decided on a threshold of 0.3 to consider a voxel as correctly detected. Meaning it was detected at least through scrolling in one direction. The largest connected component (CC) is retained, resulting in output segmentations with only one CC.

RESULTS

Dice

CONCLUSIONS

A scrolling 2D U-Net is able to perform remarkably well in 3D segmentation tasks, with minor postprocessing cleanup.

[1] Olaf Ronneberger, Philipp Fischer, and Thomas Brox, "U-Net: Convolutional Networks for Biomedical Image Segmentation," arXiv:1505.04597 [cs], May 2015, arXiv: 1505.04597.

Table 1 .

 1 Performance of the proposed method.Table1summarizes results for our test dataset. Results are reported as mean ± standard error of the mean.

		Jaccard	MI	ClDice	AvgHausD
	0.94 ± 0.01 0.88 ± 0.01 0.02 ± 0.00 0.91 ± 0.02 0.82 ± 0.21
	Recall	Precision	VER	AVER	CC gth
	0.94 ± 0.01 0.93 ± 0.02 -0.01 ± 0.03 0.04 ± 0.02 1.71 ± 0.36


Arya Yazdan-Panah

Rosana El Jurdi

Guanghui Fu

Janan Arslan

Romain Valabregue

Didier Dormont

Bruno Stankoff

Olivier Colliot

Arya Yazdan Panah

Rosana El Jurdi

Scrolling 2D U-Net for the 3D segmentation of small vessels on 7T magnetic resonance angiography

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Advancements in 7T MRI have allowed for higher spatial image resolution and better signal-to-noise ratio for the visualization of small vessels in the brain. Given the current availability of data, most segmentation techniques in the field have focused on large blood vessel segmentation as opposed to small blood vessel extraction. Thus, there is a need to train effective small vessel segmentation tools with minimal data. In this paper, we propose a scrolling 2D U-Net that integrates information regarding the 3D volume while performing 2D segmentation. This approach was compared to the nnUnetthe state-of-the-art in segmentation models.

PROPOSED METHOD 2.1. Dataset

The challenge dataset consisted of 14 MRAs acquired at 7T MRI (dimension [480,640,163]). Annotations provided as part of the challenge were extracted using semiautomatic methods, such as the Frangi filter. The dataset was split into 8, 3, and, 3 for training, validation, and testing, respectively.

Training procedure

The training procedure is outlined under Algorithm 1. We have trained the model by sliding it along anatomical directions and extracting a stack of N channels = 10 slices of the image. The model was trained for 1000 epochs using a batch size of 8, an Adam optimizer, and an initial learning rate of 0.001. The average of the Dice loss and the binary cross Entropy (BCE) was used as the loss function. Images were randomly flipped along all dimensions with a probability of 0.5.

Inference

Images are processed through the model in all N dirs = 6 anatomical directions: AP, PA, LR, RL, IS, and SI. All outputs are summed and normalized by a factor N channels × N dirs . After visual inspection, we decided on a threshold of 0.3 to consider a voxel as detected (roughly, this means it was detected at least through scrolling in one direction).

Implementation details

The segmentation model uses a modified 2D U-Net [1] with group normalization, "Leaky ReLu", and up-sampling. The model has 16 filters at the first convolution block output and performs three downsampling steps using average pooling with 10 channels as input. Algorithm 1 Training during 1 epoch, on IS direction

X ← padN (Batchimage) ▷ (B, 1, 480, 640, 163 + 2N ) Y ← padN (Batch mask ) C ← (DiceLoss() + BCE())/2 for i = [N : 510 + N [ do x ← reshape(X[..., i : i + N ]) ▷ (B, N, 480, 640) y gth ← reshape(Y [..., i : i + N ]) y pred ← model(x) L ← C(y gth , y pred )/510
▷ gradient accumulation Backpropagate L end for

Step the optimizer

RESULTS

Table 1 summarizes performances obtained on the three testing subjects. Results are reported as mean ± standard error (SE). Obviously, mean and SE can be unreliable for such small samples but we were constrained by the challenge dataset. The proposed model performs comparably to the 3D nnUnet.

Method

Dice Jaccard MI AvgHaus AVER nnUnet 0.81 ± 0.04 0.69 ± 0.07 0.77 ±0.03 0.60 ± 0.15 0.17 ± 0.11 proposed 0.81 ± 0.05 0.68 ± 0.07 0.75 ± 0.04 0.64 ± 0.1 0.21 ± 0.08 Table 1. Summary of performance metrics 4. CONCLUSIONS The proposed method performs comparably to 3D state-ofthe-art methods like nnUnet [2]. The method has the advantage of being less computationally expensive owing to its 2D architecture.

[1] Olaf Ronneberger and et al, "U-net: Convolutional networks for biomedical image segmentation," in MICCAI, 2015, pp. 234-241.

[2] Fabian Isensee and et al, "nnu-net: a self-configuring method for deep learning-based biomedical image segmentation," Nature Methods, vol. 18, no. 2, pp. 203-211, Feb 2021.


M Verset

S Palierne

D Mathon

A Autefage

Comparison of the effect of locking versus standard screws on the mechanical properties of bone-plate constructs in a comminuted diaphyseal fracture model

Comparison of the effect of locking vs standard screws on the mechanical properties of bone-plate constructs in a comminuted diaphyseal fracture model

Michaël Verset, Sophie Palierne, Didier Mathon, André Autefage, Pascal

Swider

To cite this version:

Introduction

The purpose of this study was to compare the mechanical properties of bone-plate constructs with LCP plates used with either standard screws or with locking screws on an experimental model of comminuted fracture.

Materials and Methods

A model of comminuted fracture (mid-diaphyseal 5-mm gap) was created in ovine tibiae and treated with an 8-hole, 4.5-mm, broad stainless LCP plate with 3 bicortical self-tapping screws in each fragment.

In group 1 (n = 6), 4.5 mm cortical screws were placed in a slightly inner eccentric position and in group 2 (n = 6), 5.0-mm locking screws were used. Quasistatic nondestructive tests were performed in 4-point bending, torsion, and axial compression. A laser scan micrometer continuously measured the width of the osteotomy gap. Extrinsic stiffness of bone-plate constructs (K) was calculated, as well as stiffness loss (Kr) compared with the contralateral tibia. The slope of the loadgap displacement curve (S) was calculated. Results were statistically compared between the 2 groups by a 2-way ANOVA and Tukey's post hoc tests with significance set at P <.05.

Results

In 4-point bending, there were no significant differences between group 2 (G2) and group 1 (G1) for K (G2: 321.9 ± 65.2 kN/m, G1: 286.0 ± 36.6 kN/m), for Kr (G2: 55.4 ± 10.4%, G1: 64.5 ± 3.9%), and for S (G2: 292.0 ± 84.3 N/mm, G1: 271.4 ± 36.2 N/mm). In torsion, there were no significant differences between the 2 groups for K (G2: 122.1 ± 17.1 Nm/rad, G1: 107.3 ± 16.5 Nm/rad) and for Kr (G2: 44.0 ± 10.3%, G1: 55.0 ± 10.6%). In compression, group 2 showed significantly lower K (611.1±104.0 kN/m) and higherKr (66.5±10.3%) than group 1 (1019.2± 249.5 kN/m; 47.7 ± 10.8%). There was no significant difference between the 2 groups for S (G2: 254.6 ± 72.6 N/mm, G1: 438.2 ± 228.4 N/mm).

Conclusion:

Locking screws had no statistically significant effect on the mechanical properties of LCP-plated bones in 4-point bending and torsion compared to standard screws. In axial compression, locking screws induced a significant decrease in global mechanical properties without significant decrease in local mechanical properties of the osteotomy gap.

Michaël Verset, Sophie Palierne, Didier Mathon, André Autefage, Pascal Swider. Comparison of the effect of locking vs standard screws on the mechanical properties of bone-plate constructs in a comminuted diaphyseal fracture model. 21st ECVS Annual Scientific Meeting, Jul 2012, Barcelone, Spain. pp.E5, 10.1111/j.1532-950X.2012.00988.x. hal-04090448


R Rigo-Mariani

B Sareni
email: sareni@laplace.univ-tlse.fr

X Roboam
email: roboam@laplace.univ-tlse.fr

S Astier
email: astier@laplace.univ-tlse.fr

COMPARISON OF OPTIMIZATION STATEGIES FOR POWER DISPATCHING IN SMART MICROGRIDS WITH STORAGE

Keywords: optimal power dispatching, smart grid, trust region, dynamic programming, clearing genetic algorithm

The paper refers to the off-line control of a smart microgrid based on production and consumption forecast and energy costs of grid power flows. The system consists of an association of loads, PV arrays and a high speed flywheel. Several algorithms are implemented in order to estimate the optimal power dispatching for the day ahead that minimizes global energy costs.

INTRODUCTION

The growing number of decentralized sources based on renewable energy requires smarter operations to efficiently run power grids. Distribution networks evolve to a more meshed model, and they are likely to become associations of a large number of "microgrids" combining both consumption and production [START_REF] Celli | Meshed vs. radial MV distribution network in presence of large amount of DG[END_REF]. Adding a storage device allows more flexible operations and the so called "prosumers" are able to optimize the power exchanged with the main grid in compliance with the electricity market and the forecasts for the day ahead [START_REF] Colson | A Review of Challenges to Real-Time Power Management of Microgrids[END_REF]. The studied microgrid is presented in Fig. 1 and consists of an industrial load mixed with a 175 kW photovoltaic generator and a 100 kW/100 kWh flywheel. Section 1 describes the power flow model used to represent the system. Then, an off-line optimal power dispatching is introduced with the aim of estimating for the day ahead the grid power evolution that minimizes global energy costs. Several algorithms are developed and explained in this section in order to solve the optimal dispatching problem. Finally, simulations and results obtained with the investigated algorithms are illustrated in section 2.

1-POWER FLOW MODEL AND ALGORITHMS

Voltages and currents are not considered so far. The study only refers to the optimal dispatching of power flows taking account of the power converters efficiency. Equations are automatically generated using graphs theory and incidence matrices to describe the system [START_REF] Minoux | Graphes et algorithmes[END_REF]. In Fig. 1 (right part), the black arrows between nodes represent conventions for power flows while the possible directions (mono or bi directional ways) are illustrated by red arrows. For a given consumption and production respectively at nodes 2 and 9, three degrees of freedom are required to describe the whole system and to calculate all power flows related to each link while fulfilling the power balance on each node. To perform an optimization for the day ahead knowing the forecasts for P load and P PV , the dispatching procedures take account of the purchase electricity rate C purchase and of the sold energy rate C sale . Limitations on the grid power (P grid_max and P grid_min ) are also introduced if the shape of the power flowing thought the grid has to be limited. Four strategies are investigated for determining the power flows related to the three degrees of freedoms x = [x 5 (t), x 6 (t), x 9 (t)] that minimizes the daily energy cost C(x) while fulfilling load, grid and storage constraints. With a time step of 1 hour, all the vectors are 24 points long to represent a whole day. Thus, 72 parameters have be optimized in order to find a satisfying solution (i.e 24 points for each degree of freedom).

( 1)

) ( ) ( ) ( ) ( ) ( 11 24 0 1 t x t C t x t C C h t sale purchase       x .
A trust region method [START_REF] Coleman | An Interior Trust Region Approach for Nonlinear Minimization Subject to Bounds[END_REF] and the clearing genetic algorithm [START_REF] Sareni | Fitness sharing and niching methods revisited[END_REF] with constraint handling based on penalty functions are independently used for estimating the optimal references x 5 * (t), x 6 * (t), x 9 * (t) over the day ahead. Dynamic programming is also investigated through a standard scheme and an original self-adaptive algorithm so as to explore step by step the minimization of the cost. Both procedures exploit the Bellman-Ford algorithm for the determination of the shortest path between the initial and final point while discretizing the state of charge of the flywheel with a step size of 1kWh [START_REF] Bertsekas | Dynamic programming and optimal control, 2nd Edition[END_REF]. All investigated optimization strategies will be detailed in the final paper.

2-SIMULATIONS AND RESULTS

All procedures are tested and compared with regard to the quality of the objective function C(x) and the computational time on a standard PC computer with a clock speed of 2 GHz (see Table 1). The results below are obtained with C sale = 0, C purchase = 0.13€/kWh from 6 a.m. to 10 p.m. and C purchase = 0.08€/kWh otherwise. No limit on the grid power P grid is considered in this particular example. Table 1 shows the good efficiency of the trust region method that quickly finds an optimal solution with a non-feasible starting point arbitrary chosen. Good results were also obtained with the dynamic programming especially with the self-adaptive scheme that notably decreases the computational time. The clearing genetic algorithm was outperformed for the given hypothesis but it will be shown in the final paper that it can be relevant in other cases.

Figure 1 :

 1 Figure 1: Studied smart microgrid with PV production and flywheel storage

Fig 2

 2 Fig 2 shows some results obtained for a particular transient solution x that fulfills the constraints and for the optimal solution x * found by the trust region method. In Fig 2.b. the storage unit is used in a more efficient way for x * with full cycles from 0 to 100%. It allows the objective function to decrease from 40.07 € for x to 17.43 € for x * . As shown in Fig 2.a., the optimal solution limits the power flowing through the consumption meter when the energy cost is high.

Figure 2 :

 2 Figure 2: Illustration of the optimal dispatching. Comparison a particular solution (x) with the optimal solution (x * ) : grid power flow x 1 (a) and resulting energy storage E st (b)

Table . 1

 . Comparison of optimal dispatching strategies

	Dispatching strategy	Trust Region	Basic dynamic programming	Self-adaptive dynamic programming	Clearing genetic algorithm
	Cost C(x)	17.43 €	17.52 €	17.52 €	21.80 €
	CPU Time	2 min	90 min	19 min	75 min


Emmanuelle Comets
email: emmanuelle.comets@bichat.inserm.fr

Sarah Zohar

A survey of the way pharmacokinetics are reported in published phase I clinical trials, with an emphasis on oncology

Keywords: pharmacokinetics, phase I trials, reporting, survey

Objectives: During the drug development process, phase I trials are the first occasion to study the pharmacokinetics (PK) of a drug. They are performed in healthy volunteers, or patients in oncology, and are designed to determine a safe and acceptable dose for the later phases of clinical trials. We performed a bibliographic survey to investigate the way PK is described and reported in phase I clinical trials.

Methods:

We performed a MEDLINE search to retrieve the list of papers published between 2005 and 2006 and reporting phase I clinical trials with a PK study. We used a spreadsheet to record general information concerning the study, and specific information regarding the PK, such as the sampling times, number of subjects, and method of analysis.

Results:

The search yielded 349 papers, of which 37 were excluded for various reasons. Nearly all the papers in our review concerned cancer studies, although this was not a requirement in the search.

Consistent with the selection process, 84% papers explicitly stated PK as an objective of the study.

The methods section usually included a description of the PK (88%), but 10% of the papers provided no information concerning the methods used for the PK, and in 2% the description was only partial. The analysis method was usually basic, with non-compartmental or purely descriptive methods. Observed concentrations and area under the curves were the PK variables most often reported.

The results of the PK study were frequently reported in a separate paragraph of the results section, and only 22% of the studies related the PK findings to other results from the study, such as toxicity or efficacy. In addition, important information such as the number of subjects included in the PK study or the PK sampling scheme was sometimes not reported explicitly.

Conclusion:

Concerns about the decreasing cost-effectiveness in the drug development process

prompted the regulatory authorities to recently recommend a better integration of all available information, including in particular PK, in this process. In our review we found that this information was often either missing or incomplete, which hinders that objective. We suggest several improvements to the design and the reporting of methods and results for these studies, to ensure all relevant information has been included. PK findings should also be integrated in the broader perspective of drug development, through the study of their relationship with toxicity and/or efficacy, even in early phase I stages.

Introduction

In humans, clinical drug development is initiated with what is termed phase I trials. Phase I trials follow in vitro analyses and extensive animal studies designed to select the starting dose for use in humans, i.e., the maximum recommended starting dose (MRSD) [START_REF] Reigner | Estimating the starting dose for entry into humans: principles and practice[END_REF] . Phase I represents the first translation from basic laboratory work to the clinical setting. Usually, phase I trials are performed in healthy subjects, except when the drug is intended for the treatment of malignancies. Such studies are indeed characterized by the high potential toxicity of assessed drugs at any dose required to be effective and are therefore performed in patients, often patients who have failed several previous lines of treatments. Phase I trials are aimed at obtaining reliable information on a drug's safety, tolerability, pharmacokinetics (PK), and mechanism of action. More specifically, in a healthy volunteer study, the objective is to determine the maximum safe dose under a certain PK or pharmacodynamic (PD) safety limit; for a cancer study, the objective is to determine the maximum tolerated dose (MTD), defined as the highest dose with a relatively low risk of dose-limiting toxicity (DLT), or to recommend a dose level for phase II trials [START_REF] Zhou | Choice of designs and doses for early phase trials[END_REF] .

During phase I, sufficient information about the drug's PK and pharmacological effects should be obtained to permit the design of well-controlled, scientifically valid, phase II studies [START_REF] Aarons | Role of modelling and simulation in phase I drug development[END_REF] . Phase I studies can also be performed in later stages to investigate specific questions, such as the effect of food, modifications encountered in special populations, or drug-drug interactions [START_REF]Guideline for the format and content of the human pharmacokinetics and bioavailability section of an application[END_REF] . Phase I trials therefore often include a PK study. PK studies are intended to define the time-course of drug in the body, and are required in the registration files handed in to the regulatory authority. Furthermore, the determination of dose-exposure-effect relationships is now recognised to be a crucial part of the drug development process [START_REF] Dingemanse | Integrated pharmacokinetics and pharmacodynamics in drug development[END_REF] and is facilitated by the expanding development of biomarkers and analytical techniques. Modelling and simulation tools are already extensively and successfully used in many industrial fields outside the pharmaceutical industry, which has shown that PK/PD guided approaches can streamline the development process [START_REF] Gieschke | Pharmacometrics: modelling and simulation tools to improve decision making in clinical drug development[END_REF] . The Food and Drug Administration in the US devised a strategy, reported in March 2006 in "The Critical Path Opportunities Report", highlighting specific areas to improve and speeden the development of new, effective and safe medicines [START_REF]Innovation or stagnation: critical path opportunities report[END_REF] . A significant part of that report shows how modelling and simulation techniques can be used to incorporate all the information available from different stages of the clinical development to achieve this goal, using successful attempts at integrated PK/PD development as examples [START_REF] Chien | Pharmacokinetics/pharmacodynamics and the stages of drug development: Role of modeling and simulation[END_REF] .

Over the last decades, many methodological developments have been proposed to analyse data from PK studies. A large part of these relate to the so-called population approaches, where the average parameters and their associated variabilities are estimated through nonlinear mixed-effect models [START_REF] Pillai | Non-linear mixed effects modeling -from methodology and software development to driving implementation in drug development science[END_REF] , and which are invaluable when the design of the studies is sparse or variable between subjects [START_REF] Retout | Optimization of individual and population designs using Splus[END_REF] .

Non-compartmental approaches and individual nonlinear regression have also been extended. Many software have been developed to perform these different types of analyses. However, the extent to which these methods are actually applied in clinical trials has not been evaluated.

Our aim in this paper was to evaluate, through a bibliographic survey the way PK studies of phase I clinical trials, how these trials are reported and analysed. More specifically, we were particularly interested in the following points: (i) the description of the PK study; (ii) the methods used to analyse the data; (iii) the completeness of the results reported, (iv) the concordance between presented methodology and reported results, and (v) whether the results are used to bring additional information to the study. It should be noted that the aim of this paper is not to judge the quality of publications, since there is no validated tool for that, but to see if the PK part of a published clinical trials is well reported and can be understood by all readers.

Materials and methods

Article selection

PubMed was used to retrieve the references of all the papers corresponding to the following criteria: "phase I" or "phase 1" in title, pharmacokinetic* anywhere (title, abstract, keywords). The search was limited to 'English language' and to papers published in 2005 and 2006. The full text was then retrieved.

Data abstraction form and analysis

Our objective was to determine how much information concerning the PK study was reported in the papers describing the phase I study. We built a data abstraction form, which corresponds to a checklist of items [START_REF] Boutron | Methodological differences in clinical trials evaluating nonpharmacological and pharmacological treatments of hip and knee osteoarthritis[END_REF] . The form was then structured as a spreadsheet with each column representing one item, and the readers filled one line for each paper read. The form was used to extract a large number of informations from the articles and allowed statistical analysis of the results. In addition to items describing the pharmacokinetic study of the trials, we also extracted general information concerning the paper and the trial itself. We now describe in detail the items included in the data abstraction form.

We first included in the spreadsheet a number of items describing the authors of the study: the name and address of the first author, whether a statistician, a pharmacokineticist or a pharmacologist was present amongst the co-authors (based on the affiliations). We also recorded industrial partnership, by checking whether scientists from the pharmaceutical industry were co-authors, and whether the study was sponsored by the pharmaceutical industry, based on the declarations of conflict of interest and of grants received.

The next series of items described the phase I study: the title of the paper, the journal, the year of publication, the main drug tested and whether associated drugs have been tested, the primary and secondary objectives of the study, the patient population (adults or children) and the disease for which they were treated (if any). Since a large portion of phase I clinical trials deal with anticancer drugs, we also recorded whether the tumor response was evaluated when appropriate.

Given that phase I clinical trials are designed to investigate several doses of the drug and to test the tolerability of the drug, the next series of items described the design of the study: the number of subjects included, the planned number of drug levels, the dose-escalation scheme, the total number of dose-limiting toxicities that occurred during the study, whether MTD was reached or not, and whether patient selection was performed before the analysis. If several drugs were included, double-escalation was checked when both drugs under study were escalated, either together or sequentially.

Finally, a number of items specific to the PK analysis were then included in the spreadsheet: the number of subjects included in the PK study, whether the PK sampling scheme was described, whether the PK was sampled on one or more occasions, the dose levels for which PK was sampled, the PK variables reported, and whether the relationships PK-toxicity and/or PK-efficacy were investigated.

If several drugs were included, we also recorded whether a PK interaction study was performed. In addition, we recorded the analysis method used to analyse the PK, the software, and whether a model was built to describe the PK.

The first 40 papers were read by both authors who then met to compare the results entered in the data abstraction form and to assess the reproducibility between the two readers. All discordant results were checked and corrected. The next papers were divided equally and read by only one reader.

Once the survey was completed, we used descriptive statistics (median and range) for continuous variables and frequencies and percentages for categorical variables. Chi-square tests were used to perform stratified analyses with respect to the primary objective of the study. Data management and statistical analyses were performed using R (version 2.4.0) [START_REF]R: A language and environment for statistical computing[END_REF] .

Results

General results

The flowchart resuming the selection process is given in figure 1. The keyword selection retrieved 354 papers over the two year period considered; of these, 4 appeared twice in the list, and one study was split in two papers, which were regrouped for the analysis. Twenty-four papers were excluded because they did not meet the criteria we set on the search (date of publication, language and study performed in humans). Twelve papers were found not to be drug studies, or not to include phase I data. Finally, one paper was not available from any source we had access to, including direct contact with the author. In the following, we therefore used a database of 312 papers, of which 40 were read by the two authors and the remaining by only one author.

Table I gives the main characteristics of the studies. Industrial partnerships were identified based on co-authorship, financial support, or both in 86 papers (28%). The identification of a statistician or a pharmacokineticist was based solely on the affiliation which was sometimes limited to the university or the pharmaceutical company. As seen from the table, an overwhelming majority of papers published concerned cancer studies. The main objective of the study was exploratory in 139 papers (45%, of which 126 cancer studies), regrouping first-in-man or first-in children studies, new formulations, special population including use of drugs in new pathologies, while it was confirmatory in 173 papers (55%, of which 164 cancer studies), with a majority of these being studies evaluating com-bination therapies or new administration schedules. We also recorded the specific objectives stated in the paper: in accordance with the selection criteria, the determination of the PK was an objective in 84% of the studies (n=262). Because of the number of cancer studies found in the survey, the determination of the MTD was an objective in 215 studies (69%) and the investigation of toxicities or DLT in 165 papers (53%). The studies were nearly always designed with multiple objectives in mind (in N=279 studies, 89%). We also noted that in a few papers, the objective was not clearly stated in the abstract or the introduction. The number of subjects included in the studies was typically rather low, as expected from phase I studies, but a few studies included over 100 subjects.

Because we only found 22 papers which investigated non-cancer drugs, these tables and the following results will focus only on the 290 studies in the oncology field. We performed stratified analyses on the primary objective. The primary objective did not significantly change the proportion of papers involving a pharmacokineticist or a statistician, nor the proportion of studies performed in collaboration with an industrial. The size of the study however was slightly larger in confirmatory studies (median 32, range 6-105) than in exploratory studies (median 28, range 7-91, p<0.05).

Pharmacokinetic study

Our main focus in this paper was to examine how the PK part of phase I studies was reported. The results concerning the PK study are reported in table II and III. Of the 290 studies in the oncology field, 7 did not include PK, so the two tables report the findings on 283 papers.

In the methods section, the description of the PK section was available in most papers; however, 30 papers did not describe the design of the PK study at all or mentioned only the number of samples, and in 5 the description was insufficient (sometimes describing only partly the sampling schedule or mentioning weekly samples without the timing). Sampling was performed on at least two occasions in 65% of the studies, but in several papers, although the sampling schedule on the first occasion was described, we could not determine the number of additional samples or their exact timing when reading the paper. Also the number of subjects included in the PK study could not be found in 34 papers (12%) and was not apparent in many others.

Combination therapies were investigated in 120 papers (table III), of which 113 included 2 drugs (94%) and 7 included 3 drugs or more (6%). Combination studies represented 54% (n=51) of the studies performed by academia alone, but only 37% (n=69) of the studies performed with financial support from the pharmaceutical industry (p<0.007 according to a χ 2 test). Only 38 studies however included double escalation procedures, where the different drugs are escalated separately: in most cases only the main drug is escalated. Twenty-one of these studies were self-described as PK interaction studies, but in fact 48 studies investigated the interaction between the drugs, although half considered only the effect of one drug on the PK of the other(s).

PK data was most often analysed by non-compartmental methods (NCA), alone or in combination with more advanced approaches (n=14). In 21 papers (7%), the reporting was purely descriptive, showing plots of concentrations versus time or reporting summary statistics of trough or steady-state concentrations. The method of analysis was not reported at all in 50 papers. Modelling, using population approaches or individual estimation, was used in few papers (20% of the total); in one of the papers, a model described the PD of the drug but it was unclear whether the PK had also been modelled. In the 227 papers where the analysis was described and was not descriptive (n=212), the software used for the analysis was quite often not reported (n=37, 17%); when it was reported, the software WinNonLin, which can be used for NCA and RNL analyses, was the most frequently cited (n=119, 56%). As a result of this choice of methodology, the PK variables most often reported were the area under the curve (AUC) as a measure of drug exposure, observed concentrations such as con-centration profiles,peak (C max ), steady-state (C ss , or trough (C trough ) concentrations, and parameters obtained through non-compartmental analyses such as the clearance (CL), volume of distribution or half-life.

Extensive sampling was a rule: the median number of samples collected on a single occasion was 10 (range 1, for some studies collecting only steady-state samples in routine therapeutic drug monitoring, to 26), and the total number of samples varied from 2 to 54, with a median of 15; in 30 papers the number of samples was not clearly reported or varied within subjects as the samples were collected during long-term treatment. There was no difference in the number of samples depending on whether the method was NCA or RNL versus population methods. The number of subjects reported in the PK analysis was less than the total number of subjects included in the study in 138 papers.

Reasons for excluding data included "patient's will to participate in the PK study", or "concentrations undetectable for low doses", but they were seldom stated explicitly in the paper. PK was investigated usually at several drug levels and in majority at all drug levels studied (n=207 studies), although this information was often deduced from the tables reporting the PK results and not explicitly stated, and sometimes could not be determined (n=24).

Stratified analyses were performed to investigate the differences between papers using basic methods (descriptive or NCA) versus papers using modelling approaches. Papers published by academia alone were slightly more likely to involve modelling (non-linear regression or population methods, 22 out of 76) than papers involving an industrial partner (19 out of 157, p<0.003). When modelling approaches were used, there were also slightly less samples taken (median 33 versus 39 for studies using NCA), but this trend was not significant. The same results were obtained when assuming that the papers not reporting the analytical method (n=50) had used basic approaches.

The results of the PK study were usually described in a separate paragraph of the paper, and in majority no attempt was made to relate the PK results to the toxicity of the drug, or any other result from the study (efficacy, evolution of a biomarker...). In 34 papers an attempt was made at investigating the relationship between PK and toxicity: in 7 papers, this was purely descriptive, in 19 papers, a correlation or regression between a toxicity marker and PK parameters was tested, and in 8 papers a logistic model was used. In 35 papers the relationship between PK and efficacy was studied: in 7 papers this was descriptive, in 22 papers correlation or regression was used to relate the PK with the evolution of a biomarker, while in 6 papers logistic regression or modelling was applied. In 64 papers (22%) both efficacy and toxicity were studied in relationship with PK.

Discussion

In this paper we present a survey of the papers reporting phase I studies and including PK. The guidance for PK studies in humans, issued by the FDA, contains a number of recommendations concerning the content of the report to be filed when a new drug application is submitted [START_REF]Guideline for the format and content of the human pharmacokinetics and bioavailability section of an application[END_REF] . These items are considered necessary to judge the validity of an application, and although they cannot all be reported in a scientific article, we used this guidance as a basis to establish the data abstraction form described in the present paper.

We read 312 papers published over 2005 and 2006 and reporting phase I studies with PK. We found an overwhelming majority of the papers in our survey to be published in oncology. This reflects in part a selection bias in our initial Medline search. Indeed, our search included the following criteria;

(1) "Phase I" or "Phase 1" in the title of the published paper, (2) "Phamacokinetic" anywhere, (3) paper published between 2005 and 2006, (4) clinical trials as search limit, and (5) English paper as search limit. This Medline search yielded 349 papers in which 37 papers were excluded (review, not pharmacokinetic study, etc) and 290 (83%) were oncology clinical trials or related to oncology trials.

Phase I trials in oncology may be more likely to be identified as such in the title. To check this, a second Medline search was done to estimate the selection bias, in this new search all previous criteria were similar except the first one, which was changed to "Phase I" or "Phase 1" anywhere. 590 papers were found in which, based on title and abstract, 422 (71%) papers were identified as oncology or related to oncology trials, and 67 papers (11%) excluded. Thus, the high percentage of oncology trials found in the present study is slightly overestimated but oncology trials still represent the majority of published phase I trials.

The large majority of papers found in oncology with both the initial and modified search is probably a combination of two factors. The first is a selection bias, in that phase I trials appear more likely to be identified as such in the oncology field where there may be a more standardised approach to these studies. The second is a publication bias due to the way subjects are recruited. Indeed, subjects recruited in phase I studies in oncology are usually treated at the hospital and such studies are performed in collaboration with hospital physicians and researchers, whereas studies in healthy volunteers are mostly performed within the pharmaceutical industry or increasingly outsourced to clinical research centers, and publications are not necessarily encouraged. A possible work-around the publication bias could be to survey Final Study Reports, submitted to and reviewed by Health Authorities. However, these reports are not readily accessible, and existing registries such as www.clinicalstudyresults.org, developed by the pharmaceutical industry to provide greater access to the results of its clinical studies, or www.cancer.gov/clinicaltrials/ct-types-list, include mainly results from later phases in drug development.

An industrial partner was involved either as co-author or through financial support, in two-thirds of the papers surveyed, meaning one-third were purely academia. Support by the pharmaceutical in-dustry, in the form of a grant or by furnishing the drug under study, was reported in 40% of the papers surveyed; however, it was not always clear in the other 60% whether the drug had been furnished, so this figure may be slightly underestimated. Even so, we expected that most of phase I studies should involve an industrial partner, since these studies are usually initiated by the pharmaceutical industry. Again, this could be related to a publishing bias whereby phase I studies performed in the pharmaceutical industry are not all published.

We also attempted to determine the percentage of studies involving statisticians or pharmacokineticists/pharmacologists, but these items proved tricky to establish. Based on the affiliation reported alone, we found few departments of statistics or biostatistics explicitly stated, and similarly departments of pharmacology or pharmacokinetics were seldom found. It is also possible that, when the papers were published in collaboration with the pharmaceutical industry, the statistical and PK analyses were performed in-house. Although this does not necessarily mean that statisticians or pharmacokineticists were not part of the project, we did note that the PK analysis was usually quite basic.

In the remaining part of the paper, we focus on oncology trials owing to the limited number of papers outside this area.

The PK study was usually described in separate paragraphs, both in the methods and the results section, and rarely related to the other findings in the study such as toxicity. The description of the methods was usually appropriate, except in a few papers where, as pointed out, sampling times were not described. However the description of the results were not as standardised, and information was often missing. This made it more difficult to find the relevant results, especially for nonpharmacokineticists. PK results were also frequently reported separately for the different doses, with no attempt made at providing an overall information which could be extrapolated and used in later phases of drug development.

Modelling was seldom applied, and in majority PK data was analysed using non-compartmental approaches. When modelling was used, it did not appear to have an influence on the design of the studies. A large number of studies involve measurements on at least two occasions (such as after 2 doses), as per the recommendation of the regulatory authorities on population PK [START_REF]Guidance for Industry: Population Pharmacokinetics[END_REF] , but within subject variability was never investigated even though this may be a key point for future drug administration. Modelling approaches can be useful for decision purposes [START_REF] Aarons | Role of modelling and simulation in phase I drug development[END_REF] . PK or PK/PD models developed using the data from phase I studies can be applied to guide dosage decisions for phase II, evaluate alternative formulations or drug delivery systems, or investigate alternative dosage regimen for multiple dose studies [START_REF] Chien | Pharmacokinetics/pharmacodynamics and the stages of drug development: Role of modeling and simulation[END_REF] . As such, modelling can also play a pivotal role in the drug approval process [START_REF] Bhattaram | Impact of pharmacometrics on drug approval and labeling decisions: A survey of 42 new drug applications[END_REF] . Another point to note is that the number of subjects for which the PK was studied was often not reported, and in many papers the number of subjects with PK as reported in tables appeared less than the total number of subjects, while the reason for that was seldom reported. Sometimes the subjects were considered not be evaluable for PK and their data apparently discarded, but modelling approaches could have been used to incorporate all the available data even when NCA methods fail.

In a recent survey performed by the FDA concerning the low rate of success of new drug applications, the department of Pharmacometrics states that "one of the major criticisms against drug development is its negligence to employ prior knowledge to drive drug development decisions such as trial design and analysis", leading to the proposal to improve the process [START_REF]Innovation or stagnation: critical path opportunities report[END_REF] . Integrated PK/PD development has been proposed as a more effective use of knowledge and decision making to be used prospectively during drug development [START_REF] Bruno | Population pharmacokinetics and pharmacokineticpharmacodynamic relationships for docetaxel[END_REF][START_REF] Chien | Pharmacokinetics/pharmacodynamics and the stages of drug development: Role of modeling and simulation[END_REF] . In the present survey, we found that seldom was any attempt made to relate the PK findings to the other results of the study. In particular, hardly ever was the occurrence of toxicities related to individual PK parameters such as exposure, even though most of the studies were cancer studies where toxicities are dose-limiting and interindividual variability was usually high. It is of course possible that such a study was performed at another time during the development of the drug, or that it was reported in a separate paper or in the registration file submitted to the regulatory authorities.

In 2001, Margolin et al. performed a detailed review of clinical trials reporting high-dose chemotherapy regimens [START_REF] Margolin | Methodologic guidelines for the design of high-dose chemotherapy regimens[END_REF] . As we do in the present survey, they noted the lack of correlative pharmacologic analyses, and suggested guidelines for the design of such trials. Their suggestions included statistical input into the design, execution, analysis, interpretation, and reporting of these studies, as well as more homogeneity in the analysis methods [START_REF] Margolin | Methodologic guidelines for the design of high-dose chemotherapy regimens[END_REF] . Recent regulatory recommendations have in addition emphasised the benefits of better integrating all information during the drug development process, explicitly stating PK [START_REF]Innovation or stagnation: critical path opportunities report[END_REF] . Despite these recommendations however, our study shows that there is still room for improvement both in the reporting and in the analysis of the PK studies.

Several factors may explain the relative paucity of usable data we found in this survey. First, the amount of information contained in a phase I clinical trial requires a concise report of the different sections, including PK. However, some informations like the number of samples, the description of the sampling schedule, as well as the reasons why some data was discarded, need to be explicitly stated.

In our study, we often had to deduce such information from tables or figures, and sometimes we could not find it. Second, the emphasis of the publication may be more on the clinical findings than on the PK; however, in 84% of the papers we surveyed, the determination of the PK of the drug was explicitly stated as a primary objective and therefore probably deserves a more complete reporting. Third, the editorial process may also contribute to a standardisation of the reporting, however additional material can now often published online. Fourth, PK studies are also published in their own right and more sophisticated modelling may well have been performed after the reporting of the initial study, justifying a separate publication. However even in this case it is important to provide in the initial report a complete description of the data obtained, including the reasons for which some subjects are not included in the PK analysis.

Quality is a major issue in randomized clinical trials. For instance, the CONSORT (Consolidated Standards of Reporting Trials) statements were developed, encompassing various initiatives developed by the CONSORT Group to alleviate the problems arising from inadequate reporting of clinical trials. Although it can be argued that all randomized controlled trials are not comparable, the CON-SORT Statement developed a 22-item checklist to evaluate how the trial was designed, analyzed, and interpreted. In the same way a similar initiative could be realistic in the context of phase I clinical trials. Recently, Strevel et al. have proposed a quality score to assess the quality of abstract reporting for phase I cancer trials [START_REF] Strevel | Improving the quality of abstract reporting for phase I cancer trials[END_REF] . Similar quality scales or checklists could be developed to assess the reporting of the PK studies in published papers.

Recommendations

In order to improve the reporting of PK analyses in this field, we suggest improvements in three main areas: design, reporting of the methods used, and reporting of the results.

First, the design of the PK study should include the study of several, if possible all, doses of the drug; regulatory authorities are particularly interested in examining "the changes in kinetic parameters within the recommended dosing range", and suggest examining when appropriate the influences of demographic characteristics [START_REF]Guideline for the format and content of the human pharmacokinetics and bioavailability section of an application[END_REF] . Although the FDA guideline on human PK in new drug application did not explicitly mention assessing the within-subject variability when it was issued in 1987, a later guideline from the same agency and dealing with population pharmacokinetics devotes a subsection to the "importance of sampling individuals on more than one occasion" [START_REF]Guidance for Industry: Population Pharmacokinetics[END_REF] , which contributes to the changes in pharmacokinetic parameters. We therefore recommend to evaluate the within-subject variability in pharmacokinetic parameters, as well as the between-subject variability. In the present review, we found that most of the studies we collected administered multiple doses of the drug, which make this recommendation quite feasible, and a large number already followed it. Finally, around half of the papers we reviewed studied a combination of drugs; for these studies, sampling of all drugs should be performed to provide information on individual exposure for each drug.

Second, we recommend that the methods section describe: the number and timing of samples taken; the range of doses; the method used for the analysis of the PK, clearly stating whether modelling was used; the software and the estimation algorithm; if modelling was used, the method for data weighting and the model used for residual error; a motivated choice of the PK parameters of interest; whether PK was related to other findings in the study, and if so, a description of the analysis; for studies with several drugs, whether drug interaction was studied and the method used, as well as whether both drugs were escalated separately and the combined scalation scheme.

Third, the results section should contain the following information: the number of subjects included at each dose level; the number of subjects with PK and the justification for not including certain subjects in the pharmacokinetic analysis, the distribution of the estimated parameters at each dose level; the respective contribution of within and between subject variability for each parameter; a table summarising PK parameters across all doses if the PK was found to be dose-proportional, and PK parameters estimated using an appropriate model if it was not.

In addition to these three key areas, we would like to suggest trying to relate the PK findings to the other results in the study, such as toxicity and/or efficacy. Integration of PK results in the broader picture of drug development is currently lacking [START_REF]Innovation or stagnation: critical path opportunities report[END_REF] , as is a less discontinuous development process, and this has not much progressed since previous reviews [START_REF] Margolin | Methodologic guidelines for the design of high-dose chemotherapy regimens[END_REF] . For instance, when biomarkers of the toxicity or the efficacy are measured, modelling can be used to detect a concentration-effect relationship. The relationship of PK parameters of exposure such as the AUC with binary endpoints can be studied using logistic regression; boxplots of measures of exposure versus the endpoint can be plotted and descriptive statistics and correlation tests can be given. In all cases, the information could be tabulated and stored for use in later studies, possibly in Web appendices. Finally, phase I clinical trials often include extensive sampling to describe the entire PK profile, which could be more often modelled to provide a first description of the PK and guide subsequent developments.

Conclusion

In conclusion, the present survey showed that there is still progress to be made concerning the reporting and the analysis of the PK part of phase I clinical trials. The impact of the present survey is limited to clinical trials in oncology because of the small number of studies that were not related to cancer drugs. We propose a number of recommendations to improve the reporting of the PK sections of these trials.

  Figures

Table II :

 II Characteristics of the pharmacokinetic study reported in the 283 papers concerning oncology trials and including PK.

	Item		N (%)
	Description of PK	yes	248 (88%)
		partial	5 (2%)
		no	30 (10%)
	Multiple occasions	yes	183 (65%)
		no	86 (30%)
		Missing	14 (5%)
	Analysis method (*)	Descriptive	21 (7%)
		Non compartmental method (NCA) 171 (60%)
		Non linear regression (RNL)	47 (17%)
		Population approach (POP)	9 (3%)
		Not reported	50 (18%)
	Model built	yes	48 (17%)
		no	235 (83%)
	Relationship PK/toxicity yes	34 (12%)
		no	246 (87%)
		not applicable	3 (1%)
	Relationship PK/efficacy yes	35 (12%)
		no	241 (85%)
		not applicable	7 (2%)
	PK variables (*)	Observed concentrations	34 (12%)
		Cmax, Css, Ctrough,...	196 (69%)
		AUC	224 (79%)
		NCA parameters	173 (61%)
		CL	158 (56%)
		PK parameters	12 (4%)
		Other	42 (15%)
		none reported	13 (5%)
	(*) multiple answers possible

Tables