Keywords: nonlinear mixed effects models, multiple responses, Fisher information matrix, population design, first order approximation, PFIM

SUMMARY

We focus on the Fisher information matrix used for design evaluation and optimization in nonlinear mixed effects multiple response models. We evaluate the appropriateness of its expression computed by linearization as proposed for a single response model. Using a pharmacokinetic-pharmacodynamic (PKPD) example, we first compare the computation of the Fisher information matrix by approximation to one derived from the observed matrix on a large simulation using the stochastic approximation expectation-maximization algorithm (SAEM). The expression of the Fisher information matrix for multiple responses is also evaluated by comparison to empirical information obtained through a replicated simulation study using the first order linearization estimation methods implemented in the NONMEM software (FO, FOCE) and the SAEM algorithm in the MONOLIX software. The predicted errors given by the approximated information matrix are close to those given by the information matrix obtained without linearization using SAEM and to the empirical ones obtained with FOCE and SAEM. The simulation study also illustrates the accuracy of both FOCE and SAEM estimation algorithms when jointly modelling multiple responses and the major limitations of the FO method. This study highlights the appropriateness of the approximated Fisher information matrix for multiple responses, which is implemented in PFIM 3.0, an extension of the R function PFIM dedicated to design evaluation and optimization. It also emphasizes the use of this computing tool for designing population multiple response studies, as for instance in PKPD studies or in PK studies including the modelling of the PK of a drug and its active metabolite.

Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model.

Caroline Bazzoli, Sylvie Retout, France Mentré

To cite this version:

Caroline Bazzoli, Sylvie Retout, France Mentré. Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model.. Statistics in Medicine, 2009, 28 [START_REF] Jonsson | Comparison of some practical sampling strategies for population pharmacokinetics studies[END_REF], pp.1940-56. [START_REF] Samson | Extension of the SAEM algorithm to leftcensored data in nonlinear mixed-effects model: application to HIV dynamics model[END_REF].1002/sim.3573. inserm-00371363

Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model Email: caroline.bazzoli@inserm.fr

Introduction

Nonlinear mixed effects models (NLMEM) are widely used to analyze various biological processes described by longitudinal data. Since the primary models developed by Sheiner et al. [START_REF] Sheiner | Population modelling in drug development[END_REF] in pharmacokinetic (PK) and pharmacodynamic (PD), NLMEM are become widely used for modelling of biological processes. NLMEM, also called the population approach, allow estimation of the mean value of the parameters in the studied population and their interindividual variability, or population characteristics. NLMEM are also now commonly used for the joint modelling of several biological responses such as the PK of parent drugs and of their active metabolite. NLMEM allow a sparse sampling design with few data points per individual in a large set of individuals. This can be particularly useful in studies in specific populations such as children or patients with serious diseases, where classical studies with a large number of samples are often limited for ethical or physiological reasons.

Estimation of the parameters in NLMEM is commonly performed by maximum likelihood.

However, due to the nonlinearity of the regression function, an analytical expression of the log-likelihood in nonlinear mixed effects models cannot be provided. To solve this issue several methods for estimating the parameters have been proposed, based on an approximation of the log-likelihood such as the First Order method (FO) or the First Order Conditional Estimate (FOCE) method proposed by Linsdstrom and Bayes [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF]. Both methods use a linearization of the structural model either around the expectation of the random effects parameter (FO) or around individual estimates of the random effects (FOCE). These methods have been implemented in the NONMEM software [START_REF] Beal | NONMEM users guide[END_REF][START_REF] Pillai | Non-linear mixed effects modeling -from methodology and software development to driving implementation in drug development science[END_REF] but also in the nlme function of Splus and R software [START_REF] Pinheiro | Mixed-Effects Models in S and S-Plus[END_REF]. Compared to FO, the FOCE method provides less biased estimates and, in the context of joint modelling of multiple responses, is more appropriate with fewer problems of convergence or of inter-individual variance estimation [START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data I: best-case performance[END_REF][START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data II: robustness of methods[END_REF]. Alternative methods have also been proposed to maximize the likelihood using a stochastic approximation of the integrals, such as the Gaussian quadrature [START_REF] Pinheiro | Approximations to the log-likelihood function in the nonlinear mixed effects models[END_REF] or the Adaptative Gaussian quadrature methods implemented in the NLMIXED procedure of SAS. Recently, the Stochastic Approximation Expectation-Maximization algorithm (SAEM) has been developed and implemented in the MONOLIX software [9,[START_REF] Samson | Extension of the SAEM algorithm to leftcensored data in nonlinear mixed-effects model: application to HIV dynamics model[END_REF]. It uses a stochastic approximation version of the standard expectation-maximization (EM) algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF][START_REF] Khun | Maximum likelihood estimation in nonlinear mixed effects models[END_REF].

The convergence and the consistence of the estimates have been proved by the authors. In this algorithm, the EM algorithm is used for finding maximum likelihood estimates of parameters in models, where the model depends on unobserved variables corresponding to the random effects in the NLMEM.

An appropriate choice of experimental design for estimating parameters in NLMEM is required. Called a population design in this framework, a design is defined as a group of elementary designs; each elementary design is composed of a set of sampling times to be performed in several individuals. Determining a population design involves identifying both the allocation of the sampling times and the whole group structure, that is to say the number of elementary designs, the number of samples per elementary design and the proportion or the number of individuals in each elementary design according to a fixed total number of samples. Simulation studies have shown that the precision of estimation of the parameters depends on the choice of the design [START_REF] Al-Banna | Experimental design and efficient parameter estimation in population pharmacokinetics[END_REF][START_REF] Jonsson | Comparison of some practical sampling strategies for population pharmacokinetics studies[END_REF] and that an appropriate choice can thus substantially improve the efficiency of studies. In the context of NLMEM with sparse designs, the challenge is then to determine the trade-off between few sampling times and informative data to obtain correct parameter estimates.

To evaluate population designs, the theory of optimum experimental design described for instance by Atkinson and Donev [START_REF] Atkinson | Optimum Experimental Designs[END_REF] or by Walter and Pronzato [START_REF] Walter | Identification of Parametric Models from experimental Data[END_REF] in classical nonlinear models, has been extended to NLMEM. This theory uses criteria based on the Fisher information matrix (M F ). It comes from the Cramer-Rao inequality; indeed, the inverse of M F is the lower bound of the variance covariance matrix of any unbiased estimators of the parameters. As the likelihood has no closed form in our framework, a linearization of the model around the expectation of the random effects has been proposed by Mentré et al. [START_REF] Mentré | Optimal design in random effect regression models[END_REF] and extended by Retout et al. [START_REF] Retout | Fisher information matrix for non-linear mixedeffects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[END_REF] to derive an approximate expression of M F . Accuracy of this approximation was first shown by simulation of an example based on a real PK study [START_REF] Retout | Fisher information matrix for non-linear mixedeffects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics[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 was confirmed by comparison of the predicted SE computed from this approximate M F to those given by an evaluation of M F without linearization obtained by stochastic approximation using the SAEM algorithm of MONOLIX [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]. The approximated expression of M F has been implemented in R functions PFIM and PFIMOPT for population design evaluation and optimization, respectively [21][START_REF] Retout | Development and implementation of the population Fisher information matrix for the evaluation of population pharmacokinetic designs[END_REF][START_REF] Retout | Optimization of individual and population designs using Splus[END_REF]. Recently, PFIM Interface 2.1, a graphical user interface version, has been developed, allowing both evaluation and optimization in the same tool [21]. However, currently, these tools only allow evaluation and optimization of population designs of single response models. For multiple response models, the same linearization method around the expectation of the random effects as for single response models has been proposed to approximate the population M F [START_REF] Gueorguieva | Optimal design for multivariate response pharmacokinetic models[END_REF][START_REF] Hooker | Simultaneous population optimal design for pharmacokinetic-pharmacodynamic experiments[END_REF][START_REF] Ogungbenro | Incorporating correlation in interindividual variability for the optimal design of multiresponse pharmacokinetic experiments[END_REF][START_REF] Ogungbenro | Optimal design for multiresponse pharmacokineticpharmacodynamic models -dealing with unbalanced designs[END_REF]. In those papers, illustrations of this development were provided using either a PKPD model or a joint PK model of a drug and its metabolite. However, the accuracy of the development of M F by linearization for multiple responses has not yet been evaluated. Even if the same linearization as in the single response is used, computation can become more complicated for multiple responses. Indeed, some parameters can be included in several responses and the information on those parameters is therefore obtained from each of those response profiles. This is usual in the PKPD context where PD response depends on the PK parameters. Moreover, as noted previously, use of the linearization around the expectation of the random effects appears to be inadequate for joint estimation of multiple response models [START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data I: best-case performance[END_REF][START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data II: robustness of methods[END_REF]. The appropriateness of its use in the context of design evaluation is thus also questionable and should be investigated.

The objective of this study was therefore to evaluate the first order approximation to compute the Fisher information matrix in NLMEM with multiple responses. To do this, we considered a PKPD simulation example associated with a population design. Then, we compared the predicted standard errors (SE), computed from the approximated expression of M F to those given by the evaluation of M F without linearization obtained by stochastic approximation using the SAEM algorithm of MONOLIX. We also performed another evaluation by comparison of those predicted SE to the empirical ones, obtained by estimation on simulated datasets using three different estimation algorithms: FO and FOCE (with NONMEM); SAEM (with MONOLIX). Based on those simulations, we also compared the performance of those three estimation methods in the same simultaneous analysis of this PKPD model.

In Section 2, we introduce the notations, describe the PKPD example and present the methodology used to evaluate M F and to compare the estimation methods. Section 3 describes the results of the evaluation and the comparison. Discussion of the results is provided in Section 4. The development of M F for multiple responses is given in detail in the Appendix.

Methods

Notation

In the nonlinear mixed effect multiple response model, an "elementary" design i ξ for one individual i is defined by i n sampling times. It is composed of several sub-designs such that

( ) 1 2
, , ,

i i i iK ξ ξ ξ ξ = K
, with ik ξ being the sub-design associated with the th k response,

1, , k K = K . ik ξ is defined by ( ) 1 2
, , ,

ik ik ik ikn t t t K
, the vector of the ik n sampling times for the observations of the th k response, so that

1 K i ik k n n = = ∑ .
For N individuals, we define a "population design" composed of the N allocated elementary

designs i ξ , 1, , = K i N .
A population design is therefore described by the N elementary designs for a total number n of observations such that

1 N i i n n = = ∑ : { } 1 , , Ξ = K N ξ ξ (1)
Usually population designs are composed of a limited number Q of groups of individuals with identical design within each group. Each of these groups is defined by an elementary design q ξ , 1, , = K q Q , which is composed, for the th k response, of qk n sampling times

( ) 1 2
, , , qk qk qk qkn t t t K to be performed in a number q N of individuals. The population design can then be written as follows:

[ ] [ ] { } 1 1 2 2 , ; , ; ; , Q Q N N N ξ ξ ξ   Ξ =   K (2)
A nonlinear mixed effects multiple response model or a multiple response population model is defined as follows. The vector of observations i Y for the i th individual is defined as the vector of the K different responses:

1 1
, , ,

T T T T i i i iK Y y y y   =   K (3)
where ik y , 1, , k K = K , is the vector of observations for the k th response. Each of these responses is associated with a known function f k which defines the nonlinear structural model. The K functions f k can be grouped in a vector of multiple response models F, such as:

( ) ( ) ( ) 1 1 2 2 ( , ) , , , , , ,
T T T T i i i i i i K i iK F f f f θ ξ θ ξ θ ξ θ ξ   =     K (4)
where i θ is the vector of all the individual parameters needed for all the response models in individual i. The vector of individual parameters i θ depends on β, the p-vector of the fixed effects parameters and on i b the vector of the p random effects for individual i. The relation between i θ and ( )

, i b β is modelled by a function g , ( )
, i i g b θ β = , which is usually additive, so that = + i i b θ β
, or exponential so that ( )

exp = i i b θ β . It is assumed that i b ∼ ( ) 0, N Ω with Ω defined as a p p × -diagonal matrix, for which, each diagonal element 2 r ω , 1, , r p = K ,
represents the variance of the th r component of the vector b i .

The statistical model is thus given by: ( ) ( )

, , i i i i Y F g b β ξ ε = + (5)
where i ε is the vector composed of the K vectors of residual errors ik ε , 1, , k K = K , associated with the K responses. We also suppose ik ε ∼ ( )

0, ik N Σ with ik Σ a ik ik n n × -diagonal matrix such that ( ) ( ) ( ) ( ) 2 int int , , , , , , ik i erk slopek ik erk slopek k i ik b diag f g b β σ σ ξ σ σ β ξ Σ = + (6)
i i er slope i b β σ σ ξ Σ the variance of i ε , over the K responses, such that i Σ is a i i n n × - diagonal matrix composed of each diagonal element of ik Σ with 1, , k K = K . slope σ and int er σ are two vectors of the K components int erk σ and slopek σ , 1, , , k K = K
respectively. Finally, conditionally on the value of b i , we assume that the errors i ε are independently distributed.

Let Ψ be the vector of population parameters to be estimated such as ( , , , , )

T T T p er slope λ ω ω σ σ = K , so that ( ) , T T T β λ Ψ = .

Population Fisher information matrix for multiple response models

The population Fisher information matrix for a population design Ξ (see Equation ( 1)), is defined as the sum of the N elementary Fisher information matrices ( )

, F i M ξ Ψ for each individual i: ( ) ( ) 1 , , N F F i i M M ξ = Ψ Ξ = Ψ ∑ (7)
In the case of a limited number Q of groups, as in Equation ( 2), it is expressed by:

( ) ( ) 1 , , Q F q F q q M N M ξ = Ψ Ξ = Ψ ∑ (8)
The expression of an elementary Fisher information matrix for multiple responses has been extended by Hooker et al. [START_REF] Hooker | Simultaneous population optimal design for pharmacokinetic-pharmacodynamic experiments[END_REF] using the same development as for single response models with a first Taylor expansion of the model as in Mentré et al. [START_REF] Mentré | Optimal design in random effect regression models[END_REF] and Retout et al. [START_REF] Retout | Development and implementation of the population Fisher information matrix for the evaluation of population pharmacokinetic designs[END_REF].

Its expression is given below for one individual i and depends on the approximated marginal expectation i E and variance i V of the observations i Y [START_REF] Gagnon | Optimal population designs for PK models with serial sampling[END_REF]. Details of this development are given in the Appendix.

( )

1 1 1 1 , 2
T i i i i F i i m l m l E E V V M V tr V V ξ - - -   ∂ ∂ ∂ ∂ Ψ = +   ∂Ψ ∂Ψ ∂Ψ ∂Ψ   with m and ( ) 1, , dim l = Ψ K . (9)
with, ( ) ( ( , 0), )

i i i E Y E F g β ξ ≅ = (10) ( ) ( ) ( ) ( ) ( ) ( ) int
, 0 , , 0 , , 0, , ,

T i i i i er slope i T F g F g Var Y V b b β ξ β ξ β σ σ ξ     ∂ ∂ ≅ = Ω + Σ         ∂ ∂     (11)
This expression has been implemented in an extension of the R function PFIM, PFIM 3.0.

This function has been developed for R 2.

PKPD simulation example

In this paper, we use a simple and typical PKPD model as an example to evaluate M F by simulation. It is derived from the one used by Hooker et al. [START_REF] Hooker | Simultaneous population optimal design for pharmacokinetic-pharmacodynamic experiments[END_REF] to illustrate the development of the Fisher information matrix for a multiple response model. The PK model for drug concentration is a one compartment with bolus input and first order elimination given as follows for the sampling time t PK :

exp PK PK PK PK C C dose Cl f (θ , t ) ( t ) V V = - (12)
where ( )

, T PK C Cl V θ =
is the vector of the PK parameters with Cl and V C , the clearance and the volume in the central compartment, respectively.

The PD model for drug effect is a simple max E model with baseline, expressed as a function of the predicted concentrations PK f , and given as follows for the sampling times PD t :

( ) ( ) max 0 50 , , , ( , )
PK PK PD PD PK PD PD PK PK PD E f t f t E C f t θ θ θ θ = + + (13)
where ( )

0 max 50
, ,

T PD E E C θ =
is the vector of the PD parameters with 0 E , max E and 50 C , the effect at baseline, the maximum effect and the concentration needed to observe half of the maximum effect, respectively.

We assumed an exponential model of the random effects for both the PK and the PD

( ) { } , , PK PD N ξ ξ   Ξ =   .
The curve profiles of the PK and the PD model for the fixed effects are displayed in Figure 1; the sampling times for each response are overlaid.

Evaluation of M F for multiple responses

Comparison of M F with and without linearization

In this section, we propose to compare the predicted SE obtained from the approximate M F for multiple responses computed by PFIM 3.0 to the SE obtained from more "exact" approaches using the SAEM estimation algorithm. This latter algorithm was used by Retout et al. [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] and Samson et al. [START_REF] Samson | The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed-effects model[END_REF] to show the appropriateness of this approximation in a single response model. This SAEM algorithm allows the observed population Fisher information matrix to be computed according to two approaches. The first approach was developed by Samson et al.

[29] and has been used to evaluate an "exact" population Fisher information matrix using the Louis's principle [START_REF] Louis | Finding the observed information matrix when using the EM algorithm[END_REF]. It does not require any linearization and can thus be considered as the "true" population Fisher information matrix. The second approach evaluates the Fisher information matrix using a linearization of the model around the conditional expectation of the individual parameters previously estimated by SAEM without any linearization.

To perform this comparison, we first computed the predicted M F for the population design associated with the PKPD example using PFIM 3.0, based on the linearization. We then simulated a dataset of PK and PD observations for 10 000 individuals in order to achieve asymptotic properties using the software R 2.4.1. To do that, we used the parameter values given in Table 1 and the sampling times shown in Figure 1, defining the PKPD example (section 2.3). For each individual i, we simulated a vector of random effects i b in ( )

0, N Ω ,
where the diagonal elements of Ω are the variance of the random effects, and we calculated the individual parameters using For estimation with the SAEM algorithm, we used an initial set of parameters with the values ( ) 0.2, 0.05, 1.2, 5,1.5 for the fixed effects, ( )

1, 1, 1, 1, 0.5 for the variance of the random effects and ( ) 0.5, 0.5 for the residual errors. The default values for the algorithm were used except for the number of Markov chains, which was set to 4, and the number of iterations with two different steps sizes, which was set to 1000 and 1000 to ensure good convergence.

The predicted SE obtained by linearization with PFIM 3.0 were designated PFIM. The notations SAEM_LO and SAEM_LI denote the predicted SE obtained with the SAEM algorithm using the Louis's principle and the linearization method, respectively.

Comparison of M F to empirical information through replicated simulation

Another objective was to compare the predicted SE of M F computed from PFIM 3.0 to the empirical SE obtained by the FO method, the FOCE method and the SAEM algorithm on simulated datasets. To do that, we simulated 1000 datasets of 100 individuals with the software R 2.4.1 using the same PKPD model and population design described previously.

Datasets were simulated using a similar method as in section 2.4.1, using the same parameter values and the same sampling times.

For each simulated data file, we estimated the population parameters for the PKPD model using first the FO method and the FOCE with interaction method implemented in NONMEM software version V and then, using the SAEM algorithm in MONOLIX (Version 2.1).

For the estimations using the FO and FOCE methods, two sets of initial parameters were defined. The first corresponded to the value of the parameters used for the simulation (Table 1). The second one was used only in the case of lack of convergence with the first set. The values of the second set of initial parameters were for the fixed effects: ( ) For each parameter of the PKPD model, we compared the predicted SE using the three evaluations of M F with PFIM, SAEM_LO and SAEM_LI, to the empirical SE obtained with the FO method, the FOCE method and the SAEM algorithm, denoted FO, FOCE, and SAEM, respectively. These empirical SE are defined as the sample estimate of the standard deviation from the parameter estimates for each method, considering only the subset of datasets fulfilling all convergence conditions.

0
We were also interested in comparing the distribution of the observed SE provided by each of the estimation methods to the empirical SE and to the predicted SE. In this case, we considered only the subset of datasets for which both the convergence and the variancecovariance matrix of estimation were obtained. For the distribution of the SE provided by the SAEM algorithm, we considered both methods of computation of the SE, the Louis's principle and the linearization.

Comparison of results for estimation methods with and without linearization

Using the previous simulations, we also compared the three methods of estimation: FO, FOCE and the SAEM algorithm. For each parameter, the relative bias as well as the relative RMSE were computed for the S datasets fulfilling convergence conditions ( )

1000 S ≤
, which, for l Ψ , the th l parameter of the population vector Ψ , are given by: ( )

0 0 1 1 s S l l l s l Bias S =   Ψ -Ψ Ψ =     Ψ   ∑ (14)
( )

2 0 0 1 1 s S l l l s l RMSE S =   Ψ -Ψ Ψ =     Ψ   ∑ ( 15
)
with ˆs l Ψ the estimated value of l Ψ for the th s simulated datasets and 0 l Ψ the true value.

Results

Comparison of M F with and without linearization

The SE predicted through the use of the SAEM algorithm on a large dataset (SAEM_LI and SAEM_LO) and those predicted by PFIM 3.0 are reported in This evaluation shows the appropriateness of the extension of the population Fisher Information matrix for multiple response models using the first order approximation.

Comparison of M F to empirical information through replicated simulation

Convergence was achieved for all datasets and the variance-covariance estimates were obtained for 997 datasets among the 1000 simulated datasets with the FO method.

Convergence with the FOCE method of NONMEM was obtained for only 853 sets. Among those 853 sets, the variance-covariance matrix of estimation was obtained for only 798 sets.

Finally, with the SAEM procedure, no problem of convergence or covariance-variance matrix was noted for any of the 1000 datasets.

For each parameter the empirical RSE obtained with the three estimation methods and the predicted RSE of SAEM_LI, SAEM_LO and PFIM are displayed in Figure 2. Concerning the FO method, the empirical RSE were much larger than the RSE of PFIM, SAEM_LI and SAEM_LO except for the PK parameters. This difference is above all important for the PD The distribution of observed RSE from the three estimation methods, including both methods of computation of M F with the SAEM algorithm, are reported as boxplots in Figure 3, part (A)

and part (B), for the mean and the variance parameters, respectively. For the FO method, the range of the observed RSE was much larger than for those obtained with the FOCE method or both SAEM procedures. However, the observed RSE of FO were concordant with the empirical ones. For all the parameters, the RSE predicted with PFIM were consistent with the distribution of the RSE observed with FOCE and the two SAEM procedures but not for FO.

The range of the observed RSE for FOCE, SAEM and the corresponding empirical RSE were also concordant. However, for most parameters, the RSE computed using the Louis's principle of SAEM had a broader distribution than by using linearization, with values for several datasets being outliers (Figure 3). This problem occurred in particular for the RSE on 50 2 C ω parameter.

In this example, the RSE predicted by PFIM, computed by the first order linearization, were thus concordant with the empirical ones and the observed RSE obtained from the simulation study.

Comparison of three estimation methods

The relative bias and relative RMSE obtained with the three estimation methods are presented in Table 3. Convergence was not achieved for 15% of the simulated datasets using the FOCE method whereas the FO method and the SAEM algorithm converged for all datasets.

Regarding the FO method, bias and RMSE were large especially for the parameters of the PD model (fixed effects, random effects and residual errors) whereas FOCE and SAEM provided reasonable bias and RMSE for all the parameters. For the fixed effects, slightly lower bias and RMSE were observed for the SAEM procedure compared to FOCE. We observed important RMSE ( )

40% > for the parameter 50 2 C
ω whatever the estimation method. This is in agreement with the large RSE obtained for that parameter previously.

Discussion

We evaluated the expression of the population Fisher information matrix for multiple response models using a linearization of the model [START_REF] Hooker | Simultaneous population optimal design for pharmacokinetic-pharmacodynamic experiments[END_REF], as for single response models. Note that our evaluation focused on the case of multiple responses, in which some parameters are involved in several responses. For cases in which parameters differ across responses, the same information would be obtained using M F either for a single response or multiple responses.

The expression of M F for multiple response models has been implemented in PFIM 3.0, an extension of the R function PFIM, dedicated to population design evaluation and optimization Regarding the comparison of the estimation methods for multiple response models, no problems of convergence were apparent with SAEM or the FO method for the 1000 simulated data files (see Section 3.3). However, it was difficult to fulfill convergence conditions with the FOCE method for which convergence was observed for only 853 of 1000 datasets. The simulation study illustrated the accuracy of the SAEM algorithm in the simultaneous approach, the parameter estimates being unbiased and with small RMSE. Similar results were observed for the FOCE method, but conclusions must be made with caution due to problems of convergence, as noted previously. Regarding FO, the large bias and RMSE already observed in the context of single response models [START_REF] Dartois | Evaluation of uncertainty parameters estimated by different population PK software and methods[END_REF] was also observed in this context of multiple response models with simultaneous estimation. The FO method produced larger RSE on all PD parameters compared to those computed with the FOCE method. This is in accordance with the recommendation to use the FOCE method instead of the FO method in this simultaneous estimation context [START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data I: best-case performance[END_REF][START_REF] Zhang | Simultaneous vs. sequential analysis for population PK/PD data II: robustness of methods[END_REF]. This apparent difference can be explained by considering the difference between FO and FOCE approximations. The FO method approximates the likelihood by linearizing the population model in its random effects about a value of zero whereas FOCE is defined by a linearization of the model around individual estimates of the random effects. The FOCE method uses more information than the FO method, which is an advantage for estimation in multiple response models. Although the estimation method FO performed badly, the same first order approximation in the computation of M F to predict SE performed well. The limitations of this first order approximation around the expectation of the random effects thus differ for design evaluation where derivatives of log-likelihood are computed and for parameter estimation.

Several studies have considered the single response, stressing the limitations of this linearization for design evaluation. Using a simple model with few random effects, Han [START_REF] Han | Optimal designs for Nonlinear regression Models with Application to HIV Dynamic Studies[END_REF] found that M F was quite different when computed by linearization than by an adaptative Gaussian quadrature method. Merlé et al. [START_REF] Merlé | Impact of pharmacokinetic-pharmacodynamic model linearization on the accuracy of population information matrix and optimal design[END_REF] compared the Fisher information matrix computed by linearization to one computed by stochastic simulation and showed that the linearization seems to have no impact on the population D-optimal designs obtained but only on the true efficiency of the designs. Whether via adapative Gaussian quadrature or via stochastic simulation, the evaluation of the Fisher information matrix without linearization is computationally intensive and is also limited to a matrix of small dimensions. Finally, in the context of Bayesian design where prior distributions are used, Han et al. [START_REF] Han | Bayesian experimental design for nonlinear mixed effects models with application to HIV dynamics[END_REF] proposed an attractive solution to compute M F ; however, it is also time consuming.

( ) ( ) 2 ; , i i F i T L Y M E ξ   ∂ Ψ Ψ = -     ∂Ψ ∂Ψ   (16)
where ( )

; i i L Y Ψ
is the log likelihood of the vector of observations i Y of that individual for the population parameters Ψ . Because F is nonlinear, there is no analytical expression for the log-likelihood ( )

; i i L Y Ψ
, and we use the first-order Taylor expansion of the model

( ) ( ) ( ) , , ,
i i i i F F g b θ ξ β ξ = , around the expectation of i b , that is to say around 0: ( ) ( ) ( ( , 0), ) , , ( ( , 0), )
T i i i i F g F g b F g b b β ξ β ξ β ξ   ∂ ≅ +   ∂   . (17)
With this approximation, the Equation ( 5) can be written as:

( ( , 0), ) ( ( , 0), )

T i i i i F g Y F g b b β ξ β ξ ε   ∂ ≅ + +   ∂   , (18)
Therefore, the approximated marginal expectation i E and variance i V of i Y are given by : ( ) ( ( , 0), )

i i i E Y E F g β ξ ≅ = (10) ( ) ( ) ( ) ( ) ( ) ( ) int , 0 , , 0 , , 0, , ,
T i i i i er slope i T F g F g Var Y V b b β ξ β ξ β σ σ ξ     ∂ ∂ ≅ = Ω + Σ         ∂ ∂     (11)
The log-likelihood i L is thus approximated by: ( )

( ) ( ) ( ) ( ) 1 2 ; ln 2 ln T i i i i i i i i i L Y n V Y E V Y E π - - Ψ ≅ + + - - (19)
Note that in this approximation, for sake of simplicity we have assumed that the variance of the error model is not linked to the random effects of the individual but only to the mean parameters. Based on this expression of the log-likelihood i L , we can derive the expression of an elementary Fisher information matrix for a multiple response model. For the sake of simplicity, we have omitted the indice i for the individual in the following. M F is a block matrix depending on the approximated marginal expectation E and variance V of the observations:

( ) ( ) ( ) ( ) ( ) , , 1 , , , 2
F T A E V C E V M C E V B E V ξ   Ψ ≅     (20)
where ( ( , )) ( )

ml l m V V C E V tr V V λ β - - ∂ ∂ = ∂ ∂ with ( ) 1

where int erk σ and slopek σ qualify the model=

  for the variance of the residual error of the th k response. The case 0 slopek σ returns a homoscedastic error model, whereas the case int 0 erk σ = 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 then note ( )

  4.1 and higher versions. The implementation of the population Fisher information matrix assumes that the variance of the observations with respect to the mean parameters is constant (see Appendix). PFIM 3.0 evaluates population designs in NLMEM with multiple responses and thus returns the expected standard errors, defined as the square roots of the diagonal elements of the inverse of M F , on the population parameters with the design evaluated. To use PFIM 3.0, some prior information has to be supplied by the user such as the structural model, its parameterization and a priori values of the parameters. PFIM 3.0 can also optimize population designs with different optimization options. More details are available in an extensive document that can be freely downloaded with the function PFIM 3.0 on the PFIM website [21].

  calculated the individual PK concentrations PK PK PK f (θ , t ) predicted by the model at each time PK t of PK ξ . We also computed the individual PK concentrations at each time PD t of PD then generated. Finally, for each response, we simulated the random errors PK ε and PD ε from a normal distribution with zero mean and variance derived from Equation (6) using the parameters slopePK σ and int erPD σ , respectively. Those errors were added to the previously generated PK and PD data to form the simulated observations for the PK and the PD response respectively.Using MONOLIX (Version 2.1) with SAEM as the estimation algorithm, we estimated the parameters using this simulated dataset and we then derived the observed population Fisher information matrix with the Louis's principle procedure and the linearization method of SAEM. For these two Fisher information matrices, we then transformed the observed SE for each component of the population vector Ψ obtained with a simulation of 10000 sim N = individuals into predicted SE of a population of 100 N = individuals to be adapted to the design of the example using ( th i component of Ψ.

  .08, 0.1, 1.5, 3, 0.8 ; the values of the variance of the random effects and the variance of the residual errors were the same as for the first set of initial parameters. The initial values of the parameters and the different elements required to use the SAEM algorithm were identical to those described in section 2.4.1.

parameter 2 50 Cω

 50 with a RSE close to 200%. In contrast, the empirical RSE for FOCE and SAEM were very close to the RSE predicted with SAEM_LI, SAEM_LO and PFIM.

[

  21-23]. Using a PKPD model, we have shown the appropriateness of the predicted RSE obtained with M F computed by PFIM 3.0 by comparison to those computed without any linearization by the SAEM algorithm implemented in the MONOLIX software (see Section 3.1). The predicted RSE were indeed all very close except for a slight discrepancy in the RSE for C 50 variability.The simulation study on the PKPD example showed the concordance between the RSE predicted by PFIM and the empirical RSE computed from estimation results with the FOCE or the SAEM algorithm using simulated datasets (see Section 3.2). Regarding results using the FO algorithm, the empirical RSE were much larger than those predicted by PFIM or obtained with FOCE or SAEM, in particular regarding the variability of the PD parameters. The distribution of the RSE obtained with FO, FOCE and SAEM from the simulated data files are in accordance with their respective empirical ones.

APPENDIX:

  Technical details in the development of M F in NLMEM for multiple responsesThe population Fisher information matrix models for the individual i with design i ξ is given by

  In this paper, we consider that the variance of the observations with respect to the mean parameters is constant, i.e. ( , ) 0 ml C E V = .

Table 1 .

 1 Population parameters values used for the simulated example Parameters

Figure 1 Figure 2 .

 12 Figure 1

Figure 3

 3 Figure 3.

  parameters. We associated a proportional error model with the PK model characterized by the

	parameter slopePK σ		and a homoscedastic error model with the PD model characterized by the
	parameter int erPD σ	. Thus, the vector of population parameters Ψ is described by the vector of
	the fixed effects	β	T	=	(	0 , , , C V E β β β β Cl	max E	,	β	50 C	)	and by T
	λ	T	0 2 , , E ω ω ω ω 2 2 , C Cl V	max E 2	,	50 2 C ω σ ,	slopePK	,	σ	int	erPD
	of		N =	100	individuals. They all had 3 sampling times at 0.166, 6 and 12 for PK and 4
	sampling times for PD at 0.166, 6, 12 and 20 hours. Therefore, we had one elementary design
	(	, ξ ξ PK	PD	)	with	PK ξ =	(	0.166, 6, 12	)	and	PD ξ =	(	0.166, 6, 12, 20	)	. The population design was
	thus defined by

λ the vector composed by the variance of the random effects and by the parameters for the error models such that ( ) = . The dose was fixed to 1 and the parameter values used in this paper are given in Table

1

.

We determined a population design associated with this PKPD example. This determination was empirical, without any optimization. The population design was composed of one group

Table 2

 2

	as relative SE, i.e. SE

50

C

β between the SE predicted by PFIM and the one given by SAEM_LO. Regarding the variance parameters, RSE were also very close, except for

Table 2 .

 2 Comparison of the relative standard errors (RSE in %) predicted by the SAEM algorithm implemented in MONOLIX V2.1 with both methods of computation of the SE (noted SAEM_LO for Louis's principle and SAEM_LI for the linearization method) and predicted by PFIM for the PKPD example.

	Parameters

Table 3 .

 3 Relative bias (%) and relative RMSE (%) of parameter estimates for datasets with fulfilling convergence conditions (denoted S) using the FO and the FOCE INTER methods implemented in NONMEM V software and the SAEM algorithm implemented in the MONOLIX V2.1

	software.
													Parameters
		Methods	S	β	Cl	C V β	β	0 E	β	max E	β	50 C	2 Cl ω	2 C V ω	ω	0 2 E	ω	max E 2	50 2 C ω	σ	slopePK	σ	int erPD
	Bias (%)	method FO	1000	-0.5	0.6	34.5	72.4	131.6	-10.9	-6.9	2.3	-1.7	29.8		-0.3	23.5
		FOCE	853	0.3	1.3	1.9	-3.5	-7.3	-0.9	-1.1	-0.5	-0.1	-1.4		1.9	-0.1
		SAEM	1000	0.0	-0.1	-0.1	-0.8	0.5	-0.4	-1.0	0.3	-0.1	-0.1		1.3	-0.1
	RMSE (%)	FO	1000	11.7	7.4	35.1	19.5	135.1	15.9	33.4 48.8	34.5	367.1	26.0	156.8
		FOCE	853	9.2	8.8	4.4	3.4	9.1	15.8	18.1 20.9	16.2	43.5	12.4	27.7
		SAEM	1000	5.3	5.6	5.3	3.3	4.0	17.8	16.4 21.7	16.7	41.5		7.8		6.0

In conclusion, the comparison of RSE predicted by PFIM to RSE computed without any linearization by SAEM as well as to RSE obtained from the simulation study, supports the appropriateness of the approximated Fisher information matrix for multiple response models.

Its implementation in PFIM 3.0 provides a useful computing tool for design evaluation and optimization in the development of PKPD or PK studies.

ACKNOWLEDGEMENTS

Part of this work was supported by a grant from F. Hoffmann La Roche Ltd, Basel, Switzerland. The authors would like to acknowledge Pr Marc Lavielle for his help regarding the use of the MONOLIX software as well as the IFR02 of INSERM and Hervé Le Nagard for the use of the "centre de biomodélisation".

Another alternative to the linearization method is to compute the expected Fisher information matrix using the SAEM estimation algorithm, as proposed here. A large dataset is simulated to be close to the asymptotic properties. The observed M F is then estimated for the estimation performed in this dataset using one of the two methods of deriving M F proposed for SAEM. This simulation study shows a broader distribution of the observed RSE when M F is computed by the Louis's principle [START_REF] Louis | Finding the observed information matrix when using the EM algorithm[END_REF] compared to the RSE observed with the linearization procedure. Nevertheless, correct results are obtained from both methods. Although the approach using the SAEM algorithm of MONOLIX version 2.1 can be applied for problems with a large number of random effects, it is time consuming (hours) compared to PFIM (seconds). It may be required in specific cases such as when design evaluation is used to predict the power of a test to detect covariates [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]. The appropriateness of the RSEs of PFIM combined with its fast execution emphasizes its advantage in cases of design optimization where a large number of M F often have to be computed. Moreover, in design evaluation and optimization for nonlinear models, some a priori values of the parameters are required. Often, they are not precisely known, and therefore the need to use exact methods to predict M F is questionable.

In this study, we empirically determined a population design associated with a simple PKPD example for which the dose was equal to 1. No change in the dose was envisaged because the main purpose of this work was to evaluate M F for a PKPD model associated with one population design. However, it would be interesting to study the influence of dose on the population design and thus to plan dose optimization in order to have a better understanding of the relationship between the PK and the PD model.

In the present study, we considered only the case of a diagonal Ω matrix with no correlation between the random effects of the PK and the PD parameters. This could be extended by exploring the appropriateness of M F when the PK parameters are directly correlated with the PD parameters, and thus the development of M F for a full Ω matrix. This development was performed by Mentré et al. [START_REF] Mentré | Optimal design in random effect regression models[END_REF] for the correlation of the random effects parameters of single response models and recently by Ogungbenro et al. [START_REF] Ogungbenro | Incorporating correlation in interindividual variability for the optimal design of multiresponse pharmacokinetic experiments[END_REF] for multiple response models.

Furthermore, the population Fisher information matrix implemented in PFIM 3.0 is approximated by a block diagonal matrix assuming that the variance of the observations with respect to the mean parameters is constant (see Appendix). It would therefore be interesting to investigate the influence of this assumption on the computation of the Fisher information matrix.

LIST OF FIGURES