Uncertainty Computation at Finite Distance in Nonlinear Mixed Effects Models-a New Method Based on Metropolis-Hastings Algorithm.

The standard errors (SE) of the maximum likelihood estimates (MLE) of the population parameter vector in nonlinear mixed effect models (NLMEM) are usually estimated using the inverse of the Fisher information matrix (FIM). However, at a finite distance, i.e. far from the asymptotic, the FIM can underestimate the SE of NLMEM parameters. Alternatively, the standard deviation of the posterior distribution, obtained in Stan via the Hamiltonian Monte Carlo algorithm, has been shown to be a proxy for the SE, since, under some regularity conditions on the prior, the limiting distributions of the MLE and of the maximum a posterior estimator in a Bayesian framework are equivalent. In this work, we develop a similar method using the Metropolis-Hastings (MH) algorithm in parallel to the stochastic approximation expectation maximisation (SAEM) algorithm, implemented in the saemix R package. We assess this method on different simulation scenarios and data from a real case study, comparing it to other SE computation methods. The simulation study shows that our method improves the results obtained with frequentist methods at finite distance. However, it performed poorly in a scenario with the high variability and correlations observed in the real case study, stressing the need for calibration.

Prediction of Individual Disease Progression Including Parameter Uncertainty in Rare Neurodegenerative Diseases: The Example of Autosomal-Recessive Spastic Ataxia Charlevoix Saguenay (ARSACS).

The aim of this study was to develop a model to predict individual subject disease trajectories including parameter uncertainty and accounting for missing data in rare neurological diseases, showcased by the ultra-rare disease Autosomal-Recessive Spastic Ataxia Charlevoix Saguenay (ARSACS). We modelled the change in SARA (Scale for Assessment and Rating of Ataxia) score versus Time Since Onset of symptoms using non-linear mixed effect models for a population of 173 patients with ARSACS included in the prospective real-world multicenter Autosomal Recessive Cerebellar Ataxia (ARCA) registry. We used the Multivariate Imputation Chained Equation (MICE) algorithm to impute missing covariates, and a covariate selection procedure with a pooled p-value to account for the multiply imputed data sets. We then investigated the impact of covariates and population parameter uncertainty on the prediction of the individual trajectories up to 5 years after their last visit. A four-parameter logistic function was selected. Men were estimated to have a 25% lower SARA score at disease onset and a moderately higher maximum SARA score, and time to progression (T50) was estimated to be 35% lower in patients with age of onset over 15 years. The population disease progression rate started slowly at 0.1 points per year peaking to a maximum of 0.8 points per year (at 36.8 years since onset of symptoms). The prediction intervals for SARA scores 5 years after the last visit were large (median 7.4 points, Q1-Q3: 6.4-8.5); their size was mostly driven by individual parameter uncertainty and individual disease progression rate at that time.

Joint modeling of tumor dynamics and progression-free survival in advanced breast cancer: Leveraging data from amcenestrant early phase I-II trials.

A joint modeling framework was developed using data from 75 patients of early amcenestrant phase I-II AMEERA-1-2 dose escalation and expansion cohorts. A semi-mechanistic tumor growth inhibition (TGI) model was developed. It accounts for the dynamics of sensitive and resistant tumor cells, an exposure-driven effect on tumor proliferation of sensitive cells, and a delay in the initiation of treatment effect to describe the time course of target lesion tumor size (TS) data. Individual treatment exposure overtime was introduced in the model using concentrations predicted by a population pharmacokinetic model of amcenestrant. This joint modeling framework integrated complex RECISTv1.1 criteria information, linked TS metrics to progression-free survival (PFS), and was externally evaluated using the randomized phase II trial AMEERA-3. We demonstrated that the instantaneous rate of change in TS (TS slope) was an important predictor of PFS and the developed joint model was able to predict well the PFS of amcenestrant phase II monotherapy trial using only early phase I-II data. This provides a good modeling and simulation tool to inform early development decisions.

Bayesian framework for multi-source data integration-Application to human extrapolation from preclinical studies.

In preclinical investigations, for example, in in vitro, in vivo, and in silico studies, the pharmacokinetic, pharmacodynamic, and toxicological characteristics of a drug are evaluated before advancing to first-in-man trial. Usually, each study is analyzed independently and the human dose range does not leverage the knowledge gained from all studies. Taking into account all preclinical data through inferential procedures can be particularly interesting in obtaining a more precise and reliable starting dose and dose range. Our objective is to propose a Bayesian framework for multi-source data integration, customizable, and tailored to the specific research question. We focused on preclinical results extrapolated to humans, which allowed us to predict the quantities of interest (e.g. maximum tolerated dose, etc.) in humans. We build an approach, divided into four steps, based on a sequential parameter estimation for each study, extrapolation to human, commensurability checking between posterior distributions and final information merging to increase the precision of estimation. The new framework is evaluated via an extensive simulation study, based on a real-life example in oncology. Our approach allows us to better use all the information compared to a standard framework, reducing uncertainty in the predictions and potentially leading to a more efficient dose selection.

Extending the code in the open-source saemix package to fit joint models of longitudinal and time-to-event data.

BACKGROUND AND OBJECTIVE: Joint modeling of longitudinal and time-to-event data has gained attention over recent years with extensive developments including nonlinear models for longitudinal outcomes and flexible time-to-event models for survival outcomes, possibly involving competing risks. However, in popular software such as R, the function used to describe the biomarker dynamic is mainly linear in the parameters, and the survival submodel relies on pre-implemented functions (exponential, Weibull, ...). The objective of this work is to extend the code from the saemix package (version 3.1 on CRAN) to fit parametric joint models where longitudinal submodels are not necessary linear in their parameters, with full user control over the model function. METHODS: We used the saemix package, designed to fit nonlinear mixed-effects models (NLMEM) through the Stochastic Approximation Expectation Maximization (SAEM) algorithm, and extended the main functions to joint model estimation. To compute standard errors (SE) of parameter estimates, we implemented a recently developed stochastic algorithm. A simulation study was proposed to assess (i) the performances of parameter estimation, (ii) the SE computation and (iii) the type I error when testing independence between the two submodels. Four joint models were considered in the simulation study, combining a linear or nonlinear mixed-effects model for the longitudinal submodel, with a single terminal event or a competing risk model. RESULTS: For all simulation scenarios, parameters were precisely and accurately estimated with low bias and uncertainty. For complex joint models (with NLMEM), increasing the number of chains of the algorithm was necessary to reduce bias, but earlier censoring in the competing risk scenario still challenged the estimation. The empirical SE of parameters obtained over all simulations were very close to those computed with the stochastic algorithm. For more complex joint models (involving NLMEM), some estimates of random effects variances had higher uncertainty and their SE were moderately under-estimated. Finally, type I error was controlled for each joint model. CONCLUSIONS: saemix is a flexible open-source package and we adapted it to fit complex parametric joint models that may not be estimated using standard tools. Code and examples to help users get started are freely available on Github.

Therapeutic drug monitoring and virological response at week 48 in a cohort of HIV-1-infected patients switching to dolutegravir/rilpivirine dual maintenance therapy (ANRS-MIE-BIRIDER study).

AIMS: Dolutegravir (DTG) and rilpivirine (RPV) dual therapy is now recommended as a switch option in virologically suppressed HIV patients. Literature suggests that virological failure with dual therapy could possibly relate to subtherapeutic drug concentrations. In this study, we aimed at describing the DTG and RPV trough plasma concentrations (Cmin) and plasma HIV-1 RNA viral load (VL) during maintenance dual therapy. METHODS: We performed a retrospective analysis of DTG and RPV therapeutic drug monitoring in people living with HIV/AIDS (PLWHA) with dual therapy in 9 French centres. DTG and RPV trough plasma concentrations were estimated using a Bayesian approach to predict Cmin. The relationship between the pharmacokinetics of DTG and RPV and VL > 50 copies (cp)/mL was explored using joint nonlinear mixed models. The frequency of subtherapeutic threshold (DTG Cmin below 640 ng/mL and RPV Cmin below 50 ng/mL) were compared between PLWHA presenting VL > 50 cp/mL or not during the study. RESULTS: At baseline, 209 PLWHA were enrolled in the study. At week 48, 19 people living with HIV/AIDS (9.1%) discontinued their treatment and 15 PLWHA (7.1%) exhibited VL > 50 cp/mL. Six PLWHA out of 15 (40.0%) with VL > 50 cp/mL during the follow-up had at least 1 Cmin below the respective thresholds while only 26/194 patients (13.4%) without virological replication had at least 1 concentration below the threshold (P = .015). CONCLUSION: A majority of PLWHA receiving DTG/RPV maintenance dual therapy demonstrated VL < 50 cp/mL but virological replication was more frequent in people living with HIV/AIDS with subtherapeutic Cmin.

Multivariate Exact Discrepancy: A New Tool for PK/PD Model Evaluation.

BACKGROUND: Pharmacokinetic models are evaluated using three types of metrics: those based on estimating the typical pharmacokinetic parameters, those based on predicting individual pharmacokinetic parameters and those that compare data and model distributions. In the third groups of metrics, the best-known methods are Visual Predictive Check (VPC) and Normalised Prediction Distribution Error (NPDE). Despite their usefulness, these methods have some limitations, especially for the analysis of dependent concentrations, i.e., evaluated in the same patient. OBJECTIVE: In this work, we propose an evaluation method that accounts for the dependency between concentrations. METHODS: Thanks to the study of the distribution of simulated vectors of concentrations, the method provides one probability per individual that its observations (i.e., concentrations) come from the studied model. The higher the probability, the better the model fits the individual. By examining the distribution of these probabilities for a set of individuals, we can evaluate the model as a whole. RESULTS: We demonstrate the effectiveness of our method through two examples. Our approach successfully detects misspecification in the structural model and identifies outlier kinetics in a set of kinetics. CONCLUSION: We propose a straightforward method for evaluating models during their development and selecting a model to perform therapeutic drug monitoring. Based on our preliminary results, the method is very promising but needs to be validated on a larger scale.

Impact of model misspecification on model-based tests in PK studies with parallel design: real case and simulation studies.

This article evaluates the performance of pharmacokinetic (PK) equivalence testing between two formulations of a drug through the Two-One Sided Tests (TOST) by a model-based approach (MB-TOST), as an alternative to the classical non-compartmental approach (NCA-TOST), for a sparse design with a few time points per subject. We focused on the impact of model misspecification and the relevance of model selection for the reference data. We first analysed PK data from phase I studies of gantenerumab, a monoclonal antibody for the treatment of Alzheimer's disease. Using the original rich sample data, we compared MB-TOST to NCA-TOST for validation. Then, the analysis was repeated on a sparse subset of the original data with MB-TOST. This analysis inspired a simulation study with rich and sparse designs. With rich designs, we compared NCA-TOST and MB-TOST in terms of type I error and study power. With both designs, we explored the impact of misspecifying the model on the performance of MB-TOST and adding a model selection step. Using the observed data, the results of both approaches were in general concordance. MB-TOST results were robust with sparse designs when the underlying PK structural model was correctly specified. Using the simulated data with a rich design, the type I error of NCA-TOST was close to the nominal level. When using the simulated model, the type I error of MB-TOST was controlled on rich and sparse designs, but using a misspecified model led to inflated type I errors. Adding a model selection step on the reference data reduced the inflation. MB-TOST appears as a robust alternative to NCA-TOST, provided that the PK model is correctly specified and the test drug has the same PK structural model as the reference drug.

Effects of Postpartum Supplemental Oral Ca for Dairy Cows Fed Prepartum Dietary Acidogenic Salts.

Postpartum hypocalcemia is a problem in dairy cows. Both the Jersey vs. Holstein breed and increasing parity are known risk factors. Our objectives were: (1) to evaluate a simple approach to provide dietary acidogenic salts suitable for application on small dairies and (2) to evaluate the combined effects of degree of acidification and oral Ca supplementation along with breed and parity group on periparturient Ca status of Holstein and Jersey cows. Cows were moved weekly from the far-off dry pen at 260 days pregnant to the close-up pen, where all cows received the acidogenic diets. The diet was offered as a total mixed ration and CaCl2, and our source of acidogenic salts was top-dressed in liquid form and mixed in by hand. Thirty-six cows were blocked by parity group (parity = 2 vs. parity ≥ 3) and breed (Holstein vs. Jersey) and assigned to one of two treatments (no intervention or postpartum oral Ca bolus supplementation) in an alternating fashion, based on expected date of parturition. Urinary acidification appeared complete within 3-4 days. Increased urinary Ca excretion was >93% of maximum from 7-21 days before falling to <5% of maximum by 28 days. Serum Ca concentrations 12-24 h postpartum were lower for Jerseys vs. Holsteins and for parity ≥ 3 vs. parity = 2 cows. Serum Ca over 6-48 h postpartum decreased and increased, respectively, with oral Ca supplementation for parity = 2 and parity ≥ 3 cows. Decreased prepartum urinary Ca excretion and increased colostrum yield appear to be independent risk factors of hypocalcemia for parity ≥ 3 Jerseys.

Ciprofloxacin population pharmacokinetics during long-term treatment of osteoarticular infections.

BACKGROUND: Ciprofloxacin is an antibiotic used in osteoarticular infections owing to its very good bone penetration. Very few pharmacokinetic data are available in this population. OBJECTIVES: To investigate oral ciprofloxacin population pharmacokinetics in adult patients treated for osteoarticular infections and propose guidance for more effective dosing. METHODS: A retrospective population-pharmacokinetic analysis was performed on 92 consecutive hospitalized patients in the orthopaedic department. Ciprofloxacin plasma samples were obtained on one or two occasions during treatment. Plasma concentration was measured using ultra-performance liquid chromatography system coupled with tandem mass spectrometry. Data analysis was performed using a non-linear mixed-effect approach via Monolix 2019R2. RESULTS: A total of 397 plasma samples were obtained with 11.5% and 41.6% of patients being below the therapeutic target for Gram-negative and staphylococcal infections, respectively. Ciprofloxacin pharmacokinetics were best described by a two-compartment model with a first-order absorption. Ciprofloxacin apparent plasma clearances and volumes of distribution were dependent on patients' fat-free mass according to the allometric rule. Elimination clearance was also positively related to renal function through the modification of diet in renal disease equation (MDRD) and rifampicin co-administration. When patients are co-treated with rifampicin, ciprofloxacin dosage should be increased by 50% to 60%. CONCLUSIONS: This study showed that free-fat mass was a better size predictor than total body weight for ciprofloxacin clearance and volumes terms. Moreover, both MDRD and rifampicin status were significant predictors of individual ciprofloxacin clearance. Our study suggests that individual adjustment of ciprofloxacin dose in osteoarticular infections with less-susceptible bacteria might be indicated to reach required efficacy targets.

Development of a dosing-adjustment tool for fluoroquinolones in osteoarticular infections: The Fluo-pop study.

Fluoroquinolones efficacy depend on both the drug exposure and the level of drug resistance of the bacteria responsible for the infection. Specifically for the Staphylococcus species, which is the microorganism mainly involved in osteoarticular infections (OAI), in-vitro data reported that an AUC/MIC ratio above 115 h maximizes drug efficacy. However, data on OAI patients are lacking and a simple approach to access AUCs is still a clinical issue. We conducted a prospective, single-center study in 30 OAI patients hospitalized in the Rennes University Hospital to model ofloxacin pharmacokinetics and to define a limited sampling strategy (LSS) suitable for ofloxacin and levofloxacin treatments. Modeling was conducted with the Monolix software. The final model was externally validated using levofloxacin data. Monte-Carlo simulations were used to evaluate the probability of target attainment (PTA) of different dosing regimens. Two hundred and ninety-seven (297) ofloxacin concentrations were available for the pharmacokinetic modeling. Ofloxacin pharmacokinetics was best described using a bicompartmental model with a first order elimination, and a transit compartment model absorption. CKD-EPI and sex explained half of ofloxacin pharmacokinetic variability. For LSS, the 0, 1 h and 3 h sampling scheme resulted in the best approach both for BID and TID dosages (R2 adjusted = 91.1% and 95.0%, outliers = 4.8% and 5.0%, respectively). PTA allows choosing the best drug and dosage according to various hypotheses. A simple 3-sample protocol (pre-dose, 1 h after intake and 3 h after intake) to estimate ofloxacin and levofloxacin AUC allows optimal drug dosage for the treatment of osteoarticular infections.

Conditional Non-parametric Bootstrap for Non-linear Mixed Effect Models.

PURPOSE: Non-linear mixed effect models are widely used and increasingly integrated into decision-making processes. Propagating uncertainty is an important element of this process, and while standard errors (SE) on pa- rameters are most often computed using asymptotic approaches, alternative methods such as the bootstrap are also available. In this article, we propose a modified residual parametric bootstrap taking into account the different levels of variability involved in these models. METHODS: The proposed approach uses samples from the individual conditional distribution, and was implemented in R using the saemix algorithm. We performed a simulation study to assess its performance in different scenarios, comparing it to the asymptotic approximation and to standard bootstraps in terms of coverage, also looking at bias in the parameters and their SE. RESULTS: Simulations with an Emax model with different designs and sigmoidicity factors showed a similar coverage rate to the parametric bootstrap, while requiring less hypotheses. Bootstrap improved coverage in several scenarios compared to the asymptotic method especially for the variance param-eters. However, all bootstraps were sensitive to estimation bias in the original datasets. CONCLUSIONS: The conditional bootstrap provided better coverage rate than the traditional residual bootstrap, while preserving the structure of the data generating process.

Developing Tools to Evaluate Non-linear Mixed Effect Models: 20 Years on the npde Adventure.

This article revisits 20 years of our work in developing evaluation tools adapted to non-linear mixed effect models. These hierarchical models involve a large number of assumptions concerning the structural evolution of the outcomes, the link between different outcomes, the variabilities in the parameters and model evaluation aims at assessing these various components, both to help guide the model building and to communicate on model adequacy for a given purpose. During our career, we have developed and extended simulation-based evaluation tools called normalised prediction discrepancies (npd) and normalised prediction distribution errors (npde), providing informative diagnostics through graphs and tests.

The Use of Translational Modelling and Simulation to Develop Immunomodulatory Therapy as an Adjunct to Antibiotic Treatment in the Context of Pneumonia.

The treatment of respiratory tract infections is threatened by the emergence of bacterial resistance. Immunomodulatory drugs, which enhance airway innate immune defenses, may improve therapeutic outcome. In this concept paper, we aim to highlight the utility of pharmacometrics and Bayesian inference in the development of immunomodulatory therapeutic agents as an adjunct to antibiotics in the context of pneumonia. For this, two case studies of translational modelling and simulation frameworks are introduced for these types of drugs up to clinical use. First, we evaluate the pharmacokinetic/pharmacodynamic relationship of an experimental combination of amoxicillin and a TLR4 agonist, monophosphoryl lipid A, by developing a pharmacometric model accounting for interaction and potential translation to humans. Capitalizing on this knowledge and associating clinical trial extrapolation and statistical modelling approaches, we then investigate the TLR5 agonist flagellin. The resulting workflow combines expert and prior knowledge on the compound with the in vitro and in vivo data generated during exploratory studies in order to construct high-dimensional models considering the pharmacokinetics and pharmacodynamics of the compound. This workflow can be used to refine preclinical experiments, estimate the best doses for human studies, and create an adaptive knowledge-based design for the next phases of clinical development.

Timing of antiviral treatment initiation is critical to reduce SARS-CoV-2 viral load.

We modeled the viral dynamics of 13 untreated patients infected with SARS-CoV-2 to infer viral growth parameters and predict the effects of antiviral treatments. In order to reduce peak viral load by more than 2 logs, drug efficacy needs to be greater than 80% if treatment is administered after symptom onset; an efficacy of 50% could be sufficient if treatment is initiated before symptom onset. Given their pharmacokinetic/pharmacodynamic properties, current investigated drugs may be in a range of 20-70% efficacy. They may help control virus if administered very early, but may not have a major effect in severe patients.

Timing of Antiviral Treatment Initiation is Critical to Reduce SARS-CoV-2 Viral Load.

We modeled the viral dynamics of 13 untreated patients infected with severe acute respiratory syndrome-coronavirus 2 to infer viral growth parameters and predict the effects of antiviral treatments. In order to reduce peak viral load by more than two logs, drug efficacy needs to be > 90% if treatment is administered after symptom onset; an efficacy of 60% could be sufficient if treatment is initiated before symptom onset. Given their pharmacokinetic/pharmacodynamic properties, current investigated drugs may be in a range of 6-87% efficacy. They may help control virus if administered very early, but may not have a major effect in severely ill patients.

Application of Item Response Theory to Model Disease Progression and Agomelatine Effect in Patients with Major Depressive Disorder.

INTRODUCTION: In this paper, we studied the effect over time of agomelatine, an antidepressant drug administered in patient with major depressive disorder, through item response theory (IRT), taking into account a strong placebo effect and missing not at random. We also assessed the informativeness of the HAMD-17 scale's item. MATERIALS AND METHODS: The data includes five phase III clinical trials sponsored by Servier Institute, totalling 1549 patients followed during a maximum of 1 year. At each observation, individual scores for the 17 items of the HAMD scale were recorded. The probability for each score was modelled with IRT. A non-linear mixed effects model was used to describe the evolution of the disease and was coupled with a time to event model to predict dropout. Clinical trial simulations were then used to compare placebo and active treatment. Informativeness of each item was evaluated using the Fisher information theory. RESULTS: The best model combined an IRT model, a longitudinal model for underlying depression which describes the remission and then a possible relapse, and a hazard model for dropout depending on the evolution from baseline. The drug effect was best modelled as an effect on the remission and the relapse phases. The median predicted drop in HAMD between baseline and 6 weeks was 8.8 (90% PI, 8.3-9.2) when on placebo and 13.1 (90% PI, 12.8-13.4) when treated. Nine items were found to be the most informative. CONCLUSION: The IRT framework allowed to characterise the evolution of depression with time and estimate the effect of agomelatine, as well as the link between symptoms and disease.

Recent advances in methodology for clinical trials in small populations: the InSPiRe project.

Where there are a limited number of patients, such as in a rare disease, clinical trials in these small populations present several challenges, including statistical issues. This led to an EU FP7 call for proposals in 2013. One of the three projects funded was the Innovative Methodology for Small Populations Research (InSPiRe) project. This paper summarizes the main results of the project, which was completed in 2017.The InSPiRe project has led to development of novel statistical methodology for clinical trials in small populations in four areas. We have explored new decision-making methods for small population clinical trials using a Bayesian decision-theoretic framework to compare costs with potential benefits, developed approaches for targeted treatment trials, enabling simultaneous identification of subgroups and confirmation of treatment effect for these patients, worked on early phase clinical trial design and on extrapolation from adult to pediatric studies, developing methods to enable use of pharmacokinetics and pharmacodynamics data, and also developed improved robust meta-analysis methods for a small number of trials to support the planning, analysis and interpretation of a trial as well as enabling extrapolation between patient groups. In addition to scientific publications, we have contributed to regulatory guidance and produced free software in order to facilitate implementation of the novel methods.

PFIM 4.0, an extended R program for design evaluation and optimization in nonlinear mixed-effect models.

BACKGROUND AND OBJECTIVE: Nonlinear mixed-effect models (NLMEMs) are increasingly used for the analysis of longitudinal studies during drug development. When designing these studies, the expected Fisher information matrix (FIM) can be used instead of performing time-consuming clinical trial simulations. The function PFIM is the first tool for design evaluation and optimization that has been developed in R. In this article, we present an extended version, PFIM 4.0, which includes several new features. METHODS: Compared with version 3.0, PFIM 4.0 includes a more complete pharmacokinetic/pharmacodynamic library of models and accommodates models including additional random effects for inter-occasion variability as well as discrete covariates. A new input method has been added to specify user-defined models through an R function. Optimization can be performed assuming some fixed parameters or some fixed sampling times. New outputs have been added regarding the FIM such as eigenvalues, conditional numbers, and the option of saving the matrix obtained after evaluation or optimization. Previously obtained results, which are summarized in a FIM, can be taken into account in evaluation or optimization of one-group protocols. This feature enables the use of PFIM for adaptive designs. The Bayesian individual FIM has been implemented, taking into account a priori distribution of random effects. Designs for maximum a posteriori Bayesian estimation of individual parameters can now be evaluated or optimized and the predicted shrinkage is also reported. It is also possible to visualize the graphs of the model and the sensitivity functions without performing evaluation or optimization. RESULTS: The usefulness of these approaches and the simplicity of use of PFIM 4.0 are illustrated by two examples: (i) an example of designing a population pharmacokinetic study accounting for previous results, which highlights the advantage of adaptive designs; (ii) an example of Bayesian individual design optimization for a pharmacodynamic study, showing that the Bayesian individual FIM can be a useful tool in therapeutic drug monitoring, allowing efficient prediction of estimation precision and shrinkage for individual parameters. CONCLUSION: PFIM 4.0 is a useful tool for design evaluation and optimization of longitudinal studies in pharmacometrics and is freely available at http://www.pfim.biostat.fr.