Thi Huyen

Tram Nguyen

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

France Mentre

Thu Thuy Nguyen

Individual Bayesian information matrix for predicting estimation error and shrinkage of individual parameters accounting for data below the limit of quantification

Purpose: In mixed models, the relative standard errors (RSE) and shrinkage of individual parameters can be predicted from the individual Bayesian information matrix (MBF).

We proposed an approach accounting for data below the limit of quantification (LOQ) in MBF.

Methods: MBF is the sum of the expectation of the individual Fisher information (MIF)

which can be evaluated by First-Order linearization and the inverse of random effect variance.

We expressed the individual information as a weighted sum of predicted MIF for every possible design composing of measurements above and/or below LOQ. When evaluating MIF, we derived the likelihood expressed as the product of the likelihood of observed data and the

INTRODUCTION

Since their first development in early 1970, NonLinear Mixed Effect Models (NLMEM) have rapidly gained popularity and now are widely used to analyze pharmacokinetic/pharmacodynamic (PK/PD) longitudinal data [START_REF] Mould | Basic concepts in population modeling, simulation, and model-based drug development[END_REF][START_REF] Mould | Basic concepts in population modeling, simulation, and model-based drug development-Part 2: Introduction to pharmacokinetic modeling methods[END_REF][START_REF] Upton | Basic concepts in population modeling, simulation, and model-based drug development: Part 3-Introduction to pharmacodynamic modeling methods[END_REF]. NLMEM can compensate for the lack of individual information by borrowing the strength from the data in the whole population. Therefore, they allow for more precise parameter estimation even with limited individual design. PK/PD parameters obtained in population analysis are useful to characterize or predict the PK/PD response of a population, to simulate clinical trials [START_REF] Mould | Basic concepts in population modeling, simulation, and model-based drug development[END_REF][START_REF] Mould | Basic concepts in population modeling, simulation, and model-based drug development-Part 2: Introduction to pharmacokinetic modeling methods[END_REF][START_REF] Upton | Basic concepts in population modeling, simulation, and model-based drug development: Part 3-Introduction to pharmacodynamic modeling methods[END_REF].

These population parameters are usually obtained by Maximum Likelihood approach [START_REF] Dartois | Overview of modelbuilding strategies in population PK/PD analyses: 2002-2004 literature survey[END_REF]. Once population parameters are estimated, individual parameters can then be derived, using Bayesian approach which combines 'a priori' information of population parameters estimated previously and the individual data. Individual parameter estimates in NLMEM are often obtained as the maximum a posteriori (MAP), i.e., the mode of the posterior distribution. A good precision of individual parameter estimates is needed for screening covariate effect, for sequential PK/PD analysis, for individualized treatment, etc.

The precision of parameter estimates depends on the quantity and quality of the available data, which in turns, depends on the study design (i.e. the number and timing of the samples). To evaluate and optimize a design, two approaches have been proposed. The first approach, based on clinical trial simulation, is very time-consuming and therefore, is limited in term of designs that can be evaluated. Alternatively, one can use methods based on mathematical derivation of the Fisher information Matrix (MF) [START_REF] Mentré | Optimal design in random-effects regression models[END_REF]. This approach relies on the Cramer-Rao inequality, which states that the inverse of the MF is the lower bound of the variance-covariance matrix of any unbiased estimator of the parameters. Expressions of the expected individual (MIF) and population Fisher information matrix (MPF) using first-order (FO) approximation [START_REF] Mentré | Optimal design in random-effects regression models[END_REF][START_REF] Bazzoli | Fisher information matrix for nonlinear mixed effects multiple response models: evaluation of the appropriateness of the first order linearization using a pharmacokinetic/pharmacodynamic model[END_REF] have been developed and implemented in several software programs [START_REF] Nyberg | Methods and software tools for design evaluation in population pharmacokinetics-pharmacodynamics studies[END_REF][START_REF] Mentré | Current Use and Developments Needed for Optimal Design in Pharmacometrics: A Study Performed Among DDMoRe's European Federation of Pharmaceutical Industries and Associations Members[END_REF] to evaluate and optimize designs for standard individual regression or population analysis, respectively. Beside MIF and MPF, the expected Bayesian Fisher information matrix (MBF) was also developed to evaluate the estimation error of individual parameters obtained by MAP [START_REF] Merlé | Optimal sampling times for Bayesian estimation of the pharmacokinetic parameters of nortriptyline during therapeutic drug monitoring[END_REF][START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF]. In addition to the estimation error, MBF can also be used to predict the shrinkage [START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF][START_REF] Fedorov | Mixed models: design of experiments[END_REF] a metric quantifying the informativeness of the individual data and the reliability of individual parameter estimates [START_REF] Molenberghs | Linear Mixed Models for Longitudinal Data[END_REF][START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF][START_REF] Xu | Shrinkage in nonlinear mixedeffects population models: quantification, influencing factors, and impact[END_REF]. These developments on MBF were recently implemented in PFIM, an R program for design evaluation and optimization [START_REF] Bazzoli | Design evaluation and optimisation in multiple response nonlinear mixed effect models: PFIM 3.0[END_REF].

Data below the limit of quantification (LOQ) are frequently observed in PK/PD studies. While these data can be limited in well-designed PK/PD studies, they are present at large proportion in several situations, in particular when they are used as a marker of treatment effectiveness.

For instance, the efficacy of the treatment of hepatitis virus or human immunodeficiency virus is evaluated as the ability of the treatment to clear the virus or to maintain the viral load under the detection levels. In such case, data below the LOQ cannot be avoided and can have significant impact on the estimation error and shrinkage of individual parameters [START_REF] Nguyen | Influence of a priori information, designs, and undetectable data on individual parameters estimation and prediction of hepatitis C treatment outcome[END_REF] and therefore, by the methods used to account for these data.

Several approaches have been proposed to handle data below LOQ in parameter estimation [START_REF] Beal | Ways to fit a PK model with some data below the quantification limit[END_REF][START_REF] Bergstrand | Handling data below the limit of quantification in mixed effect models[END_REF][START_REF] Thiébaut | Estimation of dynamical model parameters taking into account undetectable marker values[END_REF][START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF] and model evaluation [START_REF] Nguyen | Extension of NPDE for evaluation of nonlinear mixed effect models in presence of data below the quantification limit with applications to HIV dynamic model[END_REF]. At the study design step, the impact of these censored data can be taken into account in clinical trial simulation approach but this method is timeconsuming. When using the MF-based approach, some methods have been proposed to handle data below LOQ in MPF [START_REF] Dumont | Optimal Sampling Times for a Drug and its Metabolite using SIMCYP® Simulations as Prior Information[END_REF][START_REF] Vong | Handling below limit of quantification data in optimal trial design[END_REF] but to our knowledge, no such approach has been proposed for MBF.

The main objective of this work was to propose and evaluate a method to account for data below VK models [START_REF] Guedj | Design evaluation and optimization for models of hepatitis C viral dynamics[END_REF][START_REF] Nguyen | A pharmacokinetic-viral kinetic model describes the effect of alisporivir as monotherapy or in combination with peg-IFN on hepatitis C virologic response[END_REF], was used as an illustrative example throughout the paper.

LOQ

METHODS

Model and notation

For a given individual i (i = 1,…, N), we denote yi the ni-vector of measures obtained with the elementary sampling design ) , , , (

i i i i f y      ) , ( (1)
where i  is the vector of p individual parameters in individual i. i  is the random error which follows a normal distribution with mean 0 and variance-covariance matrix

  i n j ij i i i t diag , , 1 2 ) , ( ) , ( 
      
denoted  for simplicity in the followings, with When using standard nonlinear regression, the vector of parameters to be estimated for each

) , ( ) , ( slope inter ij i ij i t f t       
individual i is T T T i ) , (     where ) , ( slope inter    
T is a vector of length r = 1 (in the case of additive or proportional error variance model) or of length r = 2 (combined error). When using nonlinear mixed effect modelling, the vector of population parameters to be estimated is

T T p T ) , ,..., , ( 2 2 1 
 . The MAP estimate of i  is       ) ( log ) | ( log argmax )) y | ( ( argmax ˆi i i i i i p y p p        (2) Then the MAP estimate of i  is ) , ( ˆi i g     .

Individual and Bayesian Fisher information matrix

The expected individual Fisher information in standard nonlinear models is defined by

             T y IF y p y p E M       ) | ( log ) | ( log ) , ( (3) where )
| (  y p
is the probability density of the observations y given the parameters

T T T i ) , (    
. From now on, for simplicity, subscript i is omitted for individual. By deriving the following expression

))) , ( ( )) , ( (( 2 1 - | | log 2 ) 2 log( 2 ) | ( log 1       f y f y tr n n y p T         (4)
each term of the Fisher information matrix can be obtained as

                          1 1 1 , 2 1 ) , ( ) , ( )) , ( (
v u u T v v u IF tr f f M           (5)
where u = 1,…, (p+r) and v = 1,…, (p+r).

For Bayesian estimation of the random effects, the expected Bayesian Fisher information taking into account the prior information on  is defined by

                                  T y T BF y p y p E E y p y p E M             ) | ( log ) | ( log ) | ( log ) | ( log ) ( | , (6) therefore
                                  T T y BF p p E y p y p E E M             ) ( log ) ( log ) | ( log ) | ( log ) ( | (7)
We recognized in the first term of expression [START_REF] Nyberg | Methods and software tools for design evaluation in population pharmacokinetics-pharmacodynamics studies[END_REF] the expectation of the predicted individual Fisher information matrix for parameter  over the distribution of the random effects. As  follows a Gaussian distribution, the second term is the inverse of the variance of the random effects. Therefore we can write

  1 ) ), , ( ( ) (          g M E M IF BF (8) The expectation   ) ), , ( (     g M E IF
can be evaluated by first order linearization of the model around the expectation of random effects (i.e., 0) in the expression of the likelihood. It has been

shown that ) ( BF M is then approximated as 1 ) , ( ) (     M M M M IF T BF    (9)
where M is the p p  identity matrix I if the individual parameters follow a normal distribution or the

p p  diagonal matrix ) ,..., ( 1 p diag  
if the individual parameters follow a log-normal distribution [START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF].

The standard error (SE) for individual and Bayesian estimation can be derived from the square root of the diagonal terms of MIF -1 and MBF -1 respectively. The shrinkage (Sh) for MAP estimation is quantified by the ratio of the estimation variance predicted by MBF -1 and the variance-covariance matrix of the random effects  , and can be calculated as the diagonal elements of the matrix [START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF][START_REF] Fedorov | Mixed models: design of experiments[END_REF]. These developments were implemented in the R tool PFIM 4.0 [START_REF] Mentré | PFIM 4.0: new features for optimal design in nonlinear mixed effects models using R[END_REF] and its interface version PFIM Interface 4.0 (www.pfim.biostat.fr).

1 1 ) ( ) (        BF M W I

Extension to account for data below LOQ

Data below LOQ can be taken into account in maximum likelihood estimation for nonlinear models as left-censored observations [START_REF] Beal | Ways to fit a PK model with some data below the quantification limit[END_REF][START_REF] Bergstrand | Handling data below the limit of quantification in mixed effect models[END_REF][START_REF] Thiébaut | Estimation of dynamical model parameters taking into account undetectable marker values[END_REF][START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF][START_REF] Nguyen | Extension of NPDE for evaluation of nonlinear mixed effect models in presence of data below the quantification limit with applications to HIV dynamic model[END_REF]. To account for these data in the expected individual Fisher information matrix as well, we first wrote the probability density of the data y given the parameters  as

| | ) | ( ) | ( ) | (    (11)
where

          ) , ( ) , ( ) | ( j j obs j obs j t t f y y p      and              ) , ( ) , ( ) | ( ) | (
      LOQ y cens j LOQ j obs j j y p y p y p ) | ( log ) | ( log ) | ( log j y |    (12)
The contribution of all censored observations to the Fisher information is given by

        LOQ j T cens j cens j LOQ j IF y p y p M cens j j y | y | ) | ( log ) | ( log ) , (       (13)
To calculate ) , (

j y | LOQ j IF cens M   
, we derived the log-probability of a measurement yj to be below LOQ with respect to each element p  of the individual parameter vector, using

               ) , ( ) , ( ) , ( ) , ( 1 ) |
( j j p j j p cens j t t f LOQ t f t y p           (14) Thus                ) , ( ) , ( ) , ( ) , ( 1 ) | ( log j j p j j p cens j t t f LOQ t f t y p           (15) where ) ( / ) ( ) ( z z z     .
The number of sampling times j at which yj are above LOQ, denoted k, can vary from 0 to n,

k n ) (   of (n-k) measurements below LOQ. For each value of k (k = 0,…,n), there are         k n possible designs with (n-k) censored values                     k n l cens k n l obs k l k l , , 1 ) ( ) ( ) (    
. The probability for each of these designs

) (k l  is given by       cens k n l obs k l j cens j j obs j k l y p y p P ) ( ) ( ) | ( ) | ( ) ( ) (      (16)
Consequently, the expected individual Fisher information matrix for design  taking into account possible data below LOQ can be written as

                  n k k n l cens k n l cens IF obs k l IF k l IF M M P M 0 1 ) ( ) ( ) ( , ) , ( ) ( ) , (        (17)
In the case of a decreasing function f, then for a given k, there is only one possible protocol

n k k cens k n l obs k l k l t t : ) 1 ( : 1 ) ( ) ( ) (         
for which the first k measurements from t1 to tk are above LOQ while all measurements from tk+1 to tn are censored. It can be shown that

                         ) , ( ) , ( ) , ( ) , ( ) ( 1 1 )
( k k k k k l k t t f LOQ t t f LOQ P p        (18)
and the expected Fisher information matrix can be written as

         n k n k cens IF k IF k IF t M t M p M 0 : ) 1 ( : 1 , ) , ( ) , (     (19) (note that            ) , ( ) , ( 1
1 0 t t f LOQ p    and             ) , ( ) , ( 1 n n n t t f LOQ p    ).
By using this extension of MIF in the formulae of MBF in equation [START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF], we obtain the expression of the expected Bayesian Fisher information matrix taking into account data below LOQ. These developments were implemented in R.

SIMULATION STUDY

PK-VK models and studied scenarios

A simulation study was conducted to evaluate the relevance of MIF and MBF in terms of predicted RSE and shrinkage of individual parameters. For this purpose, we considered a PK-VK model inspired from [START_REF] Guedj | Design evaluation and optimization for models of hepatitis C viral dynamics[END_REF][START_REF] Nguyen | A pharmacokinetic-viral kinetic model describes the effect of alisporivir as monotherapy or in combination with peg-IFN on hepatitis C virologic response[END_REF], describing the concentration of peg-IFN and the viral kinetics during a four-week treatment. Peg-IFN was given subcutaneously once weekly at the dose of 180 µg/L. The model was written in ordinary differential equation (ODE) system as follows:

() () a a e dA kA dt dQ k A k Q dt Qt Ct V    50 0 () () () (1 ( ))
I I II Ct Et C t EC dI V T I dt dV E t p I cV dt        ( 20
)
where A is the quantity of peg-IFN at the site of absorption, Q and C are respectively the For this PK-VK model, fixed effect parameters were chosen basing on a previous PK-VK model [START_REF] Nguyen | A pharmacokinetic-viral kinetic model describes the effect of alisporivir as monotherapy or in combination with peg-IFN on hepatitis C virologic response[END_REF]. As only short-term viral load data were considered, pI and  cannot be identified and therefore, were fixed at 10 day -1 and 5.7 10 -6 day -1 .virion -1 , respectively. The vector µ of fixed effect parameters includes 7 parameters µ=(ka, ke, V, EC50, VL0, , c)= (0.9 day -1 , 0.15 day -1 , 12 L, 0.15 µg/L, 10 6 copies/mL, 0.2 day -1 , 6 day -1 ).

For the random effect model, we considered an exponential model, i.e., a log-normal distribution of individual parameters, with the same standard deviation  for all parameters (ka, ke, V, EC50, VL0, , c) and no correlation between random effects. An additive error model on the concentrations for PK response and on the log10-scale of the VK response  = (PK ,VK) was used in the simulation. Our standard scenario was defined with moderate inter-individual variability (=0.3) and residual error (PK = 0.5 µg/L, VK = 0.15 log10 copies/mL).

Three designs were considered: a rich design 5-7, with PK = (1, 3, 5, 7, 14 days) and VK = (0, 1, 3, 5, 7, 14, 28 days), a sparse design 3-5 with PK = (1, 3, 7 days) and VK = (0, 1, 3, 7, 28 days) and a very sparse design 2-3 with PK = (3, 7 days) and VK = (0, 7, 28 days). We set the limit of quantification (LOQ) at 100 copies/mL on the viral kinetic response.

Evaluation of MIF and MBF

To evaluate the prediction of MIF without or with data below LOQ, we predicted the RSE of individual parameter estimates using MIF then we compared these predictions with the RSE obtained using simulation. For the simulation study, we evaluated the typical individual profile using the fixed effect parameters and the model. We then simulated 1000 vectors of residual errors and added them to the typical profile to obtain 1000 datasets. Individual parameters were estimated for each dataset by standard nonlinear regression using R. Data below LOQ were kept at the simulated values in the "absence of data below LOQ" scenarios and were censored at LOQ values in the "presence of data below LOQ" scenarios. The likelihood for the studied multiple response model (PK and VK) written using differential equations was implemented in an R function (Supplementary material 1). This function allowed us to evaluate the contributions of observed and censored data respectively in the likelihood as given in Eq. 11 &

12. The package marqLevAlg was used for likelihood maximization [START_REF] Commenges | marqLevAlg: An algorithm for least-squares curve fitting[END_REF]. To calculate the RSE for each parameter in the simulation study, we computed the standard deviation of the 1000 individual parameter estimates and divided them by the value of the corresponding fixed effect parameter that was used to simulate data.

Similarly, for the evaluation of MBF without or with of data below LOQ, we predicted RSE and shrinkage using MBF and then compared these predictions with those obtained by simulation. Data simulation as well as calculation of MIF and MBF were performed using R 3.2.

Influence of designs, inter-individual and residual variability on individual parameter estimation

To evaluate the influence of designs and of different levels of variability on the shrinkage and RSE of individual parameters, we considered additional scenarios with higher inter-individual variability (=0.5) and/or residual errors (=(1,0.3)). In this part of the study, we presented the results for the scenarios without data below LOQ as the same trends would be observed in presence of data below LOQ.

RESULTS

Data simulation

The typical PK and VK profiles and the sampling times of different designs were shown in comparable with values obtained in literature (between 26%-68% depending on the LOQ level [START_REF] Hartwell | Peginterferon alfa and ribavirin for chronic hepatitis C in patients eligible for shortened treatment, re-treatment or in HCV/HIV co-infection: a systematic review and economic evaluation[END_REF]).

Evaluation of MIF via simulation

The predicted RSE by MIF and those obtained from simulation are reported for all parameters in Figure 2. Of note, convergence was obtained for 959/1000 and 912/1000 datasets in absence or presence of data below LOQ, respectively, and 903/1000 datasets in both cases. Overall, MIF correctly predicted RSE which were close to those obtained with simulation, except for the parameter c, where the RSE provided by MIF and simulation were 42% and 55%, respectively in absence of data below LOQ and were 44 and 66%, respectively in presence of data below LOQ. MIF tended to under-predict the RSE, especially for parameters associated with large RSE. MIF predicted an increase of RSE of the parameter  from 9% to 16% due to the presence of data below LOQ, which was in agreement with the results observed in simulation: 9% vs 14% in absence and in presence of data below LOQ, respectively.

Evaluation of MBF via simulation

RSE predicted from MBF using population parameters or obtained using clinical trial simulation are reported for all parameters in Figure 3. Based on expression of MBF (Eq. 9), the highest RSE predicted by MBF for each parameter cannot exceed the inter-individual variability, which was equal to 30%, even for very sparse design 2-3 or for the most difficultly identifiable parameters such as c. These predictions of MBF were very close to those obtained with clinical trial simulation for most scenarios and parameters. MBF also tended to under-predict the RSE for scenarios and parameters associated with largest RSE (close to 30%). For instance, the predicted RSE was lower than the first quartile of RSE obtained via simulation for the parameter ka in the design 2-3 and for the parameter c in all designs. Similarly, the shrinkage was well predicted from MBF for most parameters and scenarios (Figure 4). As expected, the RSE and shrinkage increased with the decrease of the quantity and quality of individual information, i.e., when the number of samples per subject was reduced or when data below LOQ appeared in the samples. Taking into account these censored data by the new approach adequately predicted the increase of RSE and shrinkage due to the presence of data below LOQ for , the parameter the most sensitive to data below LOQ.

Influence of designs, inter-individual and residual variability

Figure 5-A and 5-B show the RSE predicted from the MIF and MBF, respectively for all parameters in different scenarios without data below LOQ. Of note, for the very sparse design, MIF was not invertible, meaning that individual parameters were numerically non-identifiable with this design. Overall, with sparser design or higher residual errors, the RSE predicted either by MBF or MIF increased. As expected, the influence of losing informativeness, either by reducing the number of samples or increasing the error levels (), on the RSE was much lower for MBF. With higher inter-individual variability (ω), the RSE obtained with MBF increased but those predicted by MIF remained unchanged as expected. Figure 5-C shows the shrinkage of all the parameters predicted from MBF for different scenarios. As expected, the shrinkage increased with the loss of information, i.e., as the design became sparser or residual errors increased. On the contrary, as shrinkage quantifies ratio between individual information and the prior information (ω), it decreased when ω increased.

DISCUSSION

With several approaches accounting for data below LOQ in parameter estimation of NLMEM [START_REF] Beal | Ways to fit a PK model with some data below the quantification limit[END_REF][START_REF] Bergstrand | Handling data below the limit of quantification in mixed effect models[END_REF][START_REF] Thiébaut | Estimation of dynamical model parameters taking into account undetectable marker values[END_REF][START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF], these data can be easily handled in clinical trial simulation approach. However, as this approach is time-consuming, we proposed in this article a method to take into account data below LOQ in the expected individual Fisher information matrix (MIF) and individual Bayesian information matrix (MBF) evaluated using FO approximation. Extension of MIF and MBF taking into account possible data below LOQ was proposed by deriving a full likelihood distinguishing the contribution of the "observed" data and the probability of the "censored" data (data below LOQ). Such approaches are important in PK-PD studies, particularly in infectious diseases where data below LOQ are common and often used as the marker of treatment efficacy and outcome. A PK-VK model was used to evaluate the relevance and illustrate the use of these developments.

The illustrative example was a complex PK-VK model defined by an ODE system with two responses, one describing the PK of peg-IFN and the other describing VK of hepatitis C virus in response to treatment. To our knowledge, this was the first time that MIF and MBF were used to predict the RSE and shrinkage for individual parameters of a complex ODE model, without and with data below LOQ. In this model, the presence of data below LOQ influences in particular the estimation of the parameter  because this parameter defines the second phase decline of the viral kinetics. As it represents the long-term decline rate of viral load, a good estimation of this parameter was shown to be important to predict treatment outcome or treatment duration needed to achieve viral eradication, at least under peg-IFN based treatment [START_REF] Nguyen | Influence of a priori information, designs, and undetectable data on individual parameters estimation and prediction of hepatitis C treatment outcome[END_REF][START_REF] Guedj | Second-phase hepatitis C virus RNA decline during telaprevir-based therapy increases with drug effectiveness: Implications for treatment duration[END_REF][START_REF] Nguyen | HCV kinetic models and their implications in drug development[END_REF]. Of note, in our example, data below LOQ only appeared at the end of the kinetics, which is a common case in viral kinetic data. However, there are also situations that data below LOQ appear in the early phase or in the middle of the studied kinetics. Such situations can be handled with this method, as described in equations ( 16) and ( 17).

In the simulation study, we showed that the proposed extension of MIF and MBF correctly predicted the increase of RSE and/or shrinkage, especially of the parameter δ, in presence of data below LOQ. Although the RSE of this parameter remains quite satisfactory in this simulation, the increase of RSE will be more dramatically influenced by larger proportion of data below LOQ. Because of correlations in the estimates, the estimation of c and EC50 would also be influenced by data below LOQ.

In this example, MIF and MBF tended to slightly under-predict RSE, especially for parameters or designs that were associated with high estimation errors. For instance, MIF predicted the RSE for the parameter c with 15 to 20 points lower than the RSE obtained by simulation, in absence or presence of data below LOQ, respectively. However, even in the case of under-prediction, the RSE were still predicted at high values (for instance, higher than 40% for the parameter c) which were able to indicate the imprecision of estimation for these parameters. MBF also under-predicted the RSE of the parameters ka and EC50 in the very sparse design and of the parameter c in all three designs. Although the RSE and shrinkages were slightly underpredicted for some parameters or scenarios, in general, MIF and MBF computed by FO approximation performed quite well as the predicted RSE and shrinkage were close to those obtained by simulation for several parameters. However, in this simulation, we only considered a moderate variance for random effects (30% to integrate the derivatives of log-likelihood over the random effects and Monte Carlo simulation to evaluate its expectation [START_REF] Riviere | An MCMC method for the evaluation of the Fisher information matrix for non-linear mixed effect models[END_REF]. It would also be useful to combine both population FIM (MPF) and Bayesian FIM (MBF) using a compound optimality criterion [START_REF] Atkinson | DT-optimum designs for model discrimination and parameter estimation[END_REF][START_REF] Nguyen | Design optimisation for pharmacokinetic modeling of a cocktail of phenotyping drugs[END_REF] to set a balance between the two matrices and optimize individual and population parameters at the same time.

Finally, we illustrated the influence of different types of variability on the shrinkage and RSE of individual parameters, for the first time in an ODE multiple response model. In accordance with previous results for single response PK models [START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF][START_REF] Xu | Shrinkage in nonlinear mixedeffects population models: quantification, influencing factors, and impact[END_REF], we found that shrinkage had a direct relationship with the residual variability and an inverse relationship with the inter-individual variability: an increase of residual errors led to higher shrinkage but a larger inter-individual variability would lead to a reduced shrinkage. On the contrary, RSE always had a direct relationship with the two types of variability and hence, would increase with residual errors and inter-individual variability. This result once again illustrated that shrinkage was a relative measurement of the ratio between individual information and prior information contained in the population model. Therefore, a shrinkage value of a parameter should be interpreted in regard of its inter-individual variability and a common level defining a high or low shrinkage may not exist. RSE may be a better criterion evaluate the quality of individual parameter estimation.

In conclusion, MBF obtained by FO linearization is a useful and rapid method to predict standard errors and shrinkages of individual parameters, in absence or presence of data below the quantification limit. It allows for evaluation of a large number of scenarios without extensive simulations. These developments will be implemented in the next version of PFIM. and =(0.5,0.15). The designs 5-7, 3-5, 2-3 are presented in green, blue and red, respectively.

Legend to Figures

For the same design, dark and light color represent the results obtained in absence and in presence of data below LOQ, respectively. 26

  in MIF and MBF. Relative standard error (RSE) and/or shrinkage of individual parameters were predicted from MIF or MBF in absence or presence of data below LOQ and compared with those obtained in a simulation study. In the second part of this work, we studied the influence of design, different levels of inter-individual and residual variability on the precision of individual parameter estimates using the prediction of MBF and MIF. A pharmacokinetic/viral kinetic (PK-VK) model describing the concentration of pegylated-Interferon (peg-IFN) and the decline of viral load of hepatitis C virus (HCV) under treatment, inspired from previous PK-

  the observations in each individual are assumed to be independent. When j y is above LOQ, its value obs j y can be observed. When j y is below LOQ, we can only observe the censored value LOQ and we denote by cens

  and  being the probability density function and the cumulative density function respectively of the standard normal distribution. Thus we can write the log-likelihood of the individual data as

M

  accounting for data below LOQ is a weighted sum of the expected information for a design composed of a set obs k ) (  of k measurements above LOQ and a set

  cens

  quantity and concentration of peg-IFN in the central compartment. Peg-IFN is absorbed into central compartment with a first order absorption constant ka and eliminated from this compartment with an elimination constant ke. V is the distribution volume of the central compartment. The VK model considers two populations of hepatocytes, the target cells, infected cells, I. Infected cells are cleared with a rate . The free virions VI are released from the infected cells at a rate pI per cell per day and are cleared from the circulation with a rate c. In this model, peg-IFN blocks the production of new virus with an effectiveness E(t). The relationship between peg-IFN concentration and effectiveness is described by an Emax model. The initial conditions for this system are defined as (A(0) = Dose, C(0) = 0, I(0) = cVL0/pI, VI(0)= VL0), where VL0 is the baseline viral load.

For

  this evaluation, 1000 profiles of 1000 individuals were simulated as follows: first, we simulated 1000 individual parameters using the fixed effect parameters and the inter-individual variability. Second, we added simulated residual errors to the individual profiles obtained from the simulated individual parameters and the model. Data below LOQ were kept at the simulated values in the "absence of data below LOQ" scenarios and were censored at LOQ values in the "presence of data below LOQ" scenarios. Each individual parameters p i, simulated individual i (i = 1,…, 1000) were estimated as MAP[START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] by fixing the population parameters at the values used for simulation. MAP estimation was performed using MONOLIX 4.3.3 (http://www.lixoft.eu/). The observed RSE of each individual parameter estimate was calculated by dividing the standard deviation of the posterior distribution The observed shrinkage for parameter p was calculated from the empirical variance of 1000 estimated random effects p i,

Figure 1 .

 1 Figure 1. The percentages of data below LOQ for the VK response observed at day 28 were

  Figure 1. Typical PK-VK profiles with the illustrative model and sampling times corresponding

Figure 2 .

 2 Figure 2. Relative standard errors (RSE) predicted from MIF (diamond points) and obtained

Figure 3 .

 3 Figure 3. Relative standard errors (RSE) predicted from MBF (diamond points) and obtained

Figure 4 .

 4 Figure 4. Shrinkages predicted from MBF (diamond points) and obtained from clinical trial

Figure 5 .

 5 Figure 5. Influence of designs, inter-individual variability and residual errors on the relative

Figure 1 . 24 Figure 2 .

 1242 Figure 1. Typical PK-VK profiles with the illustrative model and sampling times corresponding to three evaluated designs5-6, 3-5, 2-3

Figure 3 .

 3 Figure 3. Relative standard errors (RSE) predicted from MBF (diamond points) and obtained from clinical trial simulation (boxplot) for all the parameters for different designs and the scenario in which =0.3 and =(0.5,0.15). The designs 5-7, 3-5, 2-3 are presented in green, blue and red, respectively. For the same design, dark and light color represent the results obtained in absence and in presence of data below LOQ, respectively.

Figure 4 .

 4 Figure 4. Shrinkages predicted from MBF (diamond points) and obtained from clinical trial simulation (barplot) for all the parameters for different designs and the scenario in which =0.3 and =(0.5,0.15). The designs 5-7, 3-5, 2-3 are presented in green, blue and red, respectively. For the same design, dark and light color represent the results obtained in absence and in presence of data below LOQ, respectively.

27

 27

Figure 5 .

 5 Figure 5. .Influence of designs, inter-individual variability and residual errors on the relative standard errors (RSE) and shrinkage of individual parameters for different designs: 5-7, 3-5, 2-3 (presented in green, blue and red, respectively) and different scenarios: =0.3,=(0.5,0.15) (no fill pattern bar); =0.3,=(1,0.3) (vertical hatched bar); =0.5,=(0.5,0.15) (left hatched bar); =0.5,=(1,0.3) (right hatched bar). (A) RSE predicted from MIF, (B) RSE predicted from MBF, (C) Shrinkage predicted from MBF.

  ). It was already shown that increase of inter-individual variability could deteriorate the predictions of RSE and shrinkage predicted from MBF obtained with FO approximation[START_REF] Combes | Prediction of shrinkage of individual parameters using the bayesian information matrix in non-linear mixed effect models with evaluation in pharmacokinetics[END_REF]. In such case as well as for more complex nonlinear models or studies with discrete data, alternatives to FO approximation are needed. A perspective of this work is to use Markov Chains Monte Carlo approach to evaluate the

	expectation		( ( , ), )	

IF E M g

    over the distribution of the random effects, which may help to obtain better evaluation of MBF and therefore more accurate predictions of RSE and shrinkage. Moreover, calculation of the contribution of data below LOQ to MIF and MBF could be extended to the population Fisher information matrix by using Markov Chains Monte Carlo