Sarah Baklouti

Peggy Gandia

Didier Concordet
email: didier.concordet@envt.fr

"De-Shrinking" EBEs: The Solution for Bayesian Therapeutic Drug Monitoring

Background Therapeutic drug monitoring (TDM) aims at individualising a dosage regimen and is increasingly being performed by estimating individual pharmacokinetic parameters via empirical Bayes estimates (EBEs). However, EBEs suffer from shrinkage that makes them biased. This bias is a weakness for TDM and probably a barrier to the acceptance of drug dosage adjustments by prescribers. Objective The aim of this article is to propose a methodology that allows a correction of EBE shrinkage and an improvement in their precision. Methods As EBEs are defined, they can be seen as a special case of ridge estimators depending on a parameter usually denoted λ. After a bias correction depending on λ, we chose λ so that the individual pharmacokinetic estimations have minimal imprecision. Our estimate is by construction always better than EBE with respect to bias (i.e. shrinkage) and precision.

Results

We illustrate the performance of this approach with two different drugs: iohexol and isavuconazole. Depending on the patient's actual pharmacokinetic parameter values, the improvement given by our approach ranged from 0 to 100%. Conclusion This innovative methodology is promising since, to the best of our knowledge, no other individual shrinkage correction has been proposed.

Key Points

This study provides a methodology to suppress the individual shrinkage of the usual individual estimates of the pharmacokinetic parameters (empirical Bayes estimate). This method decreases the estimate's imprecision.

The use of this methodology is illustrated for the estimations of the individual pharmacokinetic parameters of isavuconazole and iohexol.

Introduction

Therapeutic drug monitoring (TDM) is used for some drugs to control drug exposure by choosing an appropriate dosage regimen (i.e. drug dosage and frequency of administration) [START_REF] Jang | Therapeutic drug monitoring: a patient management tool for precision medicine[END_REF].

Classical TDM checks whether or not a concentration is within a concentration range if measured at a predefined time and considers the exposure (i.e. the area under the curve) to be abnormal when the concentration is outside this range. It is very simple to use in practice. However, this classical TDM assumes a perfect correlation between the measured concentration and the associated exposure, which is not the case for all patients. Indeed, two patients could have two very similar concentrations for two samples collected at the same time and still have different pharmacokinetic profiles and hence different exposures. Interpreting a measured concentration without considering the physio-pathological characteristics of the patients on the one hand and by referring to the same range of concentrations for all patients on the other hand can lead to inappropriate dose adjustments. The use of a population pharmacokinetic (PPK) model to analyse the concentrations allows a better evaluation of the patient's exposure than classic TDM, through an estimation of their own pharmacokinetic parameters [START_REF] Donagher | Individualised medicine: why we need Bayesian dosing[END_REF][START_REF] Diaz | The mathematics of drug dose individualization should be built with random-effects linear models[END_REF][START_REF] Sheiner | Bayesian individualization of pharmacokinetics: simple implementation and comparison with non-Bayesian methods[END_REF]. These estimates are generally the empirical Bayes estimates (EBEs), computed by assuming that all the PPK model parameters are known without any imprecision [START_REF] Schumacher | Bayesian approaches in pharmacokinetic decision making[END_REF][START_REF] Merlé | Bayesian design criteria: computation, comparison, and application to a pharmacokinetic and a pharmacodynamic model[END_REF].

Unfortunately, EBEs can suffer from shrinkage [START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF]. Shrinkage is an average deviation of the individual pharmacokinetic (IPK) parameter estimate toward the population mean. Because shrinkage is a kind of bias (i.e. a systematic deviation), it should be possible to remove it. For example, Lavielle and Riba [START_REF] Lavielle | Enhanced method for diagnosing pharmacometric models: random sampling from conditional distributions[END_REF] proposed a method to limit the shrinkage when validating a PPK model. In the case of TDM, the problem is not an average but an individual deviation, socalled individual shrinkage. This individual shrinkage is specific to each patient (i.e. to their own pharmacokinetic parameter values), to the time when the blood samples are collected, and to the population's parameter values; this is why it is not usually removed [START_REF] Xu | Shrinkage in nonlinear mixed-effects population models: quantification, influencing factors, and impact[END_REF].

Instead, an average shrinkage is provided for each pharmacokinetic parameter after a PPK model analysis. Usually this shrinkage is expressed as a percentage. When this average shrinkage is large, the general recommendation is to discard or at least not trust the IPK parameters calculated for a patient from only a few concentrations [START_REF] Xu | Shrinkage in nonlinear mixed-effects population models: quantification, influencing factors, and impact[END_REF][START_REF] Bertrand | Pharmacogenetics and population pharmacokinetics: impact of the design on three tests using the SAEM algorithm[END_REF]. Thus, a large shrinkage can be a limitation of TDM based on the estimation of the IPK parameters, but it is not the only one. Indeed, it has been shown [START_REF] Mentré | Optimal design in random-effects regression models[END_REF][START_REF] Duffull | Evaluation of the pre-posterior distribution of optimized sampling times for the design of pharmacokinetic studies[END_REF][START_REF] Nyberg | Serial correlation in optimal design for nonlinear mixed effects models[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] that it is not reasonable to expect a precise estimation of the pharmacokinetic parameters for a patient with only a few concentrations obtained from the patient, even using an optimal sampling strategy [START_REF] Nguyen | Individual Bayesian information matrix for predicting estimation error and shrinkage of individual parameters accounting for data below the limit of quantification[END_REF][START_REF] Dumont | PFIM 4.0, an extended R program for design evaluation and optimization in nonlinear mixed-effect models[END_REF][START_REF] Jonge | Individualised cancer chemotherapy: strategies and performance of prospective studies on therapeutic drug monitoring with dose adaptation: a review[END_REF]. Even without any shrinkage, several combinations of pharmacokinetic parameters can lead to the drug concentrations observed in a patient. One way to quantify this imprecision is thus the length of the interval, or more generally, the volume of the prediction region of all the possible values of pharmacokinetic parameters that are consistent with the observed concentrations. This imprecision is intrinsic to the way the IPK parameters are estimated; in other words, this imprecision is intrinsic to the EBEs. Using the EBEs to estimate the IPK parameters cannot therefore lead to any significant improvements.

The purpose of this work is to propose an estimate of IPK parameters with no shrinkage and with less imprecision than the EBEs.

Materials and Methods

Models and Notations

When a PPK analysis is performed, a model describing the evolution of concentrations over time is used. This model is often written as in Eq. [START_REF] Jang | Therapeutic drug monitoring: a patient management tool for precision medicine[END_REF].

where Y ij is the jth concentration measured in individual i at the time t ij , i is a vector containing the individual pharma- cokinetic parameters for the ith individual, f and g are two known functions, ij is the jth residual error for individual i at the time t ij assumed to be independent of i , h is a known function usually assumed to be h i , , A i = A i e i , A i is a known matrix of covariates of the ith individual, is a fixed effect vector, i is the vector of the random effects involved in inter-individual variability for the ith individual, Ω is the variance-covariance matrix of the i s, and b is the vector of the parameters involved in the residual error model.

Data analysis with a PPK model as described by Eq. ( 1) provides estimates of , Ω , and b , respectively denoted by θ, Ω , and b . Of course, these estimates are imprecise. How- ever, we will not take into account this uncertainty. In other words, herein, we consider these estimates to be fixed.

Once the model has been created and validated, it can be used daily for hospitalised patients. So, let us consider the concentrations (Y j ) j=1,…,n observed at the time (t j ) j=1,..,n for a hospitalised patient with covariates A for whom TDM is requested. TDM should be based on the patient's IPK parameters 𝜑 = h 𝜂, θ, A , which are unknown when the blood sample is performed. However, a prediction of can be computed using both population parameters of the model and the patient's concentrations Y j j=1,…,n .

A classical predictor of is the EBE obtained as φ = h η, θ, A , where (Eq. 2)

In a more detailed manner (Eq. 3),

The first term (T1) of the previous equation is a least squares term. The second term (T2) corresponds to a penalty. It penalises values deviating far from the population average.

In Eq. ( 3), when n increases, T2 decreases, whereas T1 tends towards a fixed value. In this case, the patient's concentrations govern the prediction and the shrinkage should be low. By contrast, when n is low, T2 dominates T1 and the population parameters dominate the prediction, i.e. the is closer to the population mean 0 than it should be. In this case, one can expect η to be "shrunk" toward 0 . Equation (3) can be interpreted as a special case of a criterion giving a ridge estimator, defined as (Eq. 4):

(1)

Y ij = f t ij ; i + g t ij ; i ;b ij where ij ∼ iid ℕ(0, 1) i = h i , , A i where i ∼ iid ℕ(0, Ω) (2) η = argsup 𝜂 P 𝜂|Y 1 , … , Y n (3) η = arg inf 𝜂 1 n ∑ n j=1 Y j -f t j , h 𝜂, θ, A 2 g 2 t j , h 𝜂, θ, A , b + ln g 2 t j , h 𝜂, θ, A , b ���������������������������������������������� ���������������������������������������������� � (T1) + 1 n 𝜂 � Ω-1 𝜂 ��� ��� � (T2)
Actually, EBE(𝜑) = φ(1) = h η(1), θ, A . The only difference between Eqs. (3) and ( 4) is the presence of in T2 . If the predictions φ(𝜆) of were good, φ(𝜆) should be close to , which is unknown. To circumvent this difficulty (i.e. being unknown), we evaluate the distance between φ(𝜆) and using simulations by proceeding as follows:

1. A set of K vectors * 1 , * 2 , … , * K are simulated where 𝜂 * k ∼ iid ℕ 0, Ω for k ∈ 1, … , K. Then, for each k ∈ 1, … , K , 𝜑 * k = h 𝜂 * k , θ, A are computed. 2. For each k ∈ 1, … , K , Y * kj = f t j ;𝜑 * k + g t j ;𝜑 * k ; b 𝜀 * j are simulated where * j ∼ iid ℕ(0, 1), for j = 1,…,n. 3.
For a given , one computes η * k (𝜆) from Eq. ( 4) by using Y * kj instead of Y j .

In other words, for a vector A containing covariate values for a hospitalised patient, it is possible to simulate K different sets of IPK parameters * k (line 1). For each set of IPK parameters, it is possible to simulate concentrations Y * kj at the same sampling times as those of the hospitalised patient (line 2). From these concentrations Y * kj , it is possible to estimate IPK parameters η * k (𝜆) thanks to a Newton-Raphson algorithm (Eq. 4) for a given (line 3). Therefore, for K simulated patients, we have actual IPK parameters

( * k ) and predicted IPK parameters η * k (𝜆) (line 3). Once this η * k (𝜆) is obtained, one can compute φ * k (𝜆) = h η * k (𝜆), θ, A . We can see whether φ * k (𝜆) is a good predictor of * k by drawing * k as a function of φ * k (𝜆) or equivalently * k as a function of η * k (𝜆).
If it is a good predictor, the points should be concentrated around the first bisector line. Unfortunately, there may exist a systematic departure from this line, signalling a bias. This is why we look for a function Γ ∈ L , which minimises the distance between the vectors * and Γ( η * (𝜆)) (Eq. 5):

Depending on the space L where the Γ functions are to be found, there should exist at least one function Γ λ so that the minimum is reached. There exist many possible choices for L . As any continuous functions can be locally approximated by polynomials, we chose L as the set of the local polyno- mial functions. This polynomial function can be estimated, for example, by locally estimated scatterplot smoothing ( 4)

η(𝜆) = arginf 1 n n ∑ j=1 Y j -f t j , h 𝜂, θ, A 2 g 2 t j , h 𝜂, θ, A , b + ln g 2 t j , h 𝜂, θ, A , b + 𝜆 1 n 𝜂 � Ω-1 𝜂. (5) Γ 𝜆 = arginf Γ∈L K ∑ k=1 ‖ ‖ ‖ 𝜂 * k -Γ η * k (𝜆) ‖ ‖ ‖ 2 .
(LOESS) regression. Next, we look for the value λ , so that (Eq. 6):

Because ∈ R + , Eq. ( 6) can be solved by looking for lambda in an adaptive discretisation of ]0; max ]. For the practical implementation, max can be chosen close to 10.

At the end of these two optimisation steps, the best prediction with respect to the mean square error of their actual η for all patients having the same characteristics (i.e. same covariates values) and the same sampling times, is given by Γ λ η λ . Obviously, this predictor is not exact: several values from the same population correspond to a given value of η λ obtained from the hospitalised patient (blue double arrow in Fig. 1). The length of this range can be used as a prediction interval of P 𝜂 * |Γ λ( η * λ and so the uncertainty of the * estimation. A more traditional way of evaluating this uncertainty is to compute empirically a 95% prediction interval defined as the 2.5 and 97.5% percentiles of the * s corresponding to the η λ obtained from the patient of interest. The way to build these intervals (regression quantile methods) is not discussed in this article.

These ideas can be illustrated by the graph as shown in Fig. 1. To simplify, we concentrate on the component of * k linked to the clearance of the drug. We denote * Cl,k the kth simulated η for the clearance and we restrict Γ λ as a function of clearance only.

The graph in Fig. 1 can be disturbing because intuitively, we would tend to represent the predicted values as a function of the simulated/true values. However, the aim is to minimise the difference between the true value of Cl and its prediction (i.e. between * Cl and Γ λ η * Cl (𝜆) ). This presents itself as a classic nonlinear regression problem, in which we minimise the vertical distance between a point and its predictor.

A crude but global evaluation of the improvement that can be obtained using our prediction instead of the usual EBE can be computed by Eq. ( 7): This non-negative quantity evaluates how much closer our prediction Γ λ η * Cl,k λ is to the actual * Cl,k than is the EBE ( =1) given by η * Cl,k (1) . When this improvement is close to 0, the EBE and our predictor are close. As a consequence, our predictor does not improve the estimation of the individual pharmacokinetic parameter. On the contrary, any positive improvement shows that our predictor is better than

(6) λ = arginf 𝜆 K ∑ k=1 ‖ ‖ ‖ 𝜂 * k -Γ 𝜆 η * k (𝜆) ‖ ‖ ‖ 2 . (7) Improv = 1 - ∑ K k=1 � 𝜂 * Cl,k -Γ λ� η * Cl,k � λ� �� 2 ∑ K k=1 � 𝜂 * Cl,k -η * Cl,k (1)
� 2 .
the EBE. This indicator is a global evaluation that does not give any information on any specific patient located in a specific area of the green curve in Fig. 1. We also provide estimates of the shrinkage for the EBE and for the estimator that we propose and determine the prediction of the shrinkage based on the method proposed by Combes et al. [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].

Finally, from this prediction, we deduce a prediction for IPK parameters, as (Eq. 8):

Once lambda λ and gamma (Γ λ ) are defined, they are used to improve predictions made for a given patient (Eq. 8). In a more detailed manner, when we have concentration(s) (Y 1 , … , Y n ) available for a patient P, simulations allow to compute λ and Γ λ . Next, Eq. ( 4) is solved using the con- centrations (Y 1 , … , Y n ) observed in the patient and with λ computed at the previous step instead of .

From this computation, we obtain ηCl,P λ and therefore Γ λ ηCl,P λ , which is the de-shrunk estimate of its clear- ance individual random effect. By taking into account the imprecision observed during the simulations, it is also possible to define the prediction interval from which the individual pharmacokinetic parameter is estimated.

For example, we estimate a Γ λ ηCl,P λ for patient P equal to 0.2. Figure 1 shows that the true median value for

(8) φ = h Γ λ η λ , θ, A .
the individual pharmacokinetic parameter is 0.18 and the prediction interval is between 0.05 and 0.29.

Simulations Study

Datasets were simulated using R version 3.5.2.

This methodology was applied to two different models from the literature. The performance evaluation of the proposed methodology was performed on 5000 subjects simulated by the models described in the corresponding articles.

Usually, it is necessary to perform K simulations for each subject; a subject being defined by its sampling times t j and its covariate values A. Since we simulated patients with exactly the same sample times and covariate values, K = 5000 simulations could be used for all simulated patients.

For the rest of this work, we focus on clearance that controls exposure and use only one concentration to estimate individual clearance as per hospital practice.

For each drug tested, the sampling times were chosen in line with the model itself (i.e. sampling times used for the model construction) and according to a clinical practicability allowing a later clinical use of this methodology.

First Model

The first model used relates to the pharmacokinetics of iohexol in dogs [START_REF] Baklouti | Population pharmacokinetic model of iohexol in dogs to estimate glomerular filtration rate and optimize sampling time[END_REF]. This molecule allows establishment of the glomerular filtration rate (GFR). As it is an exogenous molecule eliminated mainly by glomerular filtration, its clearance represents GFR. The model was built from iohexol pharmacokinetic data in 49 dogs. Each kinetic profile includes five points (i.e. 5, 15, 60, 90, and 180 min).

The PPK model is a two-compartment model parameterised with clearances and volumes. This implies the f function of the first equation. The g function was chosen as g(t, , b) = bf (t, ).

Part of the clearance variability is explained by two covariates: the health status and the plasma creatinine concentration. Health status is a dichotomous covariate defined on the basis of the International Renal Interest Society (IRIS). The plasma creatinine concentration is a continuous covariate.

This results in the h function of Eq. ( 1) being such that ln( ) = ln(h( , , A))

= ln (Cl, V c ,V p ,Q) = ln( Cl ) + 1 1 [status] + θ 2 creat + Cl , ln( V c ) + V c , ln( V p ) + V p , ln Q .
Fig. 1 The y-axis represents the * Cl that have been used to simulate the concentrations; the x-axis contains the η * Cl λ obtained from the procedure described in section 2.1. The red line is the first bisector line. For a patient P, having a clearance individual random effect, denoted by * Cl,P , far below or above zero or, equivalently, an individual pharmacokinetic parameter far from the typical value, there is a systematic error. The systematic error is represented on the plot by the distance of the green curve from the red curve. The green curve represents Γ λ ηCl λ . The blue band is the 95% prediction interval. Its width, for a given ηCl,P λ (i.e., for a given patient) gives an evaluation of the imprecision of our method

The indicator function 1 [status] = 1 for a diseased individual and 0 otherwise (Table 1).

The variance-covariance matrix Ω of = Cl , V c , V p is assumed to be a diagonal matrix.

The dose chosen to perform the simulation is equal to 1600 mg. This corresponds to the average dose that was used to build the PPK model in the study by Baklouti et al. [START_REF] Baklouti | Population pharmacokinetic model of iohexol in dogs to estimate glomerular filtration rate and optimize sampling time[END_REF]. Similarly, we performed the simulations with sick dogs with a plasma creatinine level of 2.25 mg/dL that roughly corresponds to the average concentration reported in diseased animals [START_REF] Baklouti | Population pharmacokinetic model of iohexol in dogs to estimate glomerular filtration rate and optimize sampling time[END_REF].

We tested this methodology when a single sample is available at 180 min, corresponding to the best sampling time as identified in the study by Baklouti et al. [START_REF] Baklouti | Population pharmacokinetic model of iohexol in dogs to estimate glomerular filtration rate and optimize sampling time[END_REF].

Second Model

The second model relates the pharmacokinetics of isavuconazole [START_REF] Wu | Population pharmacokinetics of intravenous isavuconazole in solid-organ transplant recipients[END_REF]. This molecule is used in the treatment of invasive fungal infections. The model was built from 471 concentrations obtained from 79 patients.

The PPK model has two compartments. The error model for residual variability is a proportional error. Two covariates are included in the model, body mass index (BMI) and sex. BMI explains a part of peripheral distribution volume variability, and sex explains a part of clearance variability.

This results in the h function of Eq. ( 1) being such that

The indicator function 1 [sex] = 1 for a man and 0 for a woman (Table 2).

The dose in the simulation was 200 mg, which corresponds to the standard dose administered in hospital. The covariate value for the BMI was 25.4 kg/m 2 , corresponding to the average value found in the study [START_REF] Wu | Population pharmacokinetics of intravenous isavuconazole in solid-organ transplant recipients[END_REF]. The mean clearance selected was the clearance for men because the model was constructed with a majority of men.

The sampling time selected to perform simulations was 195 h because this is when steady state is reached and TDM can be proposed [START_REF] Bellmann | Pharmacokinetics of antifungal drugs: practical implications for optimized treatment of patients[END_REF].

Results

The results for simulated and predicted random effects on clearances are represented in Fig. 2. Figure 2a andb show the results for iohexol, and the Fig. 2c andd d). Therefore, the left plots correspond to the results before shrinkage correction, and the right plots show the results obtained after our corrections. For the two tested models, the λ value is different from 1 (iohexol λ = 0.017 ; isavucona- zole λ = 0.059 ), so η * Cl λ are different from the η * Cl (1) and therefore from the EBEs. The difference between the median curve (obtained using LOESS regression) in green and the first bisector line in red represents the shrinkage.

ln( ) = ln(h( , , A)) = ln (Cl, V c , V p , Q) = ln Cl × 1 -1 [sex] + 1 × 1 [sex] + Cl , ln( V c ), ln( V p + 2 × (BMI -25.38) + V p , ln Q .
Figure 2a andc shows that estimations of the IPK parameters provided by EBEs are not good: they are biased (i.e. the two clouds of points are not centred on the first bisector) and imprecise (i.e. the two clouds are thick).

Figure 2a (iohexol) shows that EBEs overestimated the IPK parameters, particularly for the lowest IPK parameter values. In this case, elimination clearance is overestimated and an unjustified dose increase is applied. Conversely, for isavuconazole, elimination clearance is underestimated for the lowest values, leading to an unjustified dose reduction. Thus, for drugs with a narrow therapeutic range, this kind of mistake can be detrimental to the patient; this means a therapeutic failure or, conversely, increased toxicity. Surprisingly, it can be seen that, for these two experimental designs (covariates, sampling times, PPK model), shrinkage does not affect the high IPK parameter values. However, these results could be radically different for another experimental design and evolve in an unpredictable way. As a consequence, for each experimental design (covariates, sampling times, PPK model), it is necessary to repeat the shrinkage correction.

The second important information observed in these plots is that the IPK parameter estimates are imprecise. Indeed, one can see that a large number of actual IPK parameters are consistent with a single estimated IPK parameter, i.e. the cloud of points is thick. If the role of λ is to reduce the thickness of the cloud of points, it seems that this decrease is modest in these two examples. Indeed, there is an incompressible limit to imprecision, mainly governed by the number of available concentrations for the patient. Moreover, the imprecision varies according to the values of η: the further η is from the population average, the greater the imprecision.

Figure 3 shows the improvement obtained with our method (iohexol in Fig. 3a and isavuconazole in Fig. 3b), with respect to EBEs, for different values of departure from median population clearance (i.e. * Cl ). The improvement is maximal for patients with extreme actual clearance (lowest and highest values). Conversely, our method does not improve the clearance estimation for patients with actual clearances close to the median population clearance; this is Fig. 3 Graphical representation of the evolution of "Improv" of iohexol and isavuconazole. Improvement is much more important for values far from the average than for values close to the mean. The improvement of the estimates following the correction of the empirical Bayes estimate is therefore not always the same. So, this meth-odology does not provide a great correction for patients who have a pharmacokinetic profile close to the average. However, for pharmacokinetic profiles far from average pharmacokinetic profiles (i.e. patients for whom therapeutic drug monitoring is usually indicated), this methodology greatly improves results logical. Finally, one would expect this improvement to be the same for extreme clearance values (small and large). This is true for isavuconazole but surprisingly is not the case for iohexol. To date, we have no explanation for this phenomenon, and it most likely depends on the experimental design used. The overall improvement, as defined by Eq. ( 7), is 21.2% for iohexol and 15.5% for isavuconazole.

Table 3 presents the results of shrinkage (for the EBE or the proposed estimator) and the prediction of shrinkage provided by the Combes et al. [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] method. The proposed estimator has a lower shrinkage than the EBE, whatever the application (i.e. the PPK model, the sampling times, etc.). Surprisingly, the method proposed by Combes et al. [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] does not provide a good prediction of EBE shrinkage, especially for isavuconazole. We proposed a method to improve the prediction of IPK parameters using the measured concentration for one patient.

The EBE is the traditional predictor; however, this individual prediction suffers from an "attraction" toward the typical population value and moves the prediction away from the actual value of the IPK parameters. This "distance" is called shrinkage [START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF]. It is usually computed over both different individuals (i.e. with different IPK parameters) and different sampling times. This number gives information on the average departure from the actual value (i.e. actual clearance, volume of distribution, etc.), but it cannot be used as it stands to correct the prediction of the IPK parameters for a particular patient [START_REF] Combes | Powers of the likelihood ratio test and the correlation test using empirical bayes estimates for various shrinkages in population pharmacokinetics[END_REF][START_REF] Xu | Further evaluation of covariate analysis using empirical bayes estimates in population pharmacokinetics: the perception of shrinkage and likelihood ratio test[END_REF].

In contrast, the method we propose allows a specific shrinkage correction for each patient by considering the sampling times, the PPK model, and the patient's covariate values. Moreover, for each patient, the level of shrinkage correction depends on the proximity of the patient's pharmacokinetic parameters to the population's pharmacokinetic parameters. Indeed, we have shown that the correction was at its maximum for extreme values. This information is particularly important as, in hospital practice, when TDM is requested, the pharmacokinetic parameters of the patients are quite extreme values (i.e. far from the mean population values). This is the case for example for patients in the intensive care unit and those with cystic fibrosis or obesity, among others [START_REF] Sutter | Pharmacokinetics in patients with cystic fibrosis: a systematic review of data published between 1999 and 2019[END_REF][START_REF] Stewart | Antibiotic use in critical illness[END_REF][START_REF] Smit | Obesity and drug pharmacology: a review of the influence of obesity on pharmacokinetic and pharmacodynamic parameters[END_REF]. Thus, this approach allows us to expect a significant improvement in the estimates of the IPK parameters. In addition to TDM, a better estimate of EBEs could help to improve PPK models during their construction and validation steps.

As stated in Eq. ( 6), was chosen to minimise the shrinkage on clearance. Thus, nothing guarantees a shrinkage decrease of other individual pharmacokinetic parameters. This is a not problem because, for drugs with a linear kinetic, clearance drives exposure and then TDM. We estimated this pharmacokinetic parameter with only one concentration both to mimic hospital practice and to guarantee a high shrinkage to be removed. This choice does not exclude the application of the same kind of methodology for all the other pharmacokinetic parameters (i.e. to choose a and Γ function spe- cific to the pharmacokinetic parameter of interest).

In this work, the method of Combes et al. [START_REF] Baklouti | Population pharmacokinetic model of iohexol in dogs to estimate glomerular filtration rate and optimize sampling time[END_REF] poorly predicted EBE shrinkage. We consider that this is mainly due to the assumed asymptotic framework. It is clear that, with only one observation per individual, the variance of P(Ψ|Y) can be quite different from the corresponding Bayesian Fisher information matrix obtained by expanding ln P(Ψ|Y) at its mode. Actually, the quality of approximation of P(Ψ|Y) by a Gaussian distribution is difficult to anticipate.

One of the limits to our methodology is the need for a validated PPK model. Indeed, all these results have been obtained conditional on θ, b, Ω . In other words, an important assumption behind the method is that these estimations are precise enough to be used as if they were the actual value of the parameters. This hopefully happens for estimates obtained for a validated model built with a large number of individuals. If not, a model misspecification could lead to biased predictions. In this case, the shrinkage correction could be negligible compared with the bias. As such, the impact of a model misspecification on the final individual parameter estimate should be investigated.

However, although validation of PPK models for TDM is beyond the scope of this paper, the implicit assumption is that the chosen model properly describes the kinetic profile of the hospitalised patient. From our experience, this is far from the case. An easy way to identify a poor PPK model is to look at the residual terms. The most often used residual error models are the proportional error model (i.e. g(t, , b) = b 1 f (t, ) ) and in some cases a combined error model (i.e. g(t, , b) = b 2 + b 1 f (t, ) ). A large b 1 (e.g. > 40%) indicates a poor PPK model and thus precludes any TDM based on IPK parameter estimations.

To conclude, this innovative methodology is promising since, to the best of our knowledge, no individual shrinkage correction has been proposed and it consistently provided better results than EBEs.

Declarations

Funding No external funding was used in the preparation of this manuscript.

Conflicts of interest Sarah Baklouti, Peggy Gandia, and Didier Concordet have no potential conflicts of interest that might be relevant to the contents of this manuscript.

Ethics approval Not applicable.

  are for isavuconazole. On all plots, the vertical axis corresponds to * Cl . The horizontal axis represents the EBEs (i.e. η * Cl (1) ) on the left plots (a and c) and η * Cl λ on the right plots (b and

Fig. 2

 2 Fig. 2 Graphical representation of * Cl as a function of η * Cl (1) (i.e. empirical Bayes estimate for clearance) (left) and of * Cl as a function of Γ λ η * Cl λ (right). The red curves represent the first bisector, the

Table 1

 1 Population pharmacokinetic parameters of iohexol published by Baklouti et al.[START_REF] Baklouti | Population pharmacokinetic model of iohexol in dogs to estimate glomerular filtration rate and optimize sampling time[END_REF]

			Table 2 Population pharmacokinetics model of	Parameters	Value
	Parameters θCl θV c θV p θQ Covariates θ1 θ2 Standard deviation of = Cl , V p , V c ωCl ωV c	Value 0.00432 L/min/kg 0.163 L/kg 0.0584 L/kg 0.00336 L/min/kg -0.345 -0.517 dL/mg 0.207 0.249	isavuconazole published by Wu et al. [19]	θCl θV c θV p θQ Covariates θ1 (sex on clearance) 2.73 L/h 4.28 L/h 57.6 L 468 L 37.4 L/h θ2 (BMI on V p ) 39.7 L Standard deviation of = Cl , V p ωCl 0.461 ωV p 0.444 Residual error Proportional 20.5%
	ωV p Residual error Proportional	0.213 6.12%		BMI body mass index, V p peripheral volume of distribu-tion

Table 3

 3 Shrinkage estimates calculated from empirical Bayes estimate and the proposed estimator.

		Empirical Bayes	Proposed esti-	Shrinkage a (%)
		estimate (%)	mator (%)
	Iohexol	40	20	24
	Isavuconazole	57	49	9.3
	a Shrinkage of the proposed estimator computed using the Combes
	et al. [14] method

Availability of data and material Not applicable.

Consent to participate Not applicable.

Consent for publication Not applicable.

Code availability

The code is available from the corresponding author on request.

Author contributions DC and SB proposed the methodology to correct shrinkage, and SB performed this methodology with iohexol and isavuconazole. DC, PG, and SB drafted the paper. All authors critically reviewed several drafts and approved the final manuscript.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License, which permits any non-commercial use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative "De-Shrinking" EBEs Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by-nc/4. 0/.