Thu Thuy Nguyen
email: thu-thuy.nguyen@inserm.fr

Caroline Bazzoli

France Mentré

Design evaluation and optimisation in crossover pharmacokinetic studies analysed by nonlinear mixed effects models

Keywords: crossover, design, Fisher information matrix, nonlinear mixed effects models, PFIM, power, sample size

teaching and research institutions in France or abroad, or from public or private research centers.

Design evaluation and optimisation in crossover pharmacokinetic studies analysed by nonlinear mixed effects models

Introduction

Pharmacokinetics (PK) is the science of the kinetics of drug in an organism, generally based on drug plasma concentrations. Bioequivalence trials are performed to compare the PK of two drug formulations and interaction trials are performed to study the effect on the PK of another drug. The most commonly used design for these trials is the crossover design, where measurements are made at several occasions. Bioequivalence is tested from the log ratio of the geometric means of two PK parameters: the area under the curve (AU C) of the concentrations and the maximal concentrations (C max ) [START_REF] Fda | Guidance for Industry -Statistical Aproaches to establishing bioequivalence[END_REF][START_REF] Emea | Note for guidance on the investigation of bioavailability and bioequivalence[END_REF]. Linear mixed effects models (LMEM) including treatment effect are usually used to analyse the log-transformed AU C and C max [START_REF] Hauschke | Bioequivalence studies in drug development[END_REF]. Bioequivalence tests are then performed on the estimates of the treatment effect. The same testing method is used for an assessment of the absence of interaction. The tested end-points AU C and C max are usually estimated by non-compartmental analysis (NCA) using the trapezoidal rule to evaluate AU C [START_REF] Gabrielson | Pharmacokinetic and pharmacodynamic data analysis: concepts and applications[END_REF]. NCA needs few hypotheses but has several limitations. Although using sparse sampling in NCA is possible [START_REF] Yeh | Estimation and significant tests of area under the curve derived from incomplete blood sampling[END_REF][START_REF] Jaki | A theoretical framework for estimation of aucs in complete and incomplete sampling designs[END_REF][START_REF] Jaki | Establishing bioequivalence in complete and incomplete data designs using aucs[END_REF], a large number of samples per subject (usually between 10 and 20, as in trials on healthy volunteers) is however needed when using this approach to be able to estimate mean and variability of AU C and C max with a good precision. NCA also does not take into account nonlinear pharmacokinetics, which can bias the bioavailability estimation [START_REF] Jusko | Nonlinear assessment of phenytoin bioavailability[END_REF] and may amplify small availability differences between drugs [START_REF] Hayashi | Estimation of rhG-CSF absorption kinetics after subcutaneous administration using a modified Wagner-Nelson method with a nonlinear elimination model[END_REF]. Alternatively, another approach to analyse PK data is the population approach based on nonlinear mixed effects models (NLMEM) [START_REF] Sheiner | Modelling of individual pharmacokinetics for computer-aided drug dosage[END_REF][START_REF] Sheiner | Pharmacokinetic/pharmacodynamic modelling in drug development[END_REF], which is more complex than NCA but benefits from the knowledge accumulated about the drug and can characterise the PK with few samples per subject, enabling us to analyse studies in patients with more power than does NCA. These models can also lead to a better understanding of the biological system and help to interpret ambiguous results. However, the use of NLMEM is rather recent and still rare in early phases of drug development and/or to analyse crossover studies in spite of its several advantages [START_REF] Panhard | Impact of modeling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with ap-plication to the interaction of tenofovir on atazanavir in HIV patients[END_REF][START_REF] Dubois | Bioequivalence tests based on individual estimates using non-compartmental or model-based analyses: evaluation of estimates of sample means and type I error for different designs[END_REF][START_REF] Dubois | Model-based analyses of bioequivalence crossover trials using the SAEM algorithm[END_REF].

The work presented here focusses on crossover PK trials analysed by NLMEM. Before the modelling step, data needs to be collected and we have consequently to define an appropriate design, which consists of determining a balance between the number of subjects and the number of samples per subject as well as the allocation of sampling times according to experimental conditions. The choice of design has an important impact on the study results, on the precision of the parameter estimates and on the power of the tests [START_REF] Kelman | Experimental design and efficient parameter estimation in population pharmacokinetics[END_REF][START_REF] Hashimoto | Designs for population pharmacodynamics: value of pharmacokinetic data and population analysis[END_REF][START_REF] Jonsson | Comparison of some practical sampling strategies for population pharmacokinetic studies[END_REF]. Indeed, a bad choice of design can lead to results which are difficult to interpret and minimise the interest of the study. The main approach for design evaluation has been for a long time based on simulations but it is a cumbersome method, and thus the number of designs which can be evaluated is limited. An alternative approach has been described in the general theory of optimum experimental design used for classical nonlinear models [START_REF] Atkinson | Optimum Experimental Designs[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF], relying on the inequality of Rao-Cramer which states that the inverse of the Fisher information matrix (M F ) is the lower bound of the variance-covariance matrix of any unbiased estimate of the parameters and its diagonal elements are the expected standard errors (SE) of the parameters. Several criteria based on M F have been developed to evaluate designs. One of the criteria widely used is the criterion of D-optimality, consisting of maximising the determinant of M F . But since there is no analytical form of the likelihood in our framework, the determination of the exact analytical expression of M F is not possible. An approximation of M F was first proposed by Mentré et al [START_REF] Mentré | Optimal design in random effect regression models[END_REF] for NLMEM. Its evaluation was implemented in the R function PFIM [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF] using first order linearisation of the model around the random effect expectation. Evaluation and optimisation of design have become a large field of research with more and more publications on several extensions of M F [START_REF] Hooker | Simultaneous population optimal design for pharmacokineticpharmacodynamic experiments[END_REF][START_REF] Gueorguieva | Optimal design for multivariate response pharmacokinetic models[END_REF][START_REF] Ogungbenro | Optimal design for multiresponse pharmacokinetic-pharmacodynamic models-dealing with unbalanced designs[END_REF][START_REF] Ogungbenro | Incorporating correlation in interindividual variability for the optimal design of multiresponse pharmacokinetic experiments[END_REF][START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: Evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF] as well as the creation of several software besides PFIM by different research teams, implemented in MATLAB (MathWorks Inc., Natick, MA, USA) (PopDes, PopED, WinPOPT). These softwares dedicated to evaluating and optimising designs were compared in an oral presentation by Mentré et al [START_REF] Mentré | Software for optimal design in population pharmacokinetics and pharmacodynamics: a comparison. 16th Population Approach Group in Europe[END_REF] at the meeting of Population Approach Group in Europe (www.page-meeting.org). However, these expressions of M F were not yet applicable to crossover trials with models including within subject variability (WSV) and discrete covariates changing between periods. It is important to model not only between subject variability (BSV) but also WSV when subjects have measurements at several periods, as accentuated by Karlsson and Sheiner [START_REF] Karlsson | The importance of modelling the inter-occasion variability in population pharmacokinetic studies[END_REF] and Panhard et al [START_REF] Panhard | Impact of modeling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with ap-plication to the interaction of tenofovir on atazanavir in HIV patients[END_REF]. The development of M F in NLMEM under uni-response models with additional terms for WSV has been implemented and illustrated on a specific example of two occasions [START_REF] Retout | Further developments of the Fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics[END_REF]. However, this has not yet been generalised to crossover designs for any number of periods and sequences, neither extended to include discrete covariates which can change between periods. The prediction of power for the Wald test to detect a treatment effect has been previously described by Retout et al [START_REF] Retout | Design in nonlinear mixed effects models: Optimisation using the Fedorov-Wynn algorithm and power of the wald test for binary covariates[END_REF] using the expected SE of the treatment effect but only for an interaction study in a parallel group design, not yet for an equivalence study in a crossover design. These developments would be very useful since to our knowledge, this is the first proposed statistical method to predict power of test and to choose optimal sampling times for crossover bioequivalence or absence of interaction trials analysed by NLMEM.

In this context, the first objective of this work is to extend the expression for M F in NLMEM including WSV in addition to BSV and discrete covariates fixed or changing between periods in crossover trials. Then, we also propose a computation of the expected power for the Wald test of comparison or equivalence and the number of subjects needed (NSN) for a given power. The relevance of these developments is evaluated by simulations mimicking a PK crossover two-period study of amoxicillin (an antibiotic widely used in bacterial infection treatment) in piglets. Finally, we apply these extensions to design a future PK crossover study in piglets, based on the results of a previous study inspiring our simulations, aiming to show the absence of interaction of a compound X on the PK of amoxicillin. We wish, through this motivating example, to show how our statistical methodology can be easily applied in practice, to choose optimal sparse sampling times with almost no loss of power and to compute the NSN for a given power. The example also brings us to discuss how different elements in a study design (number of subjects, of model parameters, of periods, etc.) could influence the power and to investigate/confirm the advantage of crossover design compared to parallel group design in bioequivalence trials.

We introduce some notation and present the new extension of M F as well as the computation of the power and NSN in Section 2. By simulation, we evaluate the relevance of these developments in Section 3. An application of this extension to a real example is presented in Section 4, followed by a discussion in Section 5.

Extension of population Fisher information matrix 2.1 Notations

Design

The elementary design ξ ih of individual i (i = 1, . . . , N ) at period h (h = 1, . . . , H) is defined by the number n ih of samples and their allocation in time (t ih1 , . . . , t ihn ih ). The global elementary design ξ i of individual i for H periods is defined as ξ i = (ξ i1 , . . . , ξ iH ), with a total number of samples n i = H h=1 n ih . Consequently, the population design for N individuals and for H periods can be defined as Ξ = {ξ 1 , . . . , ξ N }. Usually, population designs are composed of a limited number Q of groups of individuals with identical designs in each group. Each of these groups is composed of a global elementary design ξ q and is performed in a number N q of subjects. The population design can thus be written as

Ξ = {[ξ 1 , N 1 ]; . . . ; [ξ Q , N Q ]}.

Models

We denote by y ih the n ih -vector of observations for the individual i measured at period h and by f the known nonlinear function describing the nonlinear structural model. The NLMEM linking y ih to the sampling times ξ ih = (t ih1 , . . . , t ihn ih ) can be written as

y ih = f (φ ih , ξ ih ) + ǫ ih (1)
f maps φ ih , the P -vector of parameters for individual i at period h and the set ξ ih of n ih sampling times to a vector of which each element is a scalar f (φ ih , t ihj ), with j = 1, . . . , n ih . ǫ ih is the vector of random error which follows a normal distribution N (0, Σ ih ) with Σ ih a n ih × n ih diagonal matrix expressed as a function of (φ ih , σ inter , σ slope , ξ ih ) such as each diagonal element of Σ ih is (σ inter + σ slope f (φ ih , t ihj )) 2 , with j = 1, . . . , n ih . σ inter and σ slope are additional and proportional parameters in the variance model of the random error. The case of σ slope = 0 returns a homoscedastic error model, whereas the case of σ inter = 0 returns a constant coefficient of variation error model. The general case where the two parameters differ from 0 is called a combined error model. We denote by C the set of discrete covariates changing or not between periods, K c is the set of all possible categories k of the covariate c in C. We assume that the vector φ ih of parameters for an individual i at period h can be expressed as a function of the covariate c ih of individual i at period h, which can be additive (for normal parameters):

φ ih = µ + c∈C k∈Kc 1 c ih =k β c k + b i + κ ih (2)
or exponential (for lognormal parameters):

φ ih = µ ⊗ exp( c∈C k∈Kc 1 c ih =k β c k ) ⊗ exp(b i + κ ih ) (3)
Here, µ is the P -vector of fixed effects, β c k is the P -vector of fixed effects for category k of covariate c (β c k = 0 if k is the category of reference) on the P parameters of the structural model. b i is the P -vector of random effects of individual i, κ ih is the P -vector of random effects of individual i for period h. c ih can take any value k in the set K c . 1 c ih =k is the indicator function of the category k in the set K c , having the value 1 if c ih is in the category k, and having the value 0 if not. Consequently,

1 c ih =k β c k is equal to β c k if 1 c ih =k =
1 and to 0 if 1 c ih =k = 0. exp with a vector v as its argument is the vector such as each element is the exponential of the corresponding element of v. In (2), the P -vector φ ih is defined by addition of the vectors on the right side of the equation. In (3), the symbol ⊗ here means that each element p (p = 1, . . . , P ) of φ ih is obtained by multiplication of the p th elements (p = 1, . . . , P ) of the vectors µ and exp( c∈C k∈Kc 1 c ih =k β c k ) and exp(b i + κ ih ).

It is assumed that b i ∼ N (0, Ω) and κ ih ∼ N (0, Γ), with Ω and Γ defined as the diagonal variance-covariance matrices of size P × P . Each element ω p of Ω and γ p of Γ (p = 1, . . . , P ) represents, respectively, the variance of BSV for the p th component of b i and the variance of WSV for the p th component of κ ih . b i and κ ih are supposed to be independent. We denote by Ω * the diagonal matrix

Ω * = Ω Γ * (4)
where Γ * is a block diagonal matrix composed of H blocks Γ. The relation between φ ih and (µ, β, b i , κ ih ), with β regrouping all the fixed effects for discrete covariates, is summarised in a function g such as φ i = g(θ, ν i ) where φ i = (φ ′ i1 , ..., φ ′ iH ) ′ is the vector of parameters for individual i for H periods, θ = (µ ′ , β ′ ) ′ is the vector of all fixed effects, and

ν i = (b ′ i , κ ′ i1 , ..., κ ′ iH )
′ is the vector of all random effects for H periods. The vector of observations of individual i for H periods, y i = (y ′ i1 , ..., y ′ iH ) ′ can thus be written

y i = f (g(θ, ν i ), ξ i ) + ǫ i (5)
where

ξ i = (ξ i1 , ..., ξ iH ) and ǫ i = (ǫ ′ i1 , ..., ǫ ′ iH ) ′ ∼ N (0, Σ i (θ, ν i , σ inter , σ slope , ξ i )). Σ i is an n i × n i diagonal matrix composed of each diagonal element of Σ ih (h = 1, . . . , H).
The following assumptions are made: ǫ i |ν i are independent between subjects, ǫ ih |κ ih are independent between periods, and ǫ i and ν i are independent for each subject.

Let λ = (ω 1 , ..., ω P , γ 1 , ..., γ P , σ inter , σ slope ) ′ be the vector of variance terms, and Ψ be the vector of all population parameters to be estimated, so that Ψ = (θ ′ , λ ′ ) ′ . We also denote by dim(θ) the number of fixed effects and dim(λ) the number of variance terms.

Population Fisher information matrix

We extend the expression of M F in NLMEM for crossover trials using the first order Taylor expansion of the model as in Mentré et al [START_REF] Mentré | Optimal design in random effect regression models[END_REF]. For simplicity, we omit the index i for the individual in this section. The elementary M F (Ψ, ξ) for an individual with elementary design ξ composed of one or several periods is calculated by

M F (Ψ, ξ) = E -∂ 2 l(Ψ, y) ∂Ψ∂Ψ ′ (6)
where l(Ψ, y) is the log-likelihood of the vector of observations y of that individual for the population parameters Ψ. There is no analytical expression for l(Ψ, y) because of the nonlinearity of the structural model f . An approximation is possible using the first-order Taylor expansion of f (g(θ, ν), ξ) around the expectation of the random effects ν, i.e. 0. The statistical model can be written as

y ∼ = f (g(θ, 0), ξ) + ∂f ′ (g(θ, 0), ξ) ∂ν ν + ǫ (7)
and the log-likelihood l is then approximated by

-2l(Ψ, y) ∼ = nln2π + ln|V | + (y -E) ′ V -1 (y -E) ( 8
)
where n is the total number of observations of the individual; E and V are the approximated marginal expectation and variance of y given by

E = E(y) ∼ = f (g(θ, 0), ξ) (9)
V = V ar(y) ∼ = ∂f ′ (g(θ, 0), ξ) ∂ν Ω * ∂f (g(θ, 0), ξ) ∂ν ′ + Σ(θ, 0, σ inter , σ slope , ξ) (10)
We derive the expression of the log-likelihood l, taking into account WSV in addition to BSV and fixed effects for discretes covariates changing between periods. We obtain the elementary M F depending on E and V . Assuming that V does not or does depend on the fixed effects, the elementary M F is a block diagonal matrix or a full matrix. There is no clear consensus on what is the best approximation but here in our approach by linearisation, we assume the choice of the block diagonal expression [START_REF] Mielke | Some considerations on the Fisher information in nonlinear mixed effects models[END_REF]:

M F (Ψ, ξ) = 1 2 A(E, V ) 0 0 B(E, V ) (11)
where

(A(E, V )) rs = 2 ∂E ′ ∂θ r V -1 ∂E ′ ∂θ s ( 12
) (B(E, V )) ml = tr ∂V ∂λ m V -1 ∂V ∂λ l V -1 (13)
with r and s = 1, ..., dim(θ); m and l = 1, ..., dim(λ).

To define the population matrix, we use the fact that M F (Ψ, Ξ) for a population design Ξ is defined by the sum of the N elementary matrices M F (Ψ, ξ i ), so that

M F (Ψ, Ξ) = N i=1 M F (Ψ, ξ i ) (14)
In the case of a limited number Q of elementary designs, we have

M F (Ψ, Ξ) = Q q=1 N q M F (Ψ, ξ q ). (15)

Computation of power and number of subjects needed

In this section, we consider the Wald tests of comparison and equivalence on a discrete covariate effect fixed or changing between periods for the case of a crossover trial. Here, these tests are performed on the estimator of the effect β c k p of each category k of each covariate c for parameter p, noted β c k p . Given a value β 1 of β c k p , and from the square root of the diagonal terms of the inverse of M F , we calculate the standard error of β 1 . Then, we predict the power of the tests and also compute the NSN to achieve a given power. Note that for simplicity, we omit index c k p for β in this section.

Comparison test

The null hypothesis for the Wald test of comparison is H 0 : {β = 0} while the alternative hypothesis is H 1 : {β = 0}. Putting β as the estimator of β, then the statistic of the Wald test under the null hypothesis

H 0 is W = β/SE( β). With a type I error α, H 0 is rejected if |W | > z 1-α/2
, where z 1-α/2 is the (1α/2) quantile of the standard normal distribution. We compute the power of the Wald test of comparison under H 1 , when β = β 1 = 0. With a given design and values of population parameters, we predict SE(β 1 ) using the extension of M F . Under H 1 , W is asymptotically distributed with a normal distribution centred at β 1 /SE(β 1 ). Therefore, the power of the Wald test is

P comp = 1 -Φ z 1-α/2 - β 1 SE(β 1 ) + Φ -z 1-α/2 - β 1 SE(β 1 ) ( 16
)
where Φ is the cumulative distribution function of the standard normal distribution.

Using the predicted SE by the extension of M F , we can also derive the NSN to achieve a power P comp to detect a covariate effect as in [START_REF] Retout | Design in nonlinear mixed effects models: Optimisation using the Fedorov-Wynn algorithm and power of the wald test for binary covariates[END_REF]. For that, we first compute the SE needed on β 1 to obtain the power of P comp , called SEN(P comp ), using the following relation

SEN(P comp ) = β 1 z 1-α/2 -Φ -1 (1 -P comp ) (17)
We then compute the NSN to obtain a power of P comp , called NSN(P comp ) using

NSN(P comp ) = N × SE(β 1 ) SEN(P comp ) 2 (18)

Equivalence test

The null hypothesis for the Wald test of equivalence is H 0 : {β ≤ -δ or β ≥ +δ} while the alternative hypothesis is H 1 : {-δ < β < +δ}, where δ is the equivalence limit. In the guidelines for equivalence assessment [START_REF] Fda | Guidance for Industry -Statistical Aproaches to establishing bioequivalence[END_REF][START_REF] Emea | Note for guidance on the investigation of bioavailability and bioequivalence[END_REF], it is recommended that δ = 0.2. H 0 is composed of two unilateral hypotheses H 0,-δ : {β ≤ -δ} and H 0,+δ : {β ≥ +δ}. Equivalence between two groups on covariate effect β can be concluded if and only if the two hypotheses H 0,-δ and H 0,+δ are rejected. The equivalence Wald test using NLMEM has already been developed [START_REF] Panhard | Evaluation by simulation of tests based on non-linear mixed effects models in pharmacokinetic interaction and bioequivalence cross-over trials[END_REF][START_REF] Panhard | Impact of modeling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with ap-plication to the interaction of tenofovir on atazanavir in HIV patients[END_REF], based on the Schuirmann's two one-sided tests (TOST) [START_REF] Schuirmann | A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability[END_REF]. Denoting the estimator of β by β, the following formulas define W -δ , the statistic of the unilateral test under the null hypothesis H 0,-δ , and W +δ , the statistic of the unilateral test under the null hypothesis H 0,+δ :

W -δ = ( β + δ)/SE( β) and W +δ = ( β -δ)/SE( β). With a type I error α, H 0 is rejected if W -δ ≥ z 1-α and W +δ ≤ -z 1-α , where z 1-α is the (1 -α) quantile of the standard normal distribution.
We compute the power of the Wald test of equivalence under H 1 , when β = β 1 ∈ [-δ, +δ]. Usually we choose β 1 = 0. From a given design and values of population parameters, we predict SE(β 1 ) using the extension of M F . The power of the Wald test is then

P equi = 1 -Φ z 1-α - β 1 + δ SE(β 1 ) if β 1 ∈ [-δ, 0] (19)
P equi = Φ -z 1-α - β 1 -δ SE(β 1 ) if β 1 ∈ [0, +δ] (20)
where Φ is defined as in Section 2.3.1. As expressed in [START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF] and [START_REF] Mentré | Optimal design in random effect regression models[END_REF], the power of the equivalence test depends on the sign of β 1 . When β 1 = 0, we can use any of the two equations to obtain the power because in that case:

P equi = 1 -Φ z 1-α - δ SE(β 1 ) = Φ -z 1-α - (-δ) SE(β 1 ) (21)
Using the predicted SE by the extension of M F , we can also derive the NSN to achieve the power P equi of showing the equivalence between two groups on covariate effect β 1 . For that, we first compute the SE needed on β 1 to obtain the power of P equi , called SEN(P equi ), using

SEN(P equi ) = -β 1 -δ -z 1-α + Φ -1 (1 -P equi ) if β 1 ∈ [-δ, 0] (22)
SEN(P equi ) = -β 1 + δ z 1-α + Φ -1 (P equi ) if β 1 ∈ [0, +δ] (23)
We then compute the NSN to obtain the power of P equi , called NSN(P equi ) using ( 18) with SEN(P equi ) instead of SEN(P comp ). If β 1 < 0, we substitute ( 22) into (18); if β 1 > 0, we substitute ( 23) into [START_REF] Atkinson | Optimum Experimental Designs[END_REF]; if β 1 = 0, we can substitute either of ( 22) or ( 23) into [START_REF] Atkinson | Optimum Experimental Designs[END_REF] because in that case:

SEN(P equi ) = -δ -z 1-α + Φ -1 (1 -P equi ) = δ z 1-α + Φ -1 (P equi ) (24)

Evaluation by simulations

The simulation model and design are based on a previous PK study DAV1 of amoxicillin in piglets, conducted by the biotechnology company Da Volterra (http://www.davolterra.com) which will be given in detail in Section 4.1.

Simulations

The PK model is a one compartment oral model with first order absorption and first order linear elimination. The observed amoxicillin concentration of subject i at sampling time j of period h is modelled by a NLMEM as described in Section 2.1.2:

y ihj = f (φ ih , t ihj ) + ǫ ihj = F D ka ih Cl ih -V ih ka ih [e -ka ih t ihj -e -(Cl ih /V ih ) t ihj ] + ǫ ihj (25)
In [START_REF] Ogungbenro | Incorporating correlation in interindividual variability for the optimal design of multiresponse pharmacokinetic experiments[END_REF], D is the dose, F the bioavailability, ka the absorption rate constant, Cl the clearance of elimination of the drug and V the volume of distribution. After oral administration only, the bioavailability can not be estimated and the vector φ of PK parameters is (ka, V /F, Cl/F ). For simplicity, this will be denoted subsequently by (ka, V, Cl) where Cl and V are apparent parameters.

We simulate crossover PK trials with two periods, one sequence as in DAV1 study. In each dataset, N = 40 subjects are allocated at period 1 to treatment 1, with drug A plus placebo and then, at period 2, to treatment 2, with the same drug A plus drug B. All compounds are administrated orally once a day. The dose is fixed at 30 mg.kg -1 for all subjects. The sampling times for all subjects and both periods are similar and are taken from among the times of the DAV1 study: 0.5, 1, 1.5, 2, 4, 6, and 8 h after dosing. The vector of fixed effects µ is composed of (µ ka = 1.00 h -1 , µ V = 3.5 L.kg -1 , µ Cl = 2 L.h -1 .kg -1 ) for treatment 1 (the treatment of reference). We take into account a treatment effect β Cl on Cl, characterising the interaction between A and B. The relation between these parameters is expressed as in [START_REF] Hauschke | Bioequivalence studies in drug development[END_REF], so that at period 1,

φ i1 =   µ ka × exp(b ka i + κ ka i1 ) µ V × exp(b V i + κ V i1 ) µ Cl × exp(b Cl i + κ Cl i1 )   (26)
and at period 2,

φ i2 =   µ ka × exp(b ka i + κ ka i2 ) µ V × exp(b V i + κ V i2 ) µ Cl × exp(β Cl ) × exp(b Cl i + κ Cl i2 )   (27)
We fix the variance of BSV ω = (0.3) 2 = 0.09 and of WSV γ = (0.15) 2 = 0.0225 for the three parameters. We assume no correlation between PK parameters, so the variance-covariance matrices Ω and Γ are diagonal. The parameters of the random error model are σ inter = 0.1 ng.L -1 and σ slope = 0. For each subject i of each simulated trial at period h, we simulate a vector of random effects b i in N (0, Ω) and two vectors of random effects κ ih in N (0, Γ), one for each period and then compute φ ih while taking into account the treatment effect β Cl using [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: Evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF][START_REF] Mentré | Software for optimal design in population pharmacokinetics and pharmacodynamics: a comparison. 16th Population Approach Group in Europe[END_REF]. We calculate the concentrations predicted by the PK model at each sampling time. For each predicted concentration, we generate the vector of the random error using the value of σ inter . To get the simulated concentrations, these errors are added to the predicted concentrations. Then, we leftcensored them at 0.05 ng.L -1 , the limit of quantification (LOQ) in the DAV1 study. We consider in the above simulation procedure different values of the treatment effect parameter β Cl such as log(0.8), log(1), log(1.1), log(1.2), log(1.25) and log(1.5). The case β Cl = log(1) = 0 corresponds to the null hypothesis of the comparison test (no treatment effect) while the cases β Cl = log(0.8) ≈ -0.2 or log(1.25) ≈ 0.2 correspond to the null hypothesis of the equivalence test (the two limits of equivalence). We also derive and study an optimal sparse design of four samples by optimising the rich one using the new extension of M F with the Fedorov-Wynn algorithm [START_REF] Fedorov | Theory of Optimal Experiments[END_REF][START_REF] Wynn | Results in the theory and construction of D-optimum designs[END_REF]. For that, we use the same simulated values above for a model without treatment effect then with different values of treatment effect parameter. The optimal sparse design obtained is the same with or without treatment effect (0.5, 2, 6, 8 h after dosing). For each simulated scenario, we simulate 1000 trials, using R 2.4.0.

Evaluation methods

Our aim is to evaluate and compare the predictions of the standard errors and test powers computed by the extension of M F to the empirical ones obtained from simulations. For that, first, to obtain these predictions, with the design and the population parameter values of each simulated scenario above, we compute the standard errors denoted by SE of each population parameter from the square root of the diagonal terms of the extended M -1 F or the relative SE denoted by RSE, defined as the ratio of SE by the parameter simulated value. Then, from the SE of treatment effect parameter, we also predict the power of the Wald test of comparison and equivalence. We consider a type I error α = 0.05 and an equivalence limit δ = 0.2.

In parallel, the empirical values are calculated as follows. Each simulated data file of each scenario with different values of β Cl is analysed in MONOLIX 2.4 (http://www.monolix.org) by an NLMEM including discrete covariate. We use the SAEM algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of EM algorithm[END_REF][START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects model[END_REF][START_REF] Samson | The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed effects model[END_REF][START_REF] Panhard | Extension of the SAEM algorithm for nonlinear mixed effects models with two levels of random effects[END_REF] implemented in this software to estimate parameters. The data under LOQ are taken into account in the analysis using the extension of the SAEM algorithm proposed by Samson et al [START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: application to HIV dynamics model[END_REF], which allows of handling the left-censored data in NLMEM as an exact Maximum Likelihood estimation method. We then calculate the empirical SE, defined as the sample estimate of the standard deviation from the parameter estimates obtained with SAEM. Concerning the powers, the observed power of the Wald test of comparison and equivalence are calculated as the proportion of trials for which the null hypothesis is rejected. Finally, we compare the results obtained from simulations to the predictions.

We present in Section 3.3.1 the comparison of the predicted RSE to the empirical ones in each simulated scenario. As the treatment effect is added on Cl, here are presented only the results for the fixed effect, variances of BSV and WSV on Cl, as well as the for the parameter σ inter of the error model. Then, we focus on the comparison of the SE of treatment effect β Cl . The relevance of this power computation using the extension of M F will be investigated in Section 3.3.2.

Results

Evaluation of the predicted standard errors

The RSE predicted by the extension of M F and the empirical ones calculated from simulations are displayed in Figure 1. The predicted RSE are very close to the empirical ones in each simulated scenario. Some RSE are slightly underestimated by M F , especially for γ Cl , but the difference is not clinically relevant (< 10%). As expected, the RSE are slightly higher in the sparse design than in the rich design, but the differences are very small, except for σ inter .

We now focus on the SE of treatment effect β Cl . The distributions of the observed SE from each simulated scenario are reported as boxplots in Figure 2. The SE(β Cl ) predicted by the extension of M F are very close to the empirical ones. The SE(β Cl ) are slightly higher in the rich design compared to the sparse one but they remain very close. These values are reported in Table I. There is no important change of SE(β Cl ) when varying β Cl .

Although not all the predicted SE and RSE are reported here, some of these values are present in the online supplementary evaluation output of the rich design1 with β Cl = log(1.1). In this example and for other scenarios including the ones with optimal designs as well, the RSE predicted for the fixed effects are very close to each other and close to 5%. The minimum values of RSE of the log parameters can be approximated by the SE, in the case of a rich design, as √ ω/ √ N (ω being here the variance of BSV). Here, √ ω/ √ N are roughly 4.9% since we have fixed the same ω for the three PK parameters (ka, V, Cl). This means that the original design (seven samples) is quite rich and the derived optimal design (four samples) is a good design, allowing of estimating well the parameters with good precision.

Evaluation of the predicted powers

With each scenario, the predicted power by the extension of M F and the observed power from simulations for the comparison test as well as for the equivalence test are reported in Table II and also as barplots in Figure 3. Note that when exp(β Cl ) = 1 for the comparison test and when exp(β Cl ) = 0.8 or exp(β Cl ) = 1.25 for the equivalence test, we are under the null hypothesis so the probability computed is not the test power but the type I error. The predicted powers by the extension of M F and the observed powers from simulations are very close for both comparison and equivalence tests. The powers obtained with the optimal sparse design are also close to the ones obtained with the rich design.

Application to a pharmacokinetic study in piglets

In this section, we apply the extension of M F to design a subsequent study DAV2 of amoxicillin in piglets based on the results of the previous study DAV1 conducted by the company Da Volterra that inspired the simulations. DAV2 will be a one-way crossover trial. The piglets receive at period 1 amoxicillin plus placebo, at period 2 amoxicillin plus a compound X developed by Da Volterra. The objective of DAV2 is to show the absence of interaction of X on the AUC of amoxicillin concentrations. The direct relation AU C = dose/Cl implies the equality between the log parameter of the treatment effect on AU C and the one on Cl [START_REF] Panhard | Impact of modeling intra-subject variability on tests based on non-linear mixed-effects models in cross-over pharmacokinetic trials with ap-plication to the interaction of tenofovir on atazanavir in HIV patients[END_REF][START_REF] Dubois | Bioequivalence tests based on individual estimates using non-compartmental or model-based analyses: evaluation of estimates of sample means and type I error for different designs[END_REF]. That is why we are interested particularly in Cl. To design DAV2, we analyse the data of DAV1 then with the results obtained, we use the extension of M F to evaluate/optimise designs and to compute the NSN to show the absence of interaction of X on Cl of amoxicillin with a power of 90%.

Analysis of DAV1 study

Materials and methods

The DAV1 study is a one-way, two-period crossover trial with 16 piglets receiving amoxicillin and placebo at period 1 and only amoxicillin at period 2. Amoxicillin was administrated orally with a single dose of 30 mg.kg -1 . Plasma amoxicillin concentrations were collected at same sampling times at each period: 0.5, 1, 1.5, 2, 4, 6, 8, 10 and 12 h. Figure 4 represents the individual amoxicillin concentrations (ng.L -1 ) versus time of each period.

To describe the pharmacokinetics of amoxicillin, we analyse jointly the data of both periods for all piglets using MONOLIX 2.4 with an NLMEM and the SAEM algorithm for estimating parameters. At 8, 10 and 12 h, respectively, 21.9%, 68.8% and 87.55% of observations of both periods are below LOQ (0.05 ng.L -1 ). The very high proportions of data below LOQ especially at 10 h and 12 h justify the choice of excluding these two samples in the simulations previously. Here, MONOLIX 2.4 allows of taking into account the data below LOQ [START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: application to HIV dynamics model[END_REF]. A one compartment oral model with first order absorption, first order linear elimination and the presence of a lag time adequately describes the amoxicillin concentrations in this study. The fixed effects are then (tlag, ka, CL/F, V /F ) but denoted for simplicity by (tlag, ka, Cl, V ). A treatment effect β Cl is added on log(Cl) and is estimated. We consider log parameters and assume the choice of diagonal matrices Ω and Γ. To choose the error model, we compare the Bayesian Information Criteria (BIC) between the homoscedastic (σ slope = 0), proportional (σ inter = 0), and combined models. We examine the SAEM convergence graph and the goodness-of-fit plots to evaluate the chosen model.

Results

We choose a homoscedastic error model with σ slope = 0. With respect to the model of random effects, when estimating ω tlag and ω Cl , the obtained variabilities are very low and their RSE very high. Consequently, we fix ω tlag = ω Cl = 0 for Ω and do not estimate these two variances. With the chosen model, the algorithm converges and the goodness-of-fit plots are satisfactory. The population parameter estimates with this model are reported in Table III. Because the dose are in mg.kg -1 , the unit of the parameters Cl and V are also expressed by their usual units per kg. The estimation of most of these parameters is good because the RSE are satisfactory except for ω ka and γ ka of which the RSE are, respectively, 144% and 121%. We note that these variabilities are low and are consequently difficult to be precisely estimated. With β Cl = 0.06 and SE( β Cl ) = 0.16, we perform, using the formulas in Section 2.3, the Wald tests of equivalence with the equivalence limit δ = 0.2 and the type I error α = 0.05. The 90% confidence interval for treatment effect on Cl is [-0.198; 0.318] and not included in [-0.2, 0.2]. Consequently, we can not show a significant equivalence between the two treatment groups.

Designing of DAV2 study 4.2.1 Methods

We evaluate and optimise the design for the DAV2 study, using the model selected and the values of population parameters estimated from the DAV1 study, with the estimated treatment effect = 0.06 (cf. Table III). We first evaluate for the next study a rich design with two periods which is the original design of the DAV1 study (0.5, 1, 1.5, 2, 4, 6, 8 h after dosing) with same sampling times at each period. We do not consider the samples 10 h and 12 h because of the high proportions of data below LOQ as there is no way presently to incorporate this information in the computation of M F . Consequently, this population design involves a total number of 224 samples to be performed in 16 piglets. We omit subsequently the parameter tlag because its estimated value is small (roughly 0.37 h) whereas the first sampling time of the study is given to be after 0.5 h. We also derive a two-period optimal design with a constraint of four points for each period (same samples at each period) from the rich one using the extension of M F and the Fedorov-Wynn algorithm as in the simulations. Thus, the optimal design involves a total number of 128 samples to be performed in 16 piglets. We evaluate the values of SE and RSE for all parameters and examine the critetion value computed from the extension of M F , defined as the determinant of M F normalised by the inverse of the total number of parameters to be estimated with NLMEM (the size of M F ). The ratio of the two criteria expresses the relative efficiency between the two designs and evaluates the factor of the mean estimation variance decrease when changing the design.

Then, we compute, with the extension of M F , the predicted power of the Wald test assessing the equivalence on the parameter Cl between the two groups: amoxicillin + placebo (AP) group (treatment 1) versus amoxicillin + compound X (AX) group (treatment 2) from the expected SE of β Cl and with the same number of piglets (N = 16) as in DAV1. We examine different given values of β Cl : 0 as is usually done in equivalence assessment, 0.06 as in DAV1 and 0.1, for the rich design as well as for the sparse design. We predict the NSN to achieve a power of 90% for the equivalence test with an equivalence limit of 0.2. The computations are first performed for a design with two periods AP/AX then with four periods AP/AX/AP/AX.

Results

The evaluation results for the two-period design of DAV2 study are reported in Table IV with the predictions of the RSE (or SE for β Cl ) for each parameter and the criterion value of both the original and optimal design composed of (0.5, 2, 4, 6 h) at each period. The original design gives correct RSE and SE for all parameters except for ω ka and γ ka as in DAV1. We note also that the RSE and SE of the sparse design are slightly higher than those given by the rich design but they are still very close. The criterion value decreases slightly when there are less sampling times per subject. The efficiency of the sparse design compared to the rich one is calculated as the ratio of the two criterion values 32.71/39.04 = 0.82. Concequently, the choice of similar samples as DAV1 is acceptable for DAV2, giving the correct SE for most of the population parameters. The sparse design is efficient while allowing us to perform fewer samples per subject.

The results for power and NSN computation are reported in Table V. We observe a lack of power for the design with two periods and an increase of power for four periods to show equivalence between two treatment groups with 16 subjects per group. Indeed, with β Cl = 0.06 as estimated in the DAV1 study and a rich design, we predict 128 subjects in a two-period trial (involving a total number of 1792 samples to be performed) and 67 subjects in a four-period trial (1876 samples) needed to achieve a power of 90%. With the sparse design, we predict 131 subjects in a two-period trial (1048 samples) and 69 subjects in a four-period trial (1104 samples) needed to achieve a power of 90%.

Discussion

In this paper, we proposed the expression of M F for crossover trials analysed with NLMEM including WSV, in addition to BSV, and with discrete covariates which can change between periods. It is an extension of the developments proposed in [START_REF] Mentré | Optimal design in random effect regression models[END_REF], [START_REF] Retout | Further developments of the Fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics[END_REF] and [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF], using a linearisation of the model around the expectation of random effects. Simulations with various scenarios showed the correctness of the predictions of the SE and power for the comparison or equivalence Wald tests computed by M F . In spite of a slight underestimation of RSE by M F for some parameters, which respects the Rao-Cramer inequality, the discrepancy is not important (< 10%); indeed, when experimenters design population studies, it is more the amplitude of RSE than the exact value that they are looking for. This extension of M F can be concluded as relevant while saving lots of computation time compared to extensive simulations. We also derived an optimal sparse design and showed that it had power very close to that of the rich design, even with fewer samples per subject. This can be explained by the fact that there are only three parameters in the PK model, so a lower number of parameters compared to the number of samples in the sparse design (four samples per period). It is also what we obtained in our application example, which showed the importance of the number of subjects for the power of test. We notice that in the application, the covariate effect β Cl was estimated from the data of the DAV1 study between two periods although in that design, the piglets received amoxicillin plus placebo at the first period and amoxicillin alone at the second period. So β Cl expressed something other than a treatment effect, perhaps a period or a placebo effect. We chose however to include this discrete covariate in the statistical model as a treatment effect and to estimate this parameter because our objective in this application is to design the subsequent study DAV2 where there will really be a treatment effect from the compound X co-administrated with amoxicillin at the second period and similarly we assume that there will be no effect of drug X (not more than placebo). Then as presented in section 4.2, we examined not only the estimated value of β Cl as in DAV1 but also different ones to predict the power and NSN for DAV2. In this application, more than 16 piglets were needed in a one-way, two-period crossover design to show an absence of interaction of compound X on the PK of amoxicillin, which can be explained by the important standard deviation of WSV on Cl ( √ γ = 0.45). Increasing the number of periods improved the power of the study because the variance of the treatment effect parameter was reduced.

The extension of M F provides a useful tool to study the influence of the number of periods, of subjects or the effect size of a discrete covariate. We can also easily examine how the power would vary with different ratios between the BSV and WSV in a crossover design, and to compare these results to those obtained from a parallel group design. This emphasises the advantage of crossover design compared to parallel group design in bioequivalence trials except when WSV is the total variability and BSV = 0. It would also be interesting to include in the model a period or sequence effect and to consider different ratios between BSV and WSV to see its influence on the results. These extensions were implemented in the new version 3.2 of PFIM, available freely since January 2010 at www.pfim.biostat.fr, with several new features in terms of model specification and of M F expression. An output example of PFIM 3.2 for design evaluation of a two-way, two-period crossover trial is shown in the Supplementary Material2 . This new mathematical development of M F is applicable to single or multiple response models, for any number of periods and sequences, with different designs at each period. However, the present implementation of M F in PFIM was performed only for the case of trials with the same sampling times in each period, which is the design of our motivating example.

In our present work, we omitted the parameter tlag when designing the DAV2 study. But the fact that tlag was not taken into account in design optimisation has no important influence on the optimal design, all the more justified as its value is very small (≈ 0.37h), and before the first sampling time (0.5h) of the study. However, it would be useful, when the value of tlag is important, to evaluate M F in the presence of this parameter. In addition, data below LOQ were omitted for design evaluation in this extension but it would be interesting to derive an expression of M F including these left-censored data. An ad hoc method would be to take into account the percentage of data below LOQ at each sampling time, evaluated by simulation, in the evaluation of M F . It would be also interesting to consider the contribution of this information to the likelihood as proposed in [START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: application to HIV dynamics model[END_REF]. Also, here we considered only diagonal Ω and Γ matrices but one may want to allow correlations between parameters. The expression of M F for the full Ω matrix has been developed by [START_REF] Mentré | Optimal design in random effect regression models[END_REF] but has not been yet implemented in PFIM.

In this work, the development of M F was based on an approximation by linearisation of the model using the Taylor expansion. In our linearised calculation, we chose assume independence between the variance of the observations and the fixed effects as in a linear mixed model, which leads to a block diagonal expression of M F . Retout and Mentré [START_REF] Retout | Further developments of the Fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics[END_REF] and other authors [START_REF] Hooker | Simultaneous population optimal design for pharmacokineticpharmacodynamic experiments[END_REF][START_REF] Ogungbenro | Optimal design for multiresponse pharmacokinetic-pharmacodynamic models-dealing with unbalanced designs[END_REF][START_REF] Ogungbenro | Incorporating correlation in interindividual variability for the optimal design of multiresponse pharmacokinetic experiments[END_REF][START_REF] Gagnon | Optimal population designs for PK models with serial sampling[END_REF] have also proposed an expression of M F taking into consideration the dependence of the observations on the parameters of the model, leading to a complete M F with an additional off-diagonal bloc C. There is no clear consensus on what is the best approximation but we think that the block diagonal expression is better. Indeed, presently, Mielke and Schwabe [START_REF] Mielke | Some considerations on the Fisher information in nonlinear mixed effects models[END_REF] in addition showed that an approach with block diagonal M F was more reliable than the one with the full M F . In the software MONOLIX, the observed M F calculated by linearisation of the model around the individual parameter estimate has also a block diagonal expression. Also, the method presented here was shown to be relevant by simulation for the present example and in others [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: Evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF][START_REF] Guedj | Design evaluation and optimisation for models of hepatitis C viral dynamics[END_REF]. This can be perhaps explained by the fact that M F calculation is based on derivations of the model. The linearisation could perhaps introduce potential problems when the model is very nonlinear or very complex and written in differential equations. Consequently, one of the perspectives of this work would be to propose a computation of M F without model linearisation using a stochastic approach.

Studies analysed through NLMEM can be perfomed with optimal sparse simpling times with almost no loss of power. This requires the knowledge of the model and its parameters, a limitation also present in the approach by simulation. Sensitivity analyses with respect to the model and the parameter values would be necessary and interesting to quantify how the results possibly vary. However, at the designing step in a PK equivalence or interaction study of a drug, one can usually use the results of a previous study on the same drug so that the PK of the drug is already well known, and the model has already been evaluated. Also, this method allows of reducing the number of samples per subject, which can be ethically and practically very important for performing studies in patients. In conclusion, we have shown the relevance of the extension of M F in NLMEM for crossover trials. The approach by NLMEM is an appropriate alternative to NCA in the case of trials in patients. The implementation of these new developments in PFIM provides a useful tool to design bioequivalence/interaction studies as well as other kinds of longitudinal studies, avoiding extensive simulations.

Figure 2 -Figure 3 -

 23 Figure 2 -Boxplots of the observed SE for treatment effect parameter in each simulated scenario. (+) represents the value of the SE predicted by the extension of M F and (×) represents the empirical SE for the rich (white bar) or the sparse (grey bar) design.

Table I -

 I Predicted and empirical standard error SE(β Cl ) for each design and value of β Cl . Cl ) Predicted SE(β Cl ) [×10 2 ] Empirical SE(β Cl ) [×10 2 ]

	Design exp(β rich 0.8	3.401	3.290
	(0.5, 1, 1.5, 2, 4, 6, 8 h) 1	3.404	3.297
		1.1	3.405	3.300
		1.2	3.406	3.305
		1.25	3.407	3.308
		1.5	3.410	3.320
	sparse	0.8	3.443	3.357
	(0.5, 2, 6, 8 h)	1	3.454	3.372
		1.1	3.459	3.377
		1.2	3.462	3.390
		1.25	3.463	3.397
		1.5	3.467	3.423

Table II -

 II Predicted powers (%) and observed powers (%) for each design and value of β Cl .

	Design	exp(β Cl )	Comparison test	Equivalence test
			Predicted power Observed power Predicted power Observed power
	rich	0.8	100	100	5.00	5.60
	(0.5, 1, 1.5, 2, 4, 6, 8 h) 1	5.00	5.30	100.00	100.00
		1.1	79.92	82.80	98.25	98.90
		1.2	99.97	99.90	32.76	32.00
		1.25	100	100	5.00	4.30
		1.5	100	100	NC	NC
	sparse	0.8	100	100	5.00	5.60
	(0.5, 2, 6, 8 h)	1	5.00	5.10	100.00	100.00
		1.1	78.69	79.80	97.99	98.20
		1.2	99.95	99.80	32.08	34.10
		1.25	100	100	5.00	4.70
		1.5	100	100	NC	NC

(NC: not computed because β Cl is not under H 1 .)

Table V -

 V Prediction of power for the Wald test of equivalence and of NSN to achieve a power of 90% computed by PFIM 3.2 for each design of DAV2 study.

		0.07					Rich design Sparse design	Predicted SE Empirical SE
	Design		β Cl	Two periods AP/AX		Four periods AP/AX/AP/AX
		0.06		Predicted power (%) NSN for a Predicted power (%) NSN for a
				with N = 16 piglets power of 90% with N = 16 piglets power of 90%
	rich	0.05	0	41.0		68		64.3	34
	(0.5, 1, 1.5, 2, 4, 6, 8 h) 0.06	27.0		128		41.3	67
	sparse (0.5, 2, 4, 6 h) SE(β Cl )	0.03 0.04	0.10 0 0.06 0.10	19.3 40.5 26.7 19.1		209 70 131 213		27.7 63.7 40.8 27.4	124 69 126 35
		0.02
		0.01
		0.00
			0.8	1	1.1	1.2	1.25	1.5
					exp(β Cl )

Supporting information may be found in the online version of this article.

Supporting information may be found in the online version of this article.

Acknowledgments

The authors would like to thank the biotechnology company Da Volterra, 172 rue de Charonne, 75011 Paris, France (http://www.davolterra.com) for providing the data of the pharmacokinetic study DAV1 of amoxicillin in piglets. Conflict of Interest: None declared.


Van Minh Nguyen
email: vanminh.nguyen@alcatel-lucent.com

Holger Claussen
email: holger.claussen@alcatel-lucent.com

Efficient Self-optimization of Neighbour Cell Lists in Macrocellular Networks

Keywords: neighbour cell list, self-optimization, selfconfiguration, cellular networks

The neighbour cell list (NCL) in cellular networks has an important impact on the number of dropped calls and is traditionally optimized manually with the help of planning tools. In this paper, a method for automatically optimizing a NCL is presented, which consists of an initialization using a selfconfiguration phase, followed by a self-optimization phase that further refines the NCL based on measurements provided by mobile stations during the network operation. The performance of the proposed methods is evaluated for different user speeds and different NCL sizes. Besides, the convergence speed of the proposed self-optimization method is evaluated. It is shown that when about 6000 measurements are reported by mobile stations, the proposed self-optimization method attains a stable maximum performance about 99% of success rate.

I. INTRODUCTION

Handover is one of the most critical issues in cellular networks. It enables connection continuity for mobiles during move, while allowing the efficient use of resources (e.g., time and frequency reuse between cells). Most of today's cellular standards use mobile-assisted handover in which the mobile measures the pilot channel signal quality of neighbour cells and reports the measurement result to the network. If the signal quality from a neighbour cell is better than that of the serving cell by a handover margin, the network initiates a handover to that cell.

Today, commercial cellular networks use a neighbour cell list (NCL) to control the neighbour cell measurements by mobiles. In order to measure neighbour cells, the mobile needs to be provided with information on the pilot channels of potential neighbour cells (e.g., cell synchronization information, pilot channel frequency). If this information is not provided by the network, the mobiles have to spend a longer time to acquire it via scanning through all possibilities. Therefore, a list of neighbour cells as handover candidates is configured for each cell. This list contains information of the selected handover candidates and is sent to all mobiles connected to the cell. The mobiles then only need to monitor the pilot signal quality of the cells included the NCL of the serving cell.

The neighbour cell list has an important influence on the mobility performance. It must contain a sufficiently large number of potential neighbour cells to ensure that any mobile in the serving cell can find at least one handover target when its own signal deteriorates. However, as the mobile measurement capacity is limited, a long list can result in delays in finding a suitable handover target. These delays may cause call drops when the user moves at high speed.

While it might appear easy to configure a neighbour cell list by simply looking at the cell topology and selecting all cells overlapping with the cell of interest, it is not as straightforward in reality. This is mainly due to the fact that the real radio coverage is unknown, since it depends largely on the environment which can only be modelled to some degree in radio simulations. For instance, the presence of unpredicted buildings, trees, and moving objects will have effects on the coverage and may create coverage holes not predicted by the model. In addition, changes in the environment due to for example newly constructed buildings cannot easily be taken into account. Thus, an accurate prediction of the cell coverage is hard and also might not be valid anymore after some time.

In recent practice, the NCL is still manually configured at the beginning of the network deployment by means of planning tools, and is manually updated when new cells are installed. During this manual configuration the cell coverage and the neighbour relation are predicted using static information such base station locations, antenna patterns, and received signal maps (obtained through driver-tests). Due to the sensibility and the dynamic of the radio propagation condition, these static predictions of the cell coverage are more or less inaccurate, and cannot take changes in the radio conditions into account. More details of the drawbacks of the manual NCL configuration can be found in [START_REF] Olofsson | A concept for dynamic neighbor cell list planning in a cellular system[END_REF], [START_REF] Amirijoo | Neighbor cell relation list and physical cell identity self-organization in LTE[END_REF].

Therefore, the manual configuration and optimization of the neighbour cell list is a non-trivial problem that requires efforts for the network operator, resulting in considerable operational expenses. In [START_REF] Olofsson | A concept for dynamic neighbor cell list planning in a cellular system[END_REF], Olofsson et al. explained how a selfoptimization approach is able to overcome these difficulties. In addition, the Next Generation Mobile Network (NGMN) alliance recently specified neighbour cell list self-optimization as a required use case [START_REF]Use cases related to self-organising network: Overall description[END_REF], and described its benefits for the network operator [START_REF]Recommendation on SON and O&M requirements[END_REF].

A key element of the neighbour cell list configuration is reliable information on the cell coverage. In a real network, the best way to obtain accurate coverage information is to continuously perform measurements of the radio signals. The most efficient way to achieve this is to let the network entities such as mobile stations and base stations perform measurements during operation. In this paper we propose a method which uses live measurements to self-configure the initial neighbour cell list, and to continuously self-optimize the neighbour cell list of a cell during operation. The paper provides the following contributions:

• A complete optimization procedure which self-configures an initial NCL and continuously self-optimizes the NCL. • An efficient method which does not generate additional overheads by using the mobile-reported measurement available by the standard requirement. • A definition of the optimal neighbour cell list.

• Efficient algorithms which avoid static constraints and compute the NCL from live measurements. This is important as it frees the operator from a tricky manual collection of the network configuration, and allows the NCL to be quickly adapted to the changes in the network and its environment. The rest of this paper is organized as follows. In Section II, related work in this area is discussed. Section III describes the concept of the proposed method. The self-configuration method used for initialization is described in Section IV. In Section V, the ongoing self-optimization process is described. Here, during the network operation, connected mobiles measure neighbour cell pilot signals and report results to the serving base station. Then the NCL is optimized such that it maximizes the probability of finding an appropriate handover target. Section VI presents the performance evaluation using simulations based on a real network deployment scenario. The paper is concluded in Section VII.

II. RELATED WORK

The general concept of dynamic neighbour cell list planning in [START_REF] Olofsson | A concept for dynamic neighbor cell list planning in a cellular system[END_REF], [START_REF] Magnusson | Dynamic neighbor cell list planning in a microcellular network[END_REF] is one of the first attempts in the area of neighbour cell list self-optimization. In this concept, a neighbour cell is dynamically added to the NCL if the percentage of handovers from the serving cell to this neighbour cell is above a predefined threshold. Now and then some test frequencies are added to the neighbour frequency list to allow mobiles measuring cells which are not currently in the NCL. This aims at discovering potential good neighbours.

Following the motivation introduced in [1], [START_REF] Magnusson | Dynamic neighbor cell list planning in a microcellular network[END_REF], Guerzoni et al. [START_REF] Guerzoni | Automatic Neighbor Cell List Optimization for UTRA FDD Networks: Theoretical Approach and Experimental Validation[END_REF] proposed an automatic optimization algorithm. This algorithm mainly consists of identifying and ranking missing neighbours. The neighbours are identified analytically, and then ranked based on the base station coordinates and the antenna direction. In parallel, the cells currently in the NCL are assessed using key indicators such as handover success rate and received signal quality. After that, good cells from the pool of missing cells are added to the NCL while bad cells in the NCL are removed.

Realizing that the method proposed in [START_REF] Guerzoni | Automatic Neighbor Cell List Optimization for UTRA FDD Networks: Theoretical Approach and Experimental Validation[END_REF] may produce a very large pool of potential missing neighbours, and that the statistical confidence on the parameters used for ranking the cells may not be optimal, Soldani and Ore [START_REF] Soldani | Self-optimizing neighbor cell list for UTRA FDD networks using detected set reporting[END_REF] proposed an enhanced procedure for UTRA FDD networks. Here, detected set reporting is used to identify missing neighbours, then thresholds are used as a coach filter to limit the pool list.

The objective here is to efficiently identify and filter missing cells, but not to compute an optimal NCL.

Further work by Parodi et al. [START_REF] Parodi | An automatic procedure for neighbor cell list definition in cellular networks[END_REF] focused on automatically creating a NCL for newly installed cells. Slightly similar to [START_REF] Soldani | Self-optimizing neighbor cell list for UTRA FDD networks using detected set reporting[END_REF], cells located in a given range from the new cell are filtered out as potential neighbours. Then, based on prediction models, each cell in the pool is checked for overlap with the new cell. A number of the best neighbours are finally selected for the NCL. In addition, a signalling implementation for LTE networks is proposed to allow a newly installed cell to detect its neighbours.

Three self-configuration methods for WiMAX networks are proposed in [START_REF] Li | On the study of self-configuration neighbour cell list for mobile WiMAX[END_REF] and are very similar to those discussed in [START_REF] Parodi | An automatic procedure for neighbor cell list definition in cellular networks[END_REF].

In [START_REF] Amirijoo | Neighbor cell relation list and physical cell identity self-organization in LTE[END_REF], [START_REF] Amirijoo | Neighbor cell relation list and measured cell identity management in LTE[END_REF], live measurements are considered for the neighbour relation table optimization and physical cell identifier (PCI) conflict detection in LTE.

The above studies contributed valuable concepts for the self-optimization of the NCL. However, when realizing the proposed concepts, they are mostly based on static information and static thresholds which limit the performance and the ability to adapt to dynamic changes in the radio environment.

III. CONCEPT

The concept of the proposed NCL self-optimization method is conceived on the basis of the following questions:

• What is the definition of an optimal NCL? • What measurement statistics are needed and how can they be obtained? • How can potential neighbours be discovered that are not currently included in the neighbour cell list? The objective of the neighbour cell list is to assist mobile stations in finding a suitable handover target. A neighbour cell is considered as a suitable handover target for a mobile if its pilot signal quality received by the mobile is higher than a minimum requirement P min , for example P min,Ec/I0 = -20 dB in WCDMA [START_REF]Requirements for support of Radio Resource Management FDD (Release 8)[END_REF]. In addition, the information obtained by measurements of the pilot signals of neighbour cells must be of good statistical confidence. The shorter the neighbour cell list is, the more frequent the neighbour cells can be measured, which results in a higher confidence in the measurements. In reality, the maximum size of the neighbour cell list is limited by the standard, such as to 32 cells as specified by 3GPP. From the above considerations, the optimal neighbour cell list is defined as the one that has the least cells within the standardized allowable number (e.g., 32), and maximizes the probability of finding a handover target. Here the probability of finding a handover target is defined as the probability that at any location in the serving cell the mobile finds at least one suitable handover target from the neighbour cell list.

The measurement statistics needed are the indicators that are used for the handover target selection such as SINR (signal-to-interference-plus-noise ratio) in GSM, and E c /I 0 in WCDMA. Since this information is only available when the network is operational but not in the planning phase, the proposed optimization method is composed of two phases: When the network is newly deployed, neighbours are identified and measured by each base station on the downlink band. The neighbour cell list of each cell is then initialized using a self-configuration algorithm based on the collected measurements. This phase is also applied to self-configure a initial neighbour cell list of a new cell that is installed when the network is already operating.

During the network operation, connected mobiles are instructed to perform measurements using standardized scanning procedures, and report results to the serving base station. This phase of collecting measurement statistics also deals with the cell discovery. In broadband networks such as WCDMA, mobiles are capable of detecting and reporting cells even when they are currently not included in the neighbour cell list1 [START_REF]Requirements for support of Radio Resource Management FDD (Release 8)[END_REF], [START_REF]Requirements for Evolved UTRA and Evolved UTRAN[END_REF]. For multi-carriers networks like GSM, the cell discovery can be done by using the concept proposed in [START_REF] Magnusson | Dynamic neighbor cell list planning in a microcellular network[END_REF].

The collected measurement statistics are then used by a self-optimization algorithm implemented in the base station to continuously refine the neighbour cell list. The update of the current neighbour cell list with the self-configured or selfoptimized neighbour cell list can be totally automatic, with the intervention of the operator to allow the introduction of additional specific policies (for example to remove prohibitive cells is not performed fully autonomously).

IV. NEIGHBOUR CELL LIST SELF-CONFIGURATION

Just after the network planning phase, no measurement from mobiles are available, neighbour cells are identified based on measurements performed by each base station. This is done by using the existing sectorized antenna of the serving base station to scan the downlink band and measure identified cells. Then a predefined number of neighbours with the best received pilot quality are included in the neighbour cell list.

The serving base station's antenna gain supports the measurements of neighbour cells in the direction of the sector coverage, and thereby helps to detect neighbour cells which overlap with the sectorized serving cell.

This approach does not require the operator to provision information on neighbour cells, and thereby reduces the deployment costs. However, its disadvantage is that the neighbours pilot signal quality can not be measured in the whole serving cell coverage area. This is why further measurements by the mobiles are required. These measurements become available during operation and are used in the ongoing self-optimization process described next.

V. NEIGHBOUR CELL LIST SELF-OPTIMIZATION

This section details the measurement statistics retrieval and the self-optimization algorithm.

A. Neighbour Measurement Statistics Retrieval

During the scanning, the connected mobiles measure the pilot signal of neighbours currently in the neighbour cell list. In parallel, potential neighbours are detected and measured by the mobiles according to the principle explained in Section III. This means that the mobiles perform measurements of neighbours that are either detected during the scanning procedure, or are identified during scanning of test frequencies. The mobiles then send measurement reports to the serving base station.

The serving base station stores a measurement statistics table. The table columns are neighbour cells identifiers, and the table rows are measurement reports sent by mobiles. The value at row i and column j represents the pilot signal quality of cell j received by a mobile reporting the measurement i.

Whenever a measurement report is received, a new row is added. If this measurement report contains cells that are not currently present in the table columns, one new column is added for each new reported cell. After that, for the recently added row, columns corresponding to the reported cells are filled with the corresponding reported values.

Each row of the statistics table reflects the coverage of neighbour cells at the reported measured location. Any mobile may send several measurement reports during a scanning procedure, and these may correspond to different measured locations in the serving cell (due to user mobility). Therefore, only the measurement reports are relevant while the identifier of the mobiles performing the measurements is not stored.

Besides, old measurements may be not valid anymore due to potential changes in the radio environment. Therefore, at some point, older rows can be removed from the statistics table.

B. Self-Optimization Algorithm

The proposed self-optimization algorithm is presented in Algorithm 1. The objective of the algorithm is to select neighbour cells such that the resulting neighbour cell list maximizes the number of locations at which there is at least one suitable target cell from the neighbour cell list available.

The input includes the statistics table, the minimum requirement of the pilot signal quality P min , and the maximum Algorithm 1: Self-optimization algorithm Input: statistics table, P min , N max ; 1 initialisation: S = statistics table; NCL = ∅; 2 while N max > size of NCL do allowable size N max . The "for" loop in lines 3-8 cycles through the number of locations (i.e., number of table rows), m j , at which neighbour cell j is a suitable handover target.

Between the remaining neighbour cells in table S, the cell corresponding to column j * (line 10) is the best one that maximizes the number of locations covered by the neighbour cell list, and thus is added to the neighbour cell list (line 11). The case of max j (m j ) = 0 ("else" in line 16) means that all the remaining neighbour cells in table S do not cover any location of the serving cell with acceptable signal quality. Adding any of these cells will not improve the neighbour cell list. In this case, the resulted neighbour cell list is the best one that maximizes the probability of finding handover target, and so the algorithm is terminated (lines 17 and 20).

When the neighbour cell identified by column j * is included in the NCL (line 11), this column j * and the corresponding rows can be removed (by line 15, and 13, resp.). The cell identified by column j * ensures that the mobiles will find at least one suitable handover target at the locations (i.e., rows) that satisfy the "if" condition in line 13. As the objective is to find the shortest neighbour cell list that covers a maximum number of locations, only the locations that are still not covered by the current NCL are considered in the next rounds.

VI. SIMULATION

This section provides the system model and assumptions, and the simulation results.

A. Simulation Scenario

A cellular network deployment in Phoenix in the United States was used for the simulation, from which cell locations and antenna orientations were extracted. The size of the simulated scenario is 16 km × 16 km and includes 120 macrocell sectors. The propagation path loss for a typically urban environment is modelled as

L p[dB] (d) = 28 + 35 log 10 d, ( 1
)
where d is the distance from the base station in metres. Shadow fading L s[dB] in dB scale is modelled as random process with normal distribution with 8 dB standard deviation and a correlation of 0.5 between different sites. It is assumed that the total transmit power for each base station is P tx = 20 W per sector, where 10% of this power is allocated to a pilot channel P tx,pilot . Each macrocell transmits via a sectorized antenna with a gain pattern of

G [dB] (θ) = G max -min 12 θ β 2 , G s , -π < θ ≤ π, (2)
where β = 70π/180 is the angle where the gain pattern is 3 dB down from peak, the maximum gain is G max = 16 dB, and G s = 20 dB.

For any location in the considered scenario the received power from a base station j can be calculated based on (1) as

P rx[dBm] = P tx[dBm] + G [dB] (θ) -L p[dB] (d) -L s[dB] . (3)
The receiver noise power N 0 is modelled as

N 0[dB] = 10 log 10 (kT × N f × W ), (4)
where the effective noise bandwidth W = 3.84 × 10 6 Hz, kT = 1.3804 × 10 -23 × 290 W/Hz, and N f [dB] = 7 dB referring to the noise figure at the UE.

For each location in the scenario, an E c /I 0 can be calculated. The pilot E c /I 0 is defined in 3GPP as the ratio of the received energy per pseudo-noise chip for the CPICH to the total received power spectral density at the UE antenna connector. It can be calculated for a UE connected to the mth cell as

CPICH E c /I 0 i = P rx,pilot,i N 0 + j P rx,j , (5)
where P rx,pilot,i denotes the pilot power received from cell i, P rx,j is the total received power of cell j including data and pilot. The pilot minimum requirement for the UE is CPICH E c /I 0 ≥ -20 dB [START_REF]Requirements for support of Radio Resource Management FDD (Release 8)[END_REF]. In the scenario considered, this requirement results in a nearly 100% coverage of the area. Fig. 2(a) shows an example of the highest received pilot E c /I 0 , and Fig. 2(b) shows the corresponding coverage regions based on received pilot E c /I 0 for the sectors within the simulated area. It is assumed that the UE is capable of measuring CPICH of 8 identified cells (i.e., cells in the NCL) within 200 ms [START_REF]Requirements for support of Radio Resource Management FDD (Release 8)[END_REF].

For the NCL performance evaluation, mobiles are randomly placed in the simulated area, and move in the network. Each mobile moves at a given speed and changes its direction after each simulated time step (see, e.g., [START_REF] Zonoozi | User mobility modeling and characterization of mobility patterns[END_REF], [START_REF] Wang | Joint effects of radio channels and node mobility on link dynamics in wireless networks[END_REF]). User speeds between 5 m/s and 50 m/s are considered. If the signal quality received from the serving cell falls below the scanning trigger threshold, the mobile starts scanning cells comprised in the

Note that after finding a handover target, the mobile enters into the handover execution phase whose outcome depends on the handover algorithm but is not impacted by the NCL. Therefore, only the handover target searching phase is used to assess the NCL quality, the handover execution phase is not taken into the NCL performance evaluation.

B. Performance

One example of NCL optimization and handover target searching is shown in Fig. 4 for the serving cell 3. Fig. 4(a) shows the results directly after initialization with a NCL computed by the self-configuration algorithm, and Fig. 4(b) shows the results with the NCL that is computed by the selfoptimization algorithm. The optimized NCL has 21 cells. For comparison, a cell list size of 21 was also used for the selfconfiguration algorithm. In Fig. 4, cells comprised in the NCLs 4(a) and 4(b), the selfoptimization method can improve the success rate of 85% from the self-configuration to 97.3%. It is shown that three cells 7, 96 and 98 make the difference between the two NCLs. By using mobile-reported measurements, the self-optimization algorithm included these three cells, and thereby eliminates call drops close to these neighbour cells.

For more general results, the handover target searching with NCLs to be evaluated is simulated in 20 serving cells, and the average success rate of each simulated NCL is computed. The self-configuration method is compared to two other methods. The first one is a distance-based method that configures a NCL by selecting a number of the closest neighbour cells. And the second one is similar to the proposed self-configuration method except that the measurements are performed using a separate antenna with uniform gain pattern. These benchmarking methods are marked in Fig. 5(a) & 5(b) as "distancebased" and "uniform antenna gain", respectively. Note that, as discussed in Section II, self-optimization concepts in the literature are not realized with live measurements like the proposed approach, but with static configuration information. Therefore, no directly comparable benchmark is available. Instead, the self-optimization method is compared with the best self-configuration method for its performance.

Fig. 5(a) shows the performance with different user speeds v. It is shown that with a scanning speed of 8 cells/200ms as in WCDMA, the user speed does not have much impact on the success rate. It is also shown that the proposed self-configuration method is better than two benchmarking methods, and the self-optimization method performs best. Fig. 5(b) shows the performance with different list sizes. We see that all the NCLs give better performance for longer lists. This is thanks to the relatively high scanning speed in WCDMA. In some case where scanning speed is low, if the NCL is too long the connection between the mobile and the serving may be lost before finding a suitable handover target. It is also clear from Fig. 5(b) that the proposed self-configuration method is better than two benchmarking methods, and the selfoptimization again performs best. It is shown that the success rates flatten out at about 10 neighbours. From this point on, the performance of the self-optimized NCL approaches a success rate close to 100%, while the configured NCLs have lower

C. Convergence Speed of the Self-Optimization Method

A further important aspect that was evaluated is the convergence speed of the self-optimization method. This provides information on from which point the self-optimization method provides reliable performances that outperform the self-configuration method used for initialization. Fig. 5(c) shows the achievable performance for different numbers of reported measurements. It is shown that the more measurements are available, the better the quality of the self-optimized NCL is, since more reliable information about the neighbour cell coverage is available.

At about 3000 measurements, the self-optimization algorithm provides a success rate about 95%. And at 6000 available measurements, the success rate attains 99%. Since the selfconfiguration achieves an average performance of approximately 95% when the list size is sufficient large as shown in Fig. 5(b), it makes sense to switch to the NCL obtained from the self-optimization algorithm once more than approximately 3000 measurement are available. From this point on, the selfoptimized NCL provides improved performance.

VII. CONCLUSION

A method for automatically optimizing a neighbour cell list was presented. The proposed approach consists of a selfconfiguration phase for initialization, followed by a switch to a self-optimization phase once sufficient measurements are reported by connected mobile stations. For initialization, the measurements-based self-configuration method achieved the best performance of approximately 95% in handover target searching success rate with 10 or more neighbour cells. The success rate is further improved by switching to a NCL obtained by the proposed self-optimization method. It was shown that a switch to a self-optimized NCL can improve the performance once more than approximately 3000 measurements reports are collected, and a near maximum performance of around 99% success rate is achieved after more than around 6000 measurements reports are available. The proposed solutions are simple but efficient by using true measurementbased self-optimization algorithms. These suggest that the optimization of NCLs can be automated, which will help to reduce operational expenses in future cellular networks.

Fig. 1 .

 1 Fig. 1. Principle of neighbour cell list self-optimization

3 foreach column j of S do 4 m j = 0; 5 foreach 9 if

 3459 row i of S do 6 if S[i, j] ≥ P min then m j ← m j + 1; max j (m j ) > 0 then 10 j * = arg max(m j ); 11 add cell identified by column j * to the NCL; 12 foreach row i of S do 13 if S[i, j * ] ≥ P min then remove row i from S;

Fig. 3 .

 3 Fig. 2. Simulation scenario

Fig. 4 .

 4 Fig. 4. Examples of scanning from serving cell 3. User speed v = 20 m/s

Fig. 5 .

 5 Fig. 5. Performance evaluation. Simulated with 20 serving cells.

This ability does not weaken the importance of the NCL since, as mentioned in the introduction, the mobiles need longer time to measure an unidentified neighbour cell.

ACKNOWLEDGMENT

This work was partially supported by POSEIDON project within the French cluster System@tic.