Cyrielle Dumont
email: cyrielle.dumont@inserm.fr

Marylore Chenel

France Mentré

TWO-STAGE ADAPTIVE DESIGNS IN NONLINEAR MIXED EFFECTS MODELS: APPLICATION TO PHARMACOKINETICS IN CHILDREN

Keywords: Adaptive design, design optimization, nonlinear mixed effects models, Population pharmacokinetics

INTRODUCTION

In clinical pharmacokinetics (PK), the evolution of drug concentration versus time is studied. PK profiles are described by mathematical models, and nonlinear mixed effects models (NLMEM) [START_REF] Sheiner | Modelling of individual pharmacokinetics for computeraided drug dosage[END_REF][START_REF] Sheiner | Pharmacokinetic/pharmacodynamic modelling in drug development[END_REF] are used to analyze PK data obtained in group of patients or healthy volunteers. This is called the population approach. PK studies in children are often more difficult to perform than those in adults. The blood volume that can be taken in children is much more limited than in adult healthy volunteers so that the number of sampling times is limited. For ethical reasons, PK studies in children are conducted in patients who may potentially benefit from the treatment, and not in healthy volunteers, as for adults. As, most of the time, only sparse data can be obtained in children, NLMEM is an appropriate methodology for analyzing PK information in pediatric trials both from a practical and ethical point of view (EMA, 2006 ;Mentré et al., 2001 ;[START_REF] Tod | Facilitation of drug evaluation in children by population methods and modelling[END_REF].

The design of PK studies is important when the number of samples and subjects is limited, as every sample must be informative, especially when clinical constraints are strong. Indeed, the design has a large impact on the precision of population parameter estimates [START_REF] Mentré | Optimal design in random-effects regression models[END_REF]. A PK design for NLMEM consists in choosing the number of patients and for each patient the elementary design. Elementary designs are composed of several sampling times to be drawn for each individual. In the EMA guideline related to PK studies in the pediatric population (EMA, 2006), simulations or theoretical optimal design approaches, based on prior knowledge, are presented as tools to be considered for the selection of sampling times, number of subjects and number of samples per subject. To avoid simulations, which are time consuming, designs can be evaluated using the Fisher information matrix (M F ) and the optimization of its determinant, corresponding to the D-optimality criterion. The calculation of M F for NLMEM was first developed by [START_REF] Mentré | Optimal design in random-effects regression models[END_REF] and [START_REF] Retout | Fisher information matrix for non-linear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF] for uniresponse NLMEM and then extended to multiresponse population PK/pharmacodynamic (PD) models [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation and appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF][START_REF] Gueorguieva | A program for individual and population optimal design for univariate and multivariate response pharmacokinetic and pharmacodynamic models[END_REF][START_REF] Retout | Fisher information matrix for non-linear mixed-effects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF] using a first-order Taylor expansion of the population PK model around the random effect [START_REF] Mentré | Optimal design in random-effects regression models[END_REF][START_REF] Retout | Development and implementation of the population Fisher information matrix for the evaluation of population pharmacokinetic designs[END_REF]. Several approximations of M F can be used in the community of PKPD design, but it was suggested [START_REF] Mielke | Some considerations on the Fisher information in nonlinear mixed effects models[END_REF][START_REF] Mielke | Approximations of the Fisher information for the construction of efficient experimental designs in nonlinear mixed effects models[END_REF] that the block diagonal expression might be more reliable than one with the full M F , when an approach by linearization is performed. The calculation of the M F for NLMEM used in population PK is implemented in several software packages, including PFIM developed in R, dedicated to design evaluation and optimization [START_REF] Bazzoli | Design evaluation and optimisation in multiple response non linear mixed effects models: PFIM 3.0[END_REF]. PFIM evaluates and optimizes population designs in NLMEM with single and multiple responses [START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation and appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF]. The expression of M F in NLMEM was also extended to evaluate and to optimize designs in crossover trials [START_REF] Nguyen | Design evaluation and optimisation in crossover pharmacokinetic studies analysed by nonlinear mixed effects models[END_REF] and is implemented into the version PFIM 3.2, available since 2010. To evaluate M F in NLMEM, prior information is needed such as the model but also a priori values of the parameters. This is called a local design [START_REF] Burstein | Evaluation of a sparse sampling strategy for determining vancomycin pharmacokinetics in preterm neonates: application of optimal sampling theory[END_REF]. Local design based on these a priori values can lead to a sub-optimal design when the true population parameter values are different from the a priori ones. Alternatives to local design are robust design, relying on a priori distribution of parameters [START_REF] Dodds | Robust population pharmacokinetic experiment design[END_REF]Foo et al., 2012 ;[START_REF] Pronzato | Robust experiment design via maximin optimization[END_REF], or adaptive design [START_REF] Chow | Adaptive design methods in clinical trials[END_REF].

Adaptive design is increasingly used for randomized clinical trials or for dose-ranging studies [START_REF] Dragalin | A simulation study to compare new adaptive dose-ranging designs[END_REF]EMA, 2009 ;FDA, 2006 ;FDA, 2010). The working group Pharmaceutical Research and Manufacturers of America (PhRMA, [START_REF] Dragalin | Adaptive designs: terminology and classification[END_REF]) defined an adaptive design as a design of experiment in which accumulated data during the trial are used to possibly modify the aspects of the study, without compromising the validity and integrity of the clinical trial [START_REF] Gallo | Adaptive designs in clinical drug development -an executive summary of the PhRMA Working Group[END_REF].

It was already shown that adaptive designs improve dose ranging during clinical drug development [START_REF] Bornkamp | Innovative approaches for designing and analyzing adaptive dose-ranging trials[END_REF]. Adaptive designs are now widely developed for randomized clinical trials or dose-ranging studies, but are barely used in population PK/PD (Foo and Duffull, 2012 ;[START_REF] Zamuner | Adaptive-Optimal Design in PET occupancy Studies[END_REF]. However, according to a survey performed by [START_REF] Mentré | Current use and developments needed for optimal design in pharmacometrics: a study performed amongst DDMoRe's EFPIA members[END_REF], adaptive design in population PK/PD is a priority for pharmaceutical companies, graduating its importance with a median equal to 4 on a scale between 0 and 5. [START_REF] Zamuner | Adaptive-Optimal Design in PET occupancy Studies[END_REF] showed that optimal adaptive design provided a more efficient design for the study they conducted in imaging. Indeed, the optimal adaptive design allowed to minimize the number of subjects needed for the study and to maximize information. Foo and Duffull (2012) developed and evaluated a method allowing to conduct PK adaptive bridging studies. Using simulation studies, the authors showed that D-optimal adaptive method applied to bridging studies was more efficient, since it provided better estimations, when the PK profiles of the target and a priori populations differed.

An important point in adaptive design is the number of adaptations, that is cohorts of patients. In a study concerning the dose-ranging optimal design, it was shown that twostage adaptive designs are more efficient than fully adaptive designs, i.e., when adaptation is performed after each patient [START_REF] Fedorov | Optimal dose-finding designs with correlated continuous and discrete responses[END_REF]. In the framework of a clinical trial in cancerology, Chen (1997) compared three-stage adaptive designs to two-stage adaptive designs. Comparison showed that the benefit of a three-stage adaptive design compared to a two-stage design is not as important as the gain of two-stage design versus one stage. They also discussed that an adaptive design with four stages or more would not allow to improve the results. Furthermore, two-stage adaptive designs are easier to conduct in clinical trials as only one modification is made.

The first objective of the present work was to introduce two-stage adaptive design in NLMEM. In section 2, we present the model and the principle of the two-stage adaptive design, and the development of M F associated to this design is described in Section 3. The second objective was to evaluate, with a simulation approach, the impact of two-stage designs on the precision of parameter estimation when the true PK parameters are different from the a priori ones. We have performed this simulation for a population PK study in paediatrics. And finally, we have investigated the influence of the sample size ratio between the two stages of an adaptive design. Section 4 presents the simulation study and the results are given in Section 5.

Model and design

A NLMEM, or a population model, is defined as follows. The elementary design ξ i of individual i (i = 1, ..., N ), is defined by n i sampling times and their allocation in time, that is ξ i = (t i1 , ..., t in i ). The vector of observations Y i for the i th individual is defined as

Y i = f (g(β, b i ), ξ i ) + ε i (1) with g(β, b i ) = β + b i (2)
for a normal distribution of parameters (additive random effects), or

g(β, b i ) = βexp(b i ) (3)
for a lognormal distribution of parameters (exponential random effects). β represents the pvector of the fixed effects parameters and b i , the vector of the p random effects for individual i. It is assumed that b i ∼ N (0, Ω) with Ω defined as a p × p diagonal matrix, for which, each diagonal element ω 2 r , r = 1, ..., p, represents the variance of the r th component of the vector b i . The function f defines the nonlinear structural model. ε i is the vector of residual error and it is also supposed that ε i ∼ N (0, Σ i ) with Σ i a n i × n i -diagonal matrix such that:

Σ i (β, b i , σ inter , σ slope , ξ i ) = diag(σ inter + σ slope × f (g(β, b i ), ξ i )) 2 .
(4)

The parameters σ inter and σ slope are respectively the additive and proportional parts of the error model. The case σ slope = 0 corresponds to a homoscedastic error model, whereas the case σ inter = 0 corresponds to a constant coefficient of variation error model. The case where the two parameters differ from zero is called a combined error model. Finally, conditionnally on the value of b i , it is assumed that the errors ε i are independently distributed. Let λ T = (ω 2 1 , ..., ω 2 p , σ inter , σ slope ) be the vector of variance terms and let Ψ be the vector of the population parameters to be estimated, so that Ψ T = (β T , λ T ).

A one-stage population design for NLMEM is composed of N individuals to whom an elementary design ξ i , i = 1, ..., N , is allocated. A population design is therefore described by N elementary designs for a total number n tot of observations such as n tot = N i=1 n i and Ξ = {ξ 1 , ..., ξ N }. In the case where the elementary design ξ is the same for the N individuals, the total number of observations is n tot = N × n.

A two-stage population design is composed of 2 groups of N 1 and N 2 individuals (N 1 + N 2 = N ) with respectively designs Ξ 1 and Ξ 2 within each group. When the design is assumed identical within each group, the population design can then be written as follows:

Ξ = {[ξ 1 , N 1 ]; [ξ 2 , N 2 ]}.
In a two-stage adaptive design, Ξ 1 is determined from a priori parameters Ψ 0 and Ξ 2 is optimized for the parameters Ψ1 estimated with the data of the first cohort. More precisely, at the first stage, from the model and a priori parameters Ψ 0 , data are collected for N 1 individuals according to a design Ξ 1 , optimized using the a priori parameters. From data Y 1 , obtained with the first cohort, the population parameters Ψ1 are estimated. At the second stage, the design Ξ 2 of the second cohort, composed of N 2 individuals, is then optimized for the parameters Ψ1 , estimated at the end of the first stage. Data Y 2 are then collected for the N 2 individuals having the design Ξ 2 and finally, data Y 1 and Y 2 , obtained from N = N 1 + N 2 individuals included in each cohort, are simultaneously analyzed to estimate parameters Ψ2 .

Here, we use the D-optimality criterion to optimize designs which consists in maximizing the determinant of the Fisher information matrix (det(M F )). We suppose that all individuals have the same elementary design ξ for a same stage; that is ξ 1 corresponds to Ξ 1 and ξ 2 to Ξ 2 . The sizes of cohorts N 1 and N 2 are fixed, as well as the number of samples n 1 and n 2 within each elementary design. Therefore, we only optimize with respect to the sampling times in ξ 1 and ξ 2 . Figure 1 represents a two-stage adaptive design.

Fisher information matrix

In the following, the index i for the individual is omitted for sake of simplicity. The elementary Fisher information matrix M F (Ψ, ξ) for an individual with design ξ is given by

M F (Ψ, ξ) = E(- ∂ 2 L(Ψ; Y ) ∂Ψ∂Ψ T ) (5)
where L(Ψ; Y ) is the log-likelihood of the vector of observations Y of that individual for the population parameters Ψ. Because f is nonlinear, there is no analytical expression for the loglikelihood and a first-order Taylor expansion [START_REF] Mentré | Optimal design in random-effects regression models[END_REF] of the model f (g(β, b), ξ), around the expectation of b, that is to say around 0, is used. Using this linearization, the statistical model can be written as:

Y ∼ = f (g(β, 0), ξ) + ( ∂f T (g(β, b), ξ) ∂b ) b=0 b + ε. (6)
The approximated marginal expectation E and variance V of Y are given by:

E(Y ) ∼ = E = f (g(β, 0), ξ) (7) V ar(Y ) ∼ = V = ( ∂f T (g(β, b), ξ) ∂b ) b=0 Ω( ∂f (g(β, b), ξ) ∂b T ) b=0 + Σ(β, 0, σ inter , σ slope , ξ) (8)
The log-likelihood L is approximated by:

-2L(Ψ; Y ) ∼ = nln(2π) + ln(|V |) + (Y -E) T V -1 (Y -E). (9)
For a linear mixed effects model, this tacitly assumes that the random variable Y is approximately Gaussian. Hence for NLMEM the approximation is better when the level of non-linearity of the model is limited and/or when inter-individual variability is small. The elementary M F depends on the approximated marginal expectation E and variance V of the observations. Assuming that the derivative of V does or does not depend on the fixed effects, the elementary M F is a full matrix or a block diagonal matrix. There is no clear consensus on what is the best approximation, but here in our approach by linearization, 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][START_REF] Nyberg | Methods and software tools for design evaluation in population pharmacokinetics-pharmacodynamics[END_REF]. With that approximation, we have:

M F (Ψ, ξ) ∼ = 1 2     A(E, V ) 0 0 B(E, V )     (10)
where

(A(E, V )) ml = 2 ∂f T ∂β m V -1 ∂f ∂β l (11)
with m and l = 1, ..., p and

(B(E, V )) ml = tr( ∂V ∂λ m V -1 ∂V ∂λ l V -1 ) (12)
with m and l = 1, ..., dim(λ).

In the case of only one group composed of N individuals, the population Fisher information matrix for a population design Ξ is computed as:

M F (Ψ, Ξ) = N × M F (Ψ, ξ). ( 13
)
This calculation of M F is implemented in PFIM [START_REF] Bazzoli | Design evaluation and optimisation in multiple response non linear mixed effects models: PFIM 3.0[END_REF].

When neglecting the dependence implied by the adaptation in the representation of the Fisher matrix, the two-stage design is performed as follows. For the two-stage adaptive design, we assume N 1 individuals associated with a same design ξ 1 at the first stage, and N 2 individuals associated with a same design ξ 2 at the second stage. For the first stage of a two-stage adaptive design, the Fisher information matrix is expressed as follows:

M F (Ψ 0 , N 1 ξ) = N 1 M F (Ψ 0 , ξ). ( 14
)
We define ξ 1 as the design which maximizes the determinant of this matrix. Then for the second stage, the expression of M F is:

M F ( Ψ1 , N 1 ξ 1 + N 2 ξ) = N 1 M F ( Ψ1 , ξ 1 ) + N 2 M F ( Ψ1 , ξ). ( 15
)
For the second stage, ξ 2 is the design which maximizes the determinant of this matrix, using the parameters Ψ1 estimated at the first stage. This expression is implemented in version 4.0 of PFIM which will be soon available.

SIMULATION STUDY

Pharmacokinetic example

The simulation study mimicked a pediatric PK trial. The PK model is a two-compartment model with first-order absorption, exponential random effects and proportional error model.

The observed concentration of a child at sampling time t j is modelled by:

f (θ, t j ) = D(Ae -α.t j + Be -β.t j -(A + B)e -ka.t j ) ( 16
)
where

θ = (ka, CL, V 1, Q, V 2) and β = 1 2 Q V 1 + Q V 2 + CL V 1 - Q V 1 + Q V 2 + CL V 1 2 -4 Q V 2 CL V 1 , α = Q V 2 CL V 1 β , A = ka V 1 Q V 2 -α (ka -α)(β -α) , B = ka V 1 Q V 2 -β (ka -β)(α -β)
.

In equation ( 16), D denotes the dose, ka the absorption rate constant, CL the clearance of elimination of the drug, Q the inter-compartmental clearance from compartment 1 to compartment 2, V 1 and V 2 the volumes respectively of the central compartment and of the peripheral compartment. With the chosen parametrization, note that the parameters vary independently of each other.

Here we assumed that the dose administered is equal to 0.1 mg. We defined two sets of pharmacokinetic parameters: the a priori parameters Ψ 0 , guessed from adults, and the true parameters Ψ * with which all simulations are performed. They are presented in Table 1. All parameters have the same variance (ω 2 = 0.3) and we consider the same proportional error (σ = 0.2).

The optimal designs with one or two stages are obtained assuming a total number N = 60 children, having the same design ξ at a same stage composed of 5 sampling times among 15 possible sampling times: 0.083, 0.17, 0.25, 0.33, 0.5, 0.75, 1, 2, 3, 4, 5, 6, 8, 10, 12 hours after drug administration. Designs are optimized with PFIM (version 3.2), which allows us to maximize the determinant of M F using statistical design within a finite set of sampling times [START_REF] Retout | Design in nonlinear mixed effects models: Optimization using the Fedorov-Wynn algorithm and power of the Wald test for binary covariates[END_REF]. Designs are then evaluated through the criterion given by det(M F ) 1/P , P being the total number of parameters in Ψ. We denote by ξ 1 and ξ * the designs optimized respectively with the parameters Ψ 0 and Ψ * , and by ξ 2 the design optimized at the second stage. Optimal designs obtained with the parameters Ψ 0 and Ψ * are respectively ξ 1 = 0.083, 1, 2, 5, 12 hours after drug administration and ξ * = 0.083, 0.33, 0.75, 2, 12 after drug administration. Figure 2 represents mean profiles obtained with the parameters Ψ 0 and Ψ * , on which are added the times of ξ 1 and ξ * .

Simulation settings

We consider designs with a total number of N = 60 children. All simulated data are obtained with the true parameters Ψ * . We are interested in 2 one-stage designs (non adaptive designs), that is the designs ξ 1 (60-0), optimized with parameters Ψ 0 , and the optimal design ξ * (0-60), optimized with the true parameters (Figure 2). In parallel, we are interested in 3 two-stage adaptive designs: ξ 50-10 (50-10), ξ 30-30 (30-30), ξ 10-50 (10-50), varying the size of each cohort N 1 and N 2 : 50 and 10, 30 and 30, 10 and 50 children, at the first stage and at the second stage, respectively.

For each of the designs mentioned above, we simulate 100 adaptive clinical trials with R. The simulation study is performed as follows. From the same design ξ 1 , we simulate 10 data sets (Y 1 1 , ..., Y 10 1 ) for the first cohort composed of N 1 children with the parameters Ψ * . From these 10 data sets, 10 vectors of parameters Ψ1 are estimated using saemix (http://cran.rproject.org/web/packages/saemix/index.html) in R. The saemix package implements the Stochastic Approximation Expectation Maximisation (SAEM) algorithm [START_REF] Comets | SAEMIX, an R version of the SAEM algorithm, Population Approach Group in Europe[END_REF], which computes the maximum likelihood estimator of the population parameters in NLMEM, without any approximation of the model. Then, for each vector of parameters Ψ1 , a design ξ 2 is optimized with the combined Fisher information matrix described by equation ( 15). 10 simulations are then performed for each of the 10 designs ξ 2 for the N 2 children of the second cohort. Thus, we obtain 100 data sets (Y 1,1 2 , ..., Y 10,10

2

), and we combine data Y 1 and Y 2 from which we estimate 100 vectors of parameters Ψ2 using saemix. The simulation scheme of adaptive clinical trials is described in Figure 3.

Evaluation methods

Our aims are to evaluate the impact of two-stage adaptive designs on the precision of parameter estimation on the one hand, and to investigate the influence of the sample size ratio of each stage on the other hand. To evaluate a design, we compute the relative root mean square errors (RRMSE) for the 100 estimated Ψ2 of the 2 one-stage designs, ξ 1 and ξ * , and of the 3 adaptive two-stage designs, ξ 50-10 , ξ 30-30 , and ξ 10-50 , as follows:

RRM SE(Ψ p ) = E(( Ψ2p -Ψ * p ) 2 ) × 100 Ψ * p (17)
with Ψ2p the estimated p component and Ψ * p the true value. We also calculate the standardized RRMSE for each parameter and each design as the RRMSE divided by the RRMSE obtained with ξ * , ξ * being the reference because it is the design optimized with the true parameters. For each design, mean standardized RRMSE is then computed.

We are also interested in the criterion predicted by PFIM, associated with the evaluation of one-stage designs. The criterion evaluates the amount of information and is computed with the true parameters Ψ * . For the two-stage designs, we compute a mean criterion, defined as the mean of the 10 criterion values obtained for each of the 10 two-stage designs computed with Ψ * . Indeed, the design ξ 1 is identical for all data sets at the first stage but we have 10 designs ξ 2 at the second stage, optimized according to the 10 data sets at the first stage.

It is of interest to note that we could have optimized the designs with respect to the RRMSE-criterion. However this is more time consuming as there is no expression of the expected RRMSE and hence the criterion has to be evaluated by simulation. This is why most approaches use the D-optimality criterion, which indeed assumes no bias and uses an approximation for the expected M F .

In section 5.1 we present the comparison of the RRMSE and of the criteria for the onestage (non adaptive) designs, and for the two-stage design ξ 30-30 , with the same number of children at the first and second stages. Then, we focus in Section 5.2 on the comparison of RRMSE and criteria of two-stage adaptive designs, with different sizes of cohort for the first and second stages (N 1 and N 2 respectively equal to 50 and 10, 30 and 30, 10 and 50).

RESULTS

Precision of parameter estimation

For the one-stage designs, the RRMSE and standardized RRMSE of each parameter, computed from the 100 vectors of estimated parameters Ψ2 are presented in Table 2. As expected, we can see that the RRMSE are lower for the design ξ * , which corresponds to the optimal design, and are high for the design ξ 1 , optimized from a priori parameters Ψ 0 .

These results are summarized by the mean standardized RRMSE, which is equal to 1.76 for ξ 1 , that is almost the double of the one for the optimal design ξ * which equals to 1 (reference value). Regarding the criterion for ξ 1 , computed with Ψ * , its value is 26.9, much lower than the criterion value for ξ * , equal to 84.7. It underlines the low information obtained with ξ 1 .

For the two-stage design with N 1 = N 2 = 30, the 10 optimal designs ξ 2 obtained for the 10 vectors of parameters Ψ1 are displayed in Figure 4 (notice that only 6 different designs are obtained). These designs ξ 2 are closer to ξ * than to ξ 1 and two of them coincide with ξ * . The RRMSE for this two-stage design ξ 30-30 , computed from the 100 vectors of parameters Ψ2 are close to the RRMSE obtained for the optimal design ξ * (see Table 2). For some parameters (for instance fixed effects CL and V2) the optimal design ξ * provides slightly greater RRMSE than the adaptive design, whereas it provides better RRMSE for other parameters and a better overall mean RRMSE. This is because a global optimization criterion on all parameters was used and also because the adaptive design gives very good results (mean standardized RRMSE of 1.06). Indeed, the mean standardized RRMSE of ξ 30-30 is equal to 1.06 and is therefore close to the value obtained for ξ * (equal to 1). The mean (respectively min and max) criterion values for the 10 two-stage design ξ 30-30 , evaluated in Ψ * , are 72.5

(respectively 69.1 and 74.6), which is close to the one for the optimal design. We can see that two-stage design has allowed a correction of the unsatisfactory results obtained for the non adaptive design ξ 1 .

Impact of size of each cohort

For the 3 two-stage designs with various number of children (ξ 50-10 , ξ 30-30 and ξ 50-10 ), the RRMSE and standardized RRMSE, computed from the 100 vectors of parameters Ψ2 , are presented in Table 3. The RRMSE are slightly worse for the design ξ 10-50 when the number of children for the first cohort is low. Indeed, its associated mean standardized RRMSE is equal to 1.15 whereas the mean standardized RRMSE for the two-stage adaptive designs ξ 30-30 and ξ 50-10 are closer to the one for the optimal design and are respectively equal to 1.06 and 1.07. Nevertheless, this trend is not confirmed by the mean (respectively min and max) criterion values since the later for the design ξ 10-50 equals 72.2 (respectively, 45.2 and 82.0) and is in the same range as the mean criterion obtained for the design ξ 30-30 . The mean (respectively min and max) criterion values obtained for the design ξ 50-10 are 59.0 (respectively 56.5 and 60.9). The difference between the minimum criterion and the maximum criterion for the design ξ 10-50 is important, contrary to the design ξ 30-30 . The minimum and maximum criterion values for the designs ξ 30-30 and ξ 50-10 are rather close, indicating that 30 subjects in the first cohort are enough to significantly reduce the variability of the design used at the second stage. Other investigations would be needed to confirm or infirm this observation.

DISCUSSION

We have proposed an approach for two-stage design for NLMEM (with an implementation that will be available soon in the next version PFIM 4.0). In the present work, the calculation of M F was performed using a first-order linearization of the model by the Taylor expansion. The derivative of the variance of observations according to the fixed effects was assumed to be zero, which corresponds to a block diagonal expression for M F . There is no clear consensus on what is the best approximation but it was suggested [START_REF] Mielke | Some considerations on the Fisher information in nonlinear mixed effects models[END_REF][START_REF] Mielke | Approximations of the Fisher information for the construction of efficient experimental designs in nonlinear mixed effects models[END_REF] that the block diagonal expression might be more reliable than one with the full M F . In the simulation study on two examples [START_REF] Nyberg | Methods and software tools for design evaluation in population pharmacokinetics-pharmacodynamics[END_REF], when using first order approximation, the block matrix gave better results. On the other hand, we did not take into account the adaptation in the optimization of the M F of the first step, nor that estimation was performed after first step for definition of M F of the second step. Hence, the dependence implied by the adaptation was neglected [START_REF] Pronzato | Design of experiments in nonlinear models[END_REF]Fedorov and Leonov, 2013).

Using simulations, we have studied the impact of two-stage designs on the precision of parameter estimation for a population PK study in children. Saemix in R was used for parameter estimation and evaluation and optimization were performed with PFIM in R. We compared the results obtained for non adaptive designs, that is with the total number of children N = 60 having the design ξ 1 , optimized with the wrong parameters, and with the total number of children having the design ξ * , optimal for the true parameters, to the results obtained for the two-stage adaptive design, with the same number of children at each of the two stages (N 1 = N 2 = 30). Results showed that the the non adaptive design ξ 1 provides poor precisions of parameter estimates. The adaptive two-stage design allows to compensate these poor results and yields performance close to those obtained for ξ * , optimized with the true parameters.

Then, we investigated the influence of the ratio of sample size of each stage. For that, we studied two-stage designs with a number of children for each stage equal to 50/10, 30/30 or 10/50, respectively. The RRMSE were sometimes less satisfactory for the third design whose the first cohort is small. The results obtained for two-stage adaptive designs were generally satisfactory for each of the three ratios considered and close to the ones obtained for the design ξ * .

Other studies would be needed and it would also be important to increase the number of replications in the simulation study which is rather low. Indeed, we estimated 100 vectors of parameters at the end of the second stage but only 10 at the end of the first stage. The reason why we have limited the number of replications is the absence of automatic connec-tion between saemix and PFIM, although both are written in R. Consequently, a manual intervention at the end of the first stage is needed to incorporate into PFIM the vectors of parameters Ψ1 estimated with saemix and which are needed for the optimization of designs ξ 2 . One perspective of this work would be to create an automatic connection to allow an important number of replications and for the use of this approach in clinical trials.

One limitation of this work was the assumptions made to conduct this simulation study. Indeed, we chose to fix the size of each cohort. Moreover, we imposed to have the same elementary design for all children in each cohort with a fixed number of sampling times.

Although these choices represent constraints usually needed in clinical trial, other scenario should be studied. Also, further simulation studies are needed to set the balance between size of first and second cohorts and with other examples. One theoretical study could also be performed to determine the balance between the sizes of first and second cohorts in NLMEM, as already done by [START_REF] Pronzato | Design of experiments in nonlinear models[END_REF] for generalized nonlinear models.

In this work, we have assumed that the model was known and correct. Indeed, the same a priori model at the beginning of the simulation study (for the first stage) was used for the second stage. No sensitivity analysis with respect to model misspecification was performed in the present work. Another step in adaptive design would be to use model averaging approaches [START_REF] Claeskens | Model selection and model averaging[END_REF].

In conclusion, local design, usually conducted for PK studies in children, relies on several strong hypotheses. Indeed, optimal design in NLMEM depends on the PK model and its parameters. Besides the structure, an a priori value should be given for each parameter of the model. We showed that the adaptive two-stage design is a promising approach for pediatric PK studies, relatively easy to conduct in clinical trials as there is only one interim analysis.

In the example considered, two-stage adaptive designs provided satisfactory parameter estimates, close to the ones obtained for the optimal design with the true parameters ans thus allowed a correction of the poor results obtained for the non adaptive design optimized with the wrong a priori parameters. In conclusion, although adaptive designs are still barely used in population PK, results are promising and should stimulate the conduction of PK studies with this approach and help convince pharmaceutical companies of its interest [START_REF] Mentré | Current use and developments needed for optimal design in pharmacometrics: a study performed amongst DDMoRe's EFPIA members[END_REF].

Figure 1 :

 1 Figure 1: Schematic representation of a two-stage adaptive design.

Figure 2 :

 2 Figure 2: Mean PK profiles in semi-log scale for the 2 sets of parameters and optimal design with 5 samples.

Figure 3 :

 3 Figure 3: Simulation plan of 100 two-stage adaptive clinical trials.

Figure 4 :

 4 Figure 4: The 10 second-stage designs ξ 2 optimized from the 10 estimated Ψ1 .

Table 1 :

 1 Population pharmacokinetic parameter valuesParameters Ψ 0 Ψ *

	ka (h -1 )	3.0 14
	CL (L.h -1 ) 1.5 1.0
	V1 (L)	2.0 1.0
	Q (L.h -1 ) 1.0 2.0
	V2 (L)	1.5 2.0

Table 2 :

 2 RRMSE computed for each parameter from the 100 vectors of parameters Ψ2 for each studied design RRMSE % (standardized RRMSE) RRMSE in bold are RRMSE greater than 30% for the fixed effects and greater than 50% for the random effects.

	Parameters Values	ξ 1 (60-0)	ξ 30-30 (30-30)	ξ * (0-60)
	ka (h -1 )	14	160 (7.05)	27.3 (1.20)	22.7
	CL (L.h -1 )	1	7.17 (1.07)	5.77 (0.857)	6.73
	V1 (L)	1	25.2 (1.65)	23.7 (1.55)	15.3
	Q (L.h -1 )	2	27.1 (1.67)	18.8 (1.16)	16.2
	V2 (L)	2	11.4 (1.00)	9.73 (0.854)	11.4
	ω 2 ka	0.3	100 (1.21)	90.4 (1.10)	82.4
	ω 2 CL	0.3	17.0 (0.950)	18.3 (1.02)	17.9
	ω 2 V 1	0.3	48.0 (1.33)	38.3 (1.06)	36.0
	ω 2 Q	0.3	89.9 (1.50)	71.5 (1.19)	59.9
	ω 2 V 2	0.3	36.4 (1.17) 26.6 (0.853 )	31.2
	σ slope	0.2	10.4 (0.765) 10.8 (0.794)	13.6
	Mean standardized RRMSE	1.76	1.06	1.00
	a

Table 3 :

 3 RRMSE computed for each parameter from the 100 vectors of parameters Ψ2 for each studied design RRMSE % (standardized RRMSE) RRMSE in bold are RRMSE greater than 30% for the fixed effects and greater than 50% for the random effects.

	Parameters Values	ξ 10-50	ξ 30-30	ξ 50-10
	ka (h -1 )	14	30.3 (1.33)	27.3 (1.20)	32.8 (1.44)
	CL (L.h -1 )	1	(1.05) 5.77 (0.857)	4.88 (0.725)
	V1 (L)	1	19.3 (1.26)	23.7 (1.55)	21.1 (1.38)
	Q (L.h -1 )	2	21.2 (1.31)	18.8 (1.16)	19.2 (1.19 )
	V2 (L)	2	12.3 (1.08) 9.73 (0.854)	11.4 (1.00)
	ω 2 ka	0.3	86.3 (1.05) 90.4 (1.10)	82.6 (1.00)
	ω 2 CL	0.3	30.6 (1.71)	18.3 (1.02)	21.4 (1.20)
	ω 2 V 1	0.3	23.5 (0.653) 38.3 (1.06)	37.1 (1.03)
	ω 2 Q	0.3	78.2 (1.31) 71.5 (1.19)	68.0 (1.14)
	ω 2 V 2	0.3	34.1 (1.09) 26.6 (0.853 )	34.3 (1.10)
	σ slope	0.2	11.0 (0.809) 10.8 (0.794)	8.14 (0.599)
	Mean standardized RRMSE	1.15	1.06	1.07
	a

ACKNOWLEDGMENTS

The research leading to these results has received support from the Innovative Medicines Initiative Joint Undertaking under grant agreement n • 115156, resources of which are composed of financial contributions from the European Union's Seventh Framework Programme (FP7/2007-2013) and EFPIA companies' in kind contribution. The DDMoRe project is also financially supported by contributions from Academic and SME partners.


Patrice Carle

Christine Choppy

Romain Kervarc

Ariane Piel

Handling Breakdowns in Unmanned Aircraft Systems

Keywords: safety monitoring, behaviour recognition tool, industrial applications, unmanned aircraft systems

This work is devoted to activity recognition in the setting of data analysis in aeronautics. Formal methods are applied to the certification and safety analysis processes of Unmanned Aircraft Systems in breakdown situations. The behaviour of these systems in case of a failure is entirely modeled and implemented. A temporal languagethe Chronicle language -describes arrangements of events which are employed to detail undesired circumstances that would lead to breaches in safety. A C++ chronicle recognition tool is used to recognise all the possible occurrences of these situations as soon as they occur.

Introduction

One obstacle to the insertion in controlled or uncontrolled airspace of aircrafts without pilots onboard results from security issues linked to the global consistency of the system. In the framework of operation safety analysis, we assess the possibility to detect incoherent states between the different entities making up the system. We formalise these different incoherent states in order to be able to automatically recognise them, and hence offer the opportunity of a self-acting surveillance. More specifically, this ongoing work relies on a fragment of the IDEAS project in charge of the Insertion of Unmanned Aircrafts in Airspace and Security, and tackles consistency problems in breakdown handling policies for Unmanned Aircrafts (UA).

An UA is defined by the FAA (Federal Aviation Administration) as "an aircraft that is operated without the possibility of direct human intervention from within or on the aircraft". The system required to operate safely an UA is referred to as an Unmanned Aircraft System (UAS). There exists several types of such architectures, and our model is based on the one presented in Fig. 1. It is composed of three entities, the UA, the Remote Pilot Station (RPS), and the Air Traffic Control (ATC). All three interact via several communication links. The RPS pilots the Unmanned Aircraft via Telecommand (TC), and the Unmanned Aircraft sends information to the pilot through Telemetry (TM). In addition, the ATC and the pilot can communicate via radio (Voice) relayed by the UA. Hence, the dynamic data flows between the agents of the system and between different systems if several UAS are considered are very elaborate. Moreover, each agent deduces from its own observations the state of the other agents. The situation in case of a fault can therefore be very complex. Formal methods are well-suited for monitoring UAS because these systems are highly automated as well as critical so require strong risk-free guarantees.

In the framework of the IDEAS project, the behaviour of each entity in case of a failure has been specified [START_REF] Lang | IDEAS-T1.1: Architecture de système de drone et scénarios de missions[END_REF]. Our work focuses on overseeing the consistency between the three entities during these rundown story lines.

Consider for instance the following situation: an unmanned aircraft which telecommand has been poorly working and then completely lost. The conforming procedure has been followed correctly, and the ATC has activated the urgency mode corresponding to the rerouting of the aircraft. However, the telecommand is suddenly restored. The pilot, busy redirecting the UA back to its original track, forgets to warn the ATC to cancel the urgency procedure. Thus, the ATC considers the aircraft to be rerouting itself to a given diversion airport while it is actually taking another route. It is important to be able to identify this situation which is critical since the ATC, lacking the appropriate information, might not properly organise the flow of traffic to separate aircrafts and prevent collisions.

A Preliminary Semi-Formal Representation

The aim of this work is to provide a tool for automatic monitoring of consistency in UAS. It is therefore necessary to formally set up the problem. The situation is standardised into the UML norm [START_REF]OMG Unified Modeling Language TM (OMG UML)[END_REF] as follows.

As a first step, a class diagram is built describing precisely the structure of the UA as shown in Fig. 2. Each UA X (UA(X)) is seen as made up of several components: the telecommand (UA(X) TC), the telemetry (UA(X) TM), the autopilot (UA(X) Pilot), the radio relay (UA(X) Relay), and the transponder code it sends to the ATC (UA(X) Code). Similarly, the associated RPS is composed of the telecommand (RPS(X) TC), the telemetry (RPS(X) TM) and a connection to the ATC (RPS(X) -ATC connection). Note that the smooth functioning of each end of a given communication link is considered separately so as to be able to differentiate the state of each entity to identify the source of potential problems and check for the consistency between the states. For example for the telecommand, the emission of instructions from the pilot, RPS(X) TC, and the reception of these instructions by the UA, UA(X) TC, are distinguished.

UA(X )

UA(X) T C UA(X) T M UA(X) pilo t UA(X ) r elay UA(X ) Co de

RPS(X ) RPS(X) T C RPS(X)-AT C connectio n RPS(X) T M ATC ATC -UA(X)

1 1 * 1 1 1 1 1 1 1 1 1 1
Fig. 2. Class diagram describing the structure of an UA.

Secondly, so as to exploit the codified behaviours specified in the IDEAS project, these were translated into state-transition diagrams for each fault. The first breakdown studied is TC failure presented in Fig. 3: the pilot receives information from the UA via telemetry but cannot send out orders to it.

With such diagrams, the system is entirely modeled. Its life cycle is mirrored by the changes in the active states of the diagram. In order to allow direct simulation, the UML diagram is implemented in C++ using the Meta State Machine (MSM) library [START_REF] Henry | MSM library of boost[END_REF] of boost (Version 1.48.0) which provides a straightforward way to define state machines. So as to simulate the life cycle of the system, scenarios activate the transitions of the diagram.

A Recognition Tool

In Sec. 2, the simulation of the system that has to be overseen is described. This section will now deal with the issue of how to put in place this supervision of the system, the aim being to detect during the simulation critical situations previously determined by experts.

A temporal language -the chronicle language [START_REF] Carle | Behaviour recognition using chronicles[END_REF] -is used to formally describe these system behaviours to oversee the consistency between the three entities of the UAS. It formally describes arrangements of events: a chronicle is either a single event, the conjunction of two chronicles, the disjunction of two chronicles, the sequence of two chronicles or the absence of a chronicle during another chronicle. For instance, let E, F and G be single events, chronicle (E F) -[G] corresponds to the sequence E F (event E followed by event F) without event G occurring in between E and F.

This formal language enables us to precisely detail interesting (in our case undesired) arrangements of events by writing chronicles. Different examples of applications of this concept may be seen in [START_REF] Carle | Detecting behaviours within HLA distributed simulations with added analysis components[END_REF], and [START_REF] Doherty | A temporal logic-based planning and execution monitoring framework for unmanned aircraft systems[END_REF] offers another approach dealing with a temporal logic for unmanned aircrafts, but in an execution, rather than safety, monitoring framework. Subsequently, once appropriate chronicles are established, all that is left to do in order to supervise the consistency of the system is to watch for potential recognitions of these chronicles in the event flow.

Table 1. Chronicles overseeing the consistency of the system in case of TC failure.

RPS rushed decision

( ((to_UA(X)_nominal_flight | | to_UA(X)_transitory_mode) T T ) -[to_UA(X)_rerouting_mode] ) & to_ATC_rerouting_mode_confirmed
where T T is a given delay.

The ATC is in rerouting mode confirmed but UA(X) is not in rerouting mode. ATC late

(to_TC_unrecovered T T ) -[to_Urgency_service || to_UA(X)_nominal_TC]
where T T is a given delay.

UA(X) TC changed to TC unrecovered and since then TT minutes have elapsed without the Urgency service of the ATC being an active state. ATC incoherent to_UA(X)_rerouting_mode

((to_ATC_rerouting_mode_inferred D) -[to_ATC_rerouting_mode_confirmed || to_UA(X)_nominal_flight])
where D is a given delay.

UA(X) changed to rerouting mode. Following that, the ATC inferred the rerouting mode but it was not confirmed by the pilot even after a delay of D minutes.

In our work, the set of events considered to build these chronicles are the entrances in the different possibly active states of the diagram of Fig. 3. Table 1 provides a few examples of chronicles overseeing the consistency of the system in the case of TC failure. The first chronicle "RPS rushed decision" corresponds to the example presented in Sect. 1. Chronicle "ATC late" could be triggered by the following circumstances. Consider an UA which has been periodically losing its TC link, for instance because of a disruptive environment, and the urgency transponder code associated to TC failure has been coming on and off. Both the pilot and the ATC are tired out and less alert, so when the TC is suddenly lost for good, no one reacts. After a significant delay predefined in its flight plan, the UA starts to reroute to a diversion airport. If neither the pilot nor the ATC realise the change in situation other aircrafts might not be sufficiently separated from the UA which cannot be sent orders to anymore. Recognising the chronicle would have prevented the possibly dangerous outcome.

For a given chronicle and a given stream of events, we are interested in establishing the list of all its recognitions gradually as events flow. This requirement prevents using simple finite state automata [START_REF] Bertrand | Chronicle modelling using automata and colored Petri nets[END_REF]. With this design, a recognition tool called Chronicle Recognition System (CRS/ONERA) has been developed by the ONERA in the late 1990s [START_REF] Carle | La reconnaissance d'intentions comme dynamique des organisations[END_REF]. It is designed on the basis of duplicating automata so as to comply with performance and interoperability requirements.

In the framework of this Ph.D. thesis, it is interesting to develop a new recognition tool that directly results from the set semantics of the chronicle language, since the recognitions produced by this tool are therefore adequate by construction and no adequacy proof is necessary. This tool, called Chronicle Recognition Library (CRL), is implemented in C++. Chronicles are plugged into the program, and then, gradually as events flow in, the program gives the set of all the recognitions of each chronicle, specifying for each recognition which events lead to it.

In this implementation, the events are represented as triplets (name, date, order) where date corresponds to the date of occurrence, and order is distinct for each event, allowing to sort and order the buffer of events to be processed. Indeed, two events can occur at the same date but their order is necessarily different. For instance, two sensors can observe the same system and hence generate events at the same instants but these will be received by the observing system in a given order.

(E || F)&G (E || F)&G # # G G G G G G G G G G y y s s s s s s s s s s (E || F) & & L L L L L L L L L L L
x x r r r r r r r r r r r A recognition of a given chronicle is a set of events leading to it. The recognition method will be illustrated by chronicle (E || F)&G on the buffer of events {(e, d e , 1), (h, d h , 2), (g, d g , 3), (f, d f , 4)}. The program processes events from the event buffer in their order. The recognition process is defined by induction for each class of chronicle. A key notion in this recognition system is that of subchronicles shown in Fig. 4.

The chronicle to be recognised is seen as a binary tree, the nodes and leaves of which are respectively operators and single events. The process computes two sets for every node -and hence for every sub-chronicle -after each event: the overall recognition set of the sub-chronicle, and its new recognitions set. In Table 2, both these sets -the new recognitions in green -are presented for each sub-chronicle and after each event of the buffer. The new recognitions set is emptied after each processing of an event. It is used to make sure to only add new recognitions to the recognition set -and hence avoid doubles -and to determine if there will be new recognitions in the ancestors of the node.

The recognition set of the chronicle is the recognition set of the root of the tree. In our example, after having processed all four events, there are two recognitions and the final recognition set is {{(e, d e , 1), (g, d g , 3)}, {(g, d g , 3), (f, d f , 4)}}. Note that event (g, d g , 3) is in both recognitions. The current work on the development of CRL is to take into account delay constraints, e.g. TT in chronicle ATC late. Mirroring Allen's interval logic [START_REF] Allen | Maintaining knowledge about temporal intervals[END_REF], delays are added to the chronicle language to be able to express duration of chronicles (lasts δ, at most δ, and at least δ -where δ is a delay) and delays after chronicles (then δ). It is clear that the system must be observed at each occurrence of an event. However, for chronicle C then δ, this is not enough and the question of when to examine the system is raised. Indeed, it cannot be constantly looked over. This is a general issue once continuous time and discrete events have to coexist. Different approaches have been put forward to address this problem, notably [START_REF] Bennett | A unifying semantics for time and events[END_REF].

A first necessary step for this ongoing work is therefore to define a function which, for both a given chronicle and the past events that have been observed, returns the next time at which the system has to be checked for a new recognition. This can be seen as a sort of validity contract, which can be directly linked to distributed simulations since these carry messages with validity domains. This question of validity is indeed strongly present in distributed simulations for aeronautics but also in certain formal methods for dynamic systems.

Conclusion and perspectives

In this paper, a temporal language is used to monitor safety for Unmanned Aircraft Systems in case of breakdowns, which is a problem of particularly great importance in the field of civil aviation. The behaviour of the UAS is completely modeled in a UML diagram which is implemented in C++, and a C++ library allows direct analysis of simulation data.

Chronicles and the associated recognition model shall be extended to improve the expressivity of the framework, in particular with delay-related constructs inspired from interval logics. The UAS application previously exposed tackles the crucial problem of the insertion of UAS in controlled or uncontrolled airspace and the associated certification and safety analysis processes. As presented in this paper, chronicle recognition can be used as a crucial tool to analyse simulation as the design stage of an UAS. It is also well-suited for the on-line analysis of real-time situations to generate alarms and avoid risky situations. This first step deals with physical failures in one aircraft but there are other issues, as separation and collision avoidance which are the cornerstone of flight safety. A model for this is currently being developed in collaboration with UAS engineers. Chronicles, already applied to a wide variety of different fields, are shown to be an adequate generic means to represent knowledge in a multi-agent system, and, as such, benefit to a large spectrum of applications. Moreover, the insertion of UAS into general airspace raises concerns that cannot be solved without a formal representation which chronicles seem to fulfil fairly.

Fig. 1 .

 1 Fig. 1. Architecture of an Unmanned Aircraft System.

Fig. 3 .

 3 Fig. 3. State-transition diagram describing the system behaviour in case of TC failure.

Fig. 4 .

 4 Fig. 4. Syntactic tree corresponding to chronicle (E || F)&G labelled with sub-chronicles.

Table 2 .

 2 Evolution of the recognition sets for chronicle (E || F)&G.

		(e, de, 1)	(h, d h , 2)	(g, dg, 3)	(f, d f , 4)
	E	{{(e, de, 1)}} {{(e, de, 1)}}	{{(e, de, 1)}} {}	{{(e, de, 1)}} {}	{{(e, de, 1)}} {}
	F	{} {}	{} {}	{} {}	{{(f, d f , 4)}} {{(f, d f , 4)}}
	E || F	{{(e, de, 1)}} {{(e, de, 1)}}	{{(e, de, 1)}} {}	{{(e, de, 1)}} {}	{{(e, de, 1)}, {(f, d f , 4)}} {{(f, d f , 4)}}
	G	{} {}	{} {}	{{(g, dg, 3)}} {{(g, dg, 3)}}	{{(g, dg, 3)}} {}
	(E || F) &G	{} {}	{} {}	{{(e, de, 1), (g, dg, 3)}} {{(e, de, 1), (g, dg, 3)}}	{{(e, de, 1), (g, dg, 3)}, {(g, dg, 3), (f, d f , 4)}} {{(g, dg, 3), (f, d f , 4)}}

Acknowledgements.

The authors would like to thank J. Bourrely for his help and support for this application, and T. Lang and C. Le Tallec for their useful insights on Unmanned Aircrafts.


Thomas Gustavsson

Dimitra Markovitsi

Ignacio Vayá

Paula Bonancía

M Consuelo Jiménez

Miguel A Miranda

Drug/protein interactions studied by time-resolved fluorescence spectroscopy

Keywords: drugs, flurbiprofen, albumin, time-resolved fluorescence, femtosecond, fluorescence upconversion

 of the interaction between flurbiprofen (FBP), a chiral non-steroidal anti-inflammatory drug, and human serum albumin (HSA), the main transport protein in the human body. We compare the results obtained for the drug-protein complex with those of various covalently linked flurbiprofentryptophan dyads having well-defined geometries. In all cases stereoselective dynamic fluorescence quenching is observed, varying greatly from one system to another. In addition, the fluorescence anisotropy decays also display a clear stereoselectivity. For the drug-protein complexes, this can be interpreted in terms of the protein microenvironment playing a significant role in the conformational relaxation of FBP, which is more restricted in the case of the (R)enantiomer.

INTRODUCTION

The transport and distribution of drugs in the body is a complex mechanisms governed by many different processes. An important role is played by specific transport proteins, responsible for the fixing and release of the drugs. Obtaining a better understanding of the factors that govern the structure and dynamics of drug-protein complexes is therefore of great interest.

Among the various techniques used to study the drug-protein interaction, fluorescence spectroscopy is very useful since it involves the first singlet excited state which is very sensitive to the local surroundings. In particular, time-resolved fluorescence informs on the excited state interactions occurring between the photo-excited drug and the protein, which in turn provides information regarding the structural dynamics of the complex.

Several femtosecond fluorescence spectroscopic measurements of drug-protein interactions have been reported in the past. [START_REF] Zhong | Femtosecond studies of protein-ligand hydrophobic binding and dynamics: Human serum albumin[END_REF][START_REF] Zhong | Femtosecond dynamics of a drug-protein complex: Daunomycin with Apo riboflavin-binding protein[END_REF][START_REF] Kamal | Ultrafast hydration dynamics in protein unfolding: Human serum albumin[END_REF][START_REF] Douhal | Femtochemistry of orange II in solution and in chemical and biological nanocavities[END_REF][START_REF] Cohen | Proton-Transfer Reaction Dynamics within the Human Serum Albumin Protein[END_REF]. These have shown the influence of the biological nanocavities not only on the diffusional/orientational motions of the drug, but also on its solvation dynamics and excited state processes such as energy, electron and proton transfer.

We report here on the interaction between flurbiprofen (FBP), a chiral non-steroidal anti-inflammatory drug, and human serum albumin (HSA), the main transport protein in the human body, using time-resolved fluorescence spectroscopy. To this aim, we have used a combination of two time-resolved fluorescence techniques; fluorescence upconversion (FU), with a time-resolution of about 100 femtoseconds, and time-correlated single photon counting (TCSPC), with a timeresolution of about 25 ps.

HSA is the most abundant protein in the bloodstream, responsible for the transport of different agents such as drugs. Therefore, the binding of a drug to HSA constitutes a key process, relevant for its distribution. It is known that FBP binds to HSA preferentially in the so-called site II. HSA contains only one tryptophan (W214), situated close to site II. Both FBP and W214 absorb and emit light in the UV spectral region, rendering the analysis of data delicate. FBP can also interact with W214 via energy-or electron-transfer after photo-excitation.

In order to investigate the FBP-Trp interaction in more detail, two covalently linked dyads, (S,S)-or (R,S)-, formed by either (S)-or (R)-FBP and (S)-TrpMe ((S)-tryptophan methyl ester) have been synthesized and employed as model systems (Figure 1) [START_REF] Vaya | Excited-State Interactions in Flurbiprofen-Tryptophan Dyads[END_REF]. Due to their low solubility in water, experiments were performed in acetonitrile.

Such drug-amino acid dyads have already been designed and studied with success in the past [START_REF] Abad | Wavelength-dependent stereodifferentiation in the fluorescence quenching of asymmetric naphthalene-based dyads by amines[END_REF][START_REF] Abad | Diastereodifferentiation of novel naphthalene dyads by fluorescence quenching and excimer formation[END_REF]. They have provided new mechanistic insight into the fundamental processes that occur between the two chromophores such as energy transfer, electron transfer, exciplex formation, etc.

STEADY-STATE FLUORESCENCE

Previous steady-state spectroscopic measurements of the dyads have shown that a strongly stereo-selective fluorescence quenching (>90%) occurs and that only a weak residual TrpMe emission can be observed (Figure 2A). This has been taken as an indication of efficient FBP -> Trp energy transfer [START_REF] Zhong | Femtosecond studies of protein-ligand hydrophobic binding and dynamics: Human serum albumin[END_REF]. For both dyads, a band centered at 450 nm, assigned to exciplex emission, was also observed (dashed lines in Figure 2A). Strong fluorescence quenching occurs also in the FBP/HSA complex, even though less important and with less marked stereo-differentiation than for the dyads (Figure 2B).

TIME-RESOLVED FLUORESCENCE OF THE DYADS

Our results show that the fluorescence decays of both dyads are much faster than those of free FBP or TrpMe. At the fluorescence maximum of FBP, around 310 nm, the fluorescence decays of the two dyads recorded by FU occur on the picosecond time scale (Figure 3A). These results are in accordance with a dynamical fluorescence quenching. A significant stereo-selectivity can be observed; the (R,S)-dyad emission decaying faster than that of the (S,S)-dyad. It should be noted that no rise in the Trp fluorescence at 340 nm could be observed, which would have supported the idea of energy transfer. The fluorescence decays of the two dyads recorded by TCSPC at 450 nm are much slower than that of TrpMe and are characterized by a rise (Figure 3B), which can be assigned to the formation of the exciplex. Since exciplex formation can be related to charge transfer, this suggests that the fluorescence quenching mechanism also involves electron transfer. Interestingly, both the formation and the decay of the exciplex emission show a clear stereo-differentiation.

TIME-RESOLVED FLUORESCENCE OF THE FBP-HSA COMPLEX

The time-resolved measurements by FU show that the time evolutions of the fluorescence of (S)-FBP/HSA and (R)-FBP/HSA are similar at early times (Figure 4A). Only after several tens of picoseconds do the traces start to deviate. Curiously, the stereo-differentiation is inversed with regards to that observed for the dyads, with (S)-FBP/HSA showing a faster decay than (R)-FBP/HSA. There is also a marked stereo-selectivity in the fluorescence anisotropy decays recorded by TCSPC at 310 nm (Figure 4B). After a fairly rapid decay (hundreds of picoseconds) the anisotropies attain a constant value which persists over several nanoseconds. Interestingly, the fluorescence anisotropy decay of (S)-FBP/HSA is faster (0.44 ± 0.03 ns) and of higher amplitude than that of (R)-FBP/HSA (0.62 ± 0.07 ns).

The constant value at longer times implies a restricted reorientational motion, which can be described by the wobblingin-a-cone model [START_REF] Schroder | Simulation of fluorescence anisotropy experiments: Probing protein dynamics[END_REF]. The different reorientational dynamics observed for the two enantiomers indicate that the conformational relaxation is more restricted for (R)-FBP within the binding pocket. This could possibly be related to the pharmaceutical activity of flurbiprofen.

CONCLUSIONS

The photophysical behaviour of the drug flurbiprofen when covalently linked to Trp in model dyads has been compared with that of its non-covalent complex with human serum albumin.

A very fast dynamic and strongly stereo-selective quenching of the FBP fluorescence is observed in the dyads. This can be attributed to efficient energy and/or electron transfer. The residual Trp fluorescence also shows a stereo-selective quenching, although slower, together with exciplex formation.

Similar trends were observed in the drug-protein complexes, although the kinetics of the involved processes are slower and the stereo-differentiation is opposite to that observed for the dyads. This shows that the relative spatial arrangement of the chromophores constitutes a critical factor in the stereodifferentiation. Moreover, the fluorescence anisotropy clearly showed that the protein microenvironment plays a significant role in the conformational relaxation of FBP, which is more restricted in the case of the (R)-enantiomer. This stereo-selectivity is possibly related to the modes of drug binding to the protein, a process of pharmacological relevance.

The present results show that time-resolved fluorescence spectroscopy constitutes a powerful tool to study drug-protein interactions. The combination of two fluorescence techniques, fluorescence upconversion and time-correlated single photon counting, allows to cover a very wide time-window and thereby distinguishing different phenomena. Finally, to our knowledge, this is the first observation of stereo-selective drug-protein interaction dynamics recorded with femtosecond time resolution.

Figure 1 .

 1 Figure 1. The model dyads (S,S)-FBP-TrpMe (A) and (R,S)-FBP-TrpMe (B).

Figure 2 .

 2 Figure 2. Steady-state fluorescence spectra of (A) (S)-FBP (black), (S,S)-FBP-TrpMe (red), and (R,S)-FBP-TrpMe (green) in acetonitrile. Dashed lines correspond to the normalized emission spectra. (B) (S)-FBP (black), (S)-FBP-HSA (red) and(R)-FBP-HSA (green) in PBS, under air, using isoabsorptive solutions at the excitation wavelength (267 nm).

Figure 3 .

 3 Figure 3. Fluorescence decays at (A) 310 nm using FU and (B) 450 nm using TCSPC of (S)-FBP (black), (S,S)-FBP-TrpMe (red), (R,S)-FBP-TrpMe (green) and (S)-TrpMe (blue) in MeOH.

Figure 4 .

 4 Figure 4. (A) Fluorescence decays (by FU) and (B) fluorescence anisotropy decays (by TCSPC) at 310 nm of (S)-FBP (black), (S)-FBP/HSA (red) and (R)-FBP/HSA (green). (Fitted curves indicated by dashed lines)