Model-based bioequivalence approach for sparse pharmacokinetic bioequivalence studies: Model selection or model averaging?

Conventional pharmacokinetic (PK) bioequivalence (BE) studies aim to compare the rate and extent of drug absorption from a test (T) and reference (R) product using non-compartmental analysis (NCA) and the two one-sided test (TOST). Recently published regulatory guidance recommends alternative model-based (MB) approaches for BE assessment when NCA is challenging, as for long-acting injectables and products which require sparse PK sampling. However, our previous research on MB-TOST approaches showed that model misspecification can lead to inflated type I error. The objective of this research was to compare the performance of model selection (MS) on R product arm data and model averaging (MA) from a pool of candidate structural PK models in MBBE studies with sparse sampling. Our simulation study was inspired by a real case BE study using a two-way crossover design. PK data were simulated using three structural models under the null hypothesis and one model under the alternative hypothesis. MB-TOST was applied either using each of the five candidate models or following MS and MA with or without the simulated model in the pool. Assuming T and R have the same PK model, our simulation shows that following MS and MA, MB-TOST controls type I error rates at or below 0.05 and attains similar or even higher power than when using the simulated model. Thus, we propose to use MS prior to MB-TOST for BE studies with sparse PK sampling and to consider MA when candidate models have similar Akaike information criterion.

Joint modeling of monocyte HLA-DR expression trajectories predicts 28-day mortality in severe SARS-CoV-2 patients.

The recent SarsCov2 pandemic has disrupted healthcare system notably impacting intensive care units (ICU). In severe cases, the immune system is dysregulated, associating signs of hyperinflammation and immunosuppression. In the present work, we investigated, using a joint modeling approach, whether the trajectories of cellular immunological parameters were associated with survival of COVID-19 ICU patients. This study is based on the REA-IMMUNO-COVID cohort including 538 COVID-19 patients admitted to ICU between March 2020 and May 2022. Measurements of monocyte HLA-DR expression (mHLA-DR), counts of neutrophils, of total lymphocytes, and of CD4+ and CD8+ subsets were performed five times during the first month after ICU admission. Univariate joint models combining survival at day 28 (D28), hospital discharge and longitudinal analysis of those biomarkers' kinetics with mixed-effects models were performed prior to the building of a multivariate joint model. We showed that a higher mHLA-DR value was associated with a lower risk of death. Predicted mHLA-DR nadir cutoff value that maximized the Youden index was 5414 Ab/C and led to an AUC = 0.70 confidence interval (95%CI) = [0.65; 0.75] regarding association with D28 mortality while dynamic predictions using mHLA-DR kinetics until D7, D12 and D20 showed AUCs of 0.82 [0.77; 0.87], 0.81 [0.75; 0.87] and 0.84 [0.75; 0.93]. Therefore, the final joint model provided adequate discrimination performances at D28 after collection of biomarker samples until D7, which improved as more samples were collected. After severe COVID-19, decreased mHLA-DR expression is associated with a greater risk of death at D28 independently of usual clinical confounders.

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.

A population pharmacokinetic model to guide clozapine dose selection, based on age, sex, ethnicity, body weight and smoking status.

AIMS: Guidance on clozapine dosing in treatment-resistant schizophrenia is based largely on data from White young adult males. This study aimed to investigate the pharmacokinetic profiles of clozapine and N-desmethylclozapine (norclozapine) across the age range, accounting for sex, ethnicity, smoking status and body weight. METHODS: A population pharmacokinetic model, implemented in Monolix, that linked plasma clozapine and norclozapine via a metabolic rate constant, was used to analyse data from a clozapine therapeutic drug monitoring service, 1993-2017. RESULTS: There were 17 787 measurements from 5960 patients (4315 male) aged 18-86 years. The estimated clozapine plasma clearance was reduced from 20.2 to 12.0 L h-1 between 20 and 80 years. Model-based dose predictions to attain a predose plasma clozapine concentration of 0.35 mg L-1 was 275 (90% prediction interval 125, 625) mg day-1 in nonsmoking, White males weighing 70 kg and aged 40 years. The corresponding predicted dose was increased by 30% in smokers, decreased by 18% in females, and was 10% higher and 14% lower in otherwise analogous Afro-Caribbean and Asian patients, respectively. Overall, the predicted dose decreased by 56% between 20 and 80 years. CONCLUSION: The large sample size and wide age range of the patients studied allowed precise estimation of dose requirements to attain predose clozapine concentration of 0.35 mg L-1 . The analysis was, however, limited by the absence of data on clinical outcome and future studies are required to determine optimal predose concentrations specifically in those aged over 65 years.

Impact of variants of concern on SARS-CoV-2 viral dynamics in non-human primates.

The impact of variants of concern (VoC) on SARS-CoV-2 viral dynamics remains poorly understood and essentially relies on observational studies subject to various sorts of biases. In contrast, experimental models of infection constitute a powerful model to perform controlled comparisons of the viral dynamics observed with VoC and better quantify how VoC escape from the immune response. Here we used molecular and infectious viral load of 78 cynomolgus macaques to characterize in detail the effects of VoC on viral dynamics. We first developed a mathematical model that recapitulate the observed dynamics, and we found that the best model describing the data assumed a rapid antigen-dependent stimulation of the immune response leading to a rapid reduction of viral infectivity. When compared with the historical variant, all VoC except beta were associated with an escape from this immune response, and this effect was particularly sensitive for delta and omicron variant (p<10-6 for both). Interestingly, delta variant was associated with a 1.8-fold increased viral production rate (p = 0.046), while conversely omicron variant was associated with a 14-fold reduction in viral production rate (p<10-6). During a natural infection, our models predict that delta variant is associated with a higher peak viral RNA than omicron variant (7.6 log10 copies/mL 95% CI 6.8-8 for delta; 5.6 log10 copies/mL 95% CI 4.8-6.3 for omicron) while having similar peak infectious titers (3.7 log10 PFU/mL 95% CI 2.4-4.6 for delta; 2.8 log10 PFU/mL 95% CI 1.9-3.8 for omicron). These results provide a detailed picture of the effects of VoC on total and infectious viral load and may help understand some differences observed in the patterns of viral transmission of these viruses.

Nonlinear multilevel joint model for individual lesion kinetics and survival to characterize intra-individual heterogeneity in patients with advanced cancer.

In advanced cancer patients, tumor burden is calculated using the sum of the longest diameters (SLD) of the target lesions, a measure that lumps all lesions together and ignores intra-patient heterogeneity. Here, we used a rich dataset of 342 metastatic bladder cancer patients treated with a novel immunotherapy agent to develop a Bayesian multilevel joint model that can quantify heterogeneity in lesion dynamics and measure their impact on survival. Using a nonlinear model of tumor growth inhibition, we estimated that dynamics differed greatly among lesions, and inter-lesion variability accounted for 21% and 28% of the total variance in tumor shrinkage and treatment effect duration, respectively. Next, we investigated the impact of individual lesion dynamics on survival. Lesions located in the liver and in the bladder had twice as much impact on the instantaneous risk of death compared to those located in the lung or the lymph nodes. Finally, we evaluated the utility of individual lesion follow-up for dynamic predictions. Consistent with results at the population level, the individual lesion model outperformed a model relying only on SLD, especially at early landmark times and in patients with liver or bladder target lesions. Our results show that an individual lesion model can characterize the heterogeneity in tumor dynamics and its impact on survival in advanced cancer patients.

Evaluation of model-based bioequivalence approach for single sample pharmacokinetic studies.

In a traditional pharmacokinetic (PK) bioequivalence (BE) study, a two-way crossover study is conducted, PK parameters (namely the area under the time-concentration curve [AUC] and the maximal concentration [ C max ]) are obtained by noncompartmental analysis (NCA), and the BE analysis is performed using the two one-sided test (TOST) method. For ophthalmic drugs, however, only one sample of aqueous humor, in one eye, per eye can be obtained in each patient, which precludes the traditional BE analysis. To circumvent this issue, the U.S. Food and Drug Administration (FDA) has proposed an approach coupling NCA with either parametric or nonparametric bootstrap (NCA bootstrap). The model-based TOST (MB-TOST) has previously been proposed and evaluated successfully for various settings of sparse PK BE studies. In this paper, we evaluate, via simulations, MB-TOST in the specific setting of single sample PK BE study and compare its performance to NCA bootstrap. We performed BE study simulations using a published PK model and parameter values and evaluated multiple scenarios, including study design (parallel or crossover), sampling times (5 or 10 spread across the dosing interval), and geometric mean ratio (of 0.8, 0.9, 1, and 1.25). Using the simulated structural PK model, MB-TOST performed similarly to NCA bootstrap for AUC. For C max , the latter tended to be conservative and less powerful. Our research suggests that MB-TOST may be considered as an alternative BE approach for single sample PK studies, provided that the PK model is correctly specified and the test drug has the same structural model as the reference drug.

Pharmacokinetic Model of Tenofovir and Emtricitabine and Their Intracellular Metabolites in Patients in the ANRS 134-COPHAR 3 Trial Using Dose Records.

Tenofovir (TFV) and emtricitabine (FTC) are part of the recommended highly active antiretroviral therapy (ART). Both molecules show a large interindividual pharmacokinetic (PK) variability. Here, we modeled the concentrations of plasma TFV and FTC and their intracellular metabolites (TFV diphosphate [TFV-DP] and FTC triphosphate [FTC-TP]) collected after 4 and 24 weeks of treatment in 34 patients from the ANRS 134-COPHAR 3 trial. These patients received daily (QD) atazanavir (300 mg), ritonavir (100 mg), and a fixed-dose combination of coformulated TFV disoproxil fumarate (300 mg) and FTC (200 mg). Dosing history was collected using a medication event monitoring system. A three-compartment model with absorption delay (Tlag) was selected to describe the PK of, respectively, TFV/TFV-DP and FTC/FTC-TP. TFV and FTC apparent clearances, 114 L/h (relative standard error [RSE] = 8%) and 18.1 L/h (RSE = 5%), respectively, were found to decrease with age. However, no significant association was found with the polymorphisms ABCC2 rs717620, ABCC4 rs1751034, and ABCB1 rs1045642. The model allows prediction of TFV-DP and FTC-TP concentrations at steady state with alternative regimens.

Assessing the Increased Variability in Individual Lesion Kinetics During Immunotherapy: Does It Exist, and Does It Matter?

PURPOSE: Several studies have raised the hypothesis that immunotherapy may exacerbate the variability in individual lesions, increasing the risk of observing divergent kinetic profiles within the same patient. This questions the use of the sum of the longest diameter to follow the response to immunotherapy. Here, we aimed to study this hypothesis by developing a model that estimates the different sources of variability in lesion kinetics, and we used this model to evaluate the impact of this variability on survival. METHODS: We relied on a semimechanistic model to follow the nonlinear kinetics of lesions and their impact on the risk of death, adjusted on organ location. The model incorporated two levels of random effects to characterize both between- and within-patient variability in response to treatment. The model was estimated on 900 patients from a phase III randomized trial evaluating programmed death-ligand 1 checkpoint inhibitor atezolizumab versus chemotherapy in patients with second-line metastatic urothelial carcinoma (IMvigor211). RESULTS: The within-patient variability in the four parameters that characterize individual lesion kinetics represented between 12% and 78% of the total variability during chemotherapy. Similar results were obtained during atezolizumab, except for the durability of the treatment effects, for which the within-patient variability was markedly larger than during chemotherapy (40% v 12%, respectively). Accordingly, the occurrence of divergent profile consistently increased over time in patients treated with atezolizumab and was equal to about 20% after 1 year of treatment. Finally, we show that accounting for the within-patient variability provided a better prediction of most at-risk patients than a model relying solely on the sum of the longest diameter. CONCLUSION: Within-patient variability provides valuable information for the assessment of treatment efficacy and the detection of at-risk patients.

Guiding safer risperidone prescribing in Alzheimer's disease with therapeutic drug monitoring.

Previous analysis of pharmacokinetic data on risperidone-treated patients with dementia predicted that 20% had concentration-to-dose (C/D) ratios of the active moiety (risperidone and 9-hydroxy(OH)-risperidone) above 14 ng/mL per mg/day, which were in turn associated with a greater risk of extrapyramidal side effects. This study aimed to further explore risperidone pharmacokinetics in a second dataset. Nonlinear mixed effects modelling, using a Bayesian approach, was applied to data from a randomized controlled trial of risperidone in people with dementia. Covariates included age and glomerular filtration rate (GFR). Age had a significant effect on risperidone clearance (ß = -1.5) and GFR on 9-OH-risperidone clearance (ß = 0.2). The model predicted that 26.2% (95% confidence interval 18.6-32.6%) had C/D ratios above 14 ng/mL per mg/day. These findings confirm the importance of age-related risperidone dose adjustments and argue strongly for therapeutic drug monitoring in the initial stages of treatment to identify those at greatest risk of toxicity.

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.

Modelling the association between biomarkers and clinical outcome: An introduction to nonlinear joint models.

Nonlinear joint models are a powerful tool to precisely analyse the association between a nonlinear biomarker and a time-to-event process, such as death. Here, we review the main methodological techniques required to build these models and to make inferences and predictions. We describe the main clinical applications and discuss the future developments of such models.

Efficient model-based bioequivalence testing.

The classical approach to analyze pharmacokinetic (PK) data in bioequivalence studies aiming to compare two different formulations is to perform noncompartmental analysis (NCA) followed by two one-sided tests (TOST). In this regard, the PK parameters area under the curve (AUC) and $C_{\max}$ are obtained for both treatment groups and their geometric mean ratios are considered. According to current guidelines by the U.S. Food and Drug Administration and the European Medicines Agency, the formulations are declared to be sufficiently similar if the $90\%$ confidence interval for these ratios falls between $0.8$ and $1.25 $. As NCA is not a reliable approach in case of sparse designs, a model-based alternative has already been proposed for the estimation of $\rm AUC$ and $C_{\max}$ using nonlinear mixed effects models. Here we propose another, more powerful test than the TOST and demonstrate its superiority through a simulation study both for NCA and model-based approaches. For products with high variability on PK parameters, this method appears to have closer type I errors to the conventionally accepted significance level of $0.05$, suggesting its potential use in situations where conventional bioequivalence analysis is not applicable.

Quantifying the relationship between SARS-CoV-2 viral load and infectiousness.

The relationship between SARS-CoV-2 viral load and infectiousness is poorly known. Using data from a cohort of cases and high-risk contacts, we reconstructed viral load at the time of contact and inferred the probability of infection. The effect of viral load was larger in household contacts than in non-household contacts, with a transmission probability as large as 48% when the viral load was greater than 1010 copies per mL. The transmission probability peaked at symptom onset, with a mean probability of transmission of 29%, with large individual variations. The model also projects the effects of variants on disease transmission. Based on the current knowledge that viral load is increased by two- to eightfold with variants of concern and assuming no changes in the pattern of contacts across variants, the model predicts that larger viral load levels could lead to a relative increase in the probability of transmission of 24% to 58% in household contacts, and of 15% to 39% in non-household contacts.

Impact of study design and statistical model in pharmacogenetic studies with gene-treatment interaction.

Gene-treatment interactions, just like drug-drug interactions, can have dramatic effects on a patient response and therefore influence the clinician decision at the patient's bedside. Crossover designs, although they are known to decrease the number of subjects in drug-interaction studies, are seldom used in pharmacogenetic studies. We propose to evaluate, via realistic clinical trial simulations, to what extent crossover designs can help quantifying the gene-treatment interaction effect. We explored different scenarios of crossover and parallel design studies comparing two symptom-modifying treatments in a chronic and stable disease accounting for the impact of a one gene and one gene-treatment interaction. We varied the number of subjects, the between and within subject variabilities, the gene polymorphism frequency and the effect sizes of the treatment, gene, and gene-treatment interaction. Each simulated dataset was analyzed using three models: (i) estimating only the treatment effect, (ii) estimating the treatment and the gene effects, and (iii) estimating the treatment, the gene, and the gene-treatment interaction effects. We showed how ignoring the gene-treatment interaction results in the wrong treatment effect estimates. We also highlighted how crossover studies are more powerful to detect a treatment effect in the presence of a gene-treatment interaction and more often lead to correct treatment attribution.

Towards safer risperidone prescribing in Alzheimer's disease.

BACKGROUND: In the treatment of psychosis, agitation and aggression in Alzheimer's disease, guidelines emphasise the need to 'use the lowest possible dose' of antipsychotic drugs, but provide no information on optimal dosing. AIMS: This analysis investigated the pharmacokinetic profiles of risperidone and 9-hydroxy (OH)-risperidone, and how these related to treatment-emergent extrapyramidal side-effects (EPS), using data from The Clinical Antipsychotic Trials of Intervention Effectiveness in Alzheimer's Disease study (Clinicaltrials.gov identifier: NCT00015548). METHOD: A statistical model, which described the concentration-time course of risperidone and 9-OH-risperidone, was used to predict peak, trough and average concentrations of risperidone, 9-OH-risperidone and 'active moiety' (combined concentrations) (n = 108 participants). Logistic regression was used to investigate the associations of pharmacokinetic biomarkers with EPS. Model-based predictions were used to simulate the dose adjustments needed to avoid EPS. RESULTS: The model showed an age-related reduction in risperidone clearance (P < 0.0001), reduced renal elimination of 9-OH-risperidone (elimination half-life 27 h), and slower active moiety clearance in 22% of patients, (concentration-to-dose ratio: 20.2 (s.d. = 7.2) v. 7.6 (s.d. = 4.9) ng/mL per mg/day, Mann-Whitney U-test, P < 0.0001). Higher trough 9-OH-risperidone and active moiety concentrations (P < 0.0001) and lower Mini-Mental State Examination (MMSE) scores (P < 0.0001), were associated with EPS. Model-based predictions suggest the optimum dose ranged from 0.25 mg/day (85 years, MMSE of 5), to 1 mg/day (75 years, MMSE of 15), with alternate day dosing required for those with slower drug clearance. CONCLUSIONS: Our findings argue for age- and MMSE-related dose adjustments and suggest that a single measure of the concentration-to-dose ratio could be used to identify those with slower drug clearance.

Population Pharmacokinetic Model of Plasma and Cellular Mycophenolic Acid in Kidney Transplant Patients from the CIMTRE Study.

BACKGROUND AND OBJECTIVE: Mycophenolate mofetil is widely used in kidney transplant recipients. Mycophenolate mofetil is hydrolysed by blood esterases to mycophenolic acid (MPA), the active drug. Although MPA therapeutic drug monitoring has been recommended to optimise the treatment efficacy by the area under the plasma concentration vs time curve, little is known regarding MPA concentrations in peripheral blood mononuclear cells, where MPA inhibits inosine monophosphate dehydrogenase. This study aimed to build a pharmacokinetic model using a population approach to describe MPA total and unbound concentrations in plasma and into peripheral blood mononuclear cells in 78 adult kidney transplant recipients receiving mycophenolate mofetil therapy combined with tacrolimus and prednisone. METHODS: Total and unbound plasma concentrations and peripheral blood mononuclear cell concentrations were assayed. A three-compartment model, two for plasma MPA and one for peripheral blood mononuclear cell MPA, with a zero-order absorption and a first-order elimination was used to describe the data. RESULTS: Mycophenolic acid average concentrations in peripheral blood mononuclear cells were well above half-maximal effective concentration for inosine monophosphate dehydrogenase and no relationship was found with the occurrence of graft rejection. Three covariates affected unbound and intracellular MPA pharmacokinetics: creatinine clearance, which has an effect on unbound MPA clearance, human serum albumin, which influences fraction unbound MPA and the ABCB1 3435 C>T (rs1045642) genetic polymorphism, which has an effect on MPA efflux transport from peripheral blood mononuclear cells. CONCLUSION: This population pharmacokinetic model demonstrated the intracellular accumulation of MPA, the efflux of MPA out of the cells being dependent on P-glycoprotein transporters. Nevertheless, further studies are warranted to investigate the relevance of MPA concentrations in peripheral blood mononuclear cells to dosing regimen optimisation.

Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy.

Treatment evaluation in advanced cancer mainly relies on overall survival and tumor size dynamics. Both markers and their association can be simultaneously analyzed by using joint models, and these approaches are supported by many softwares or packages. However, these approaches are essentially limited to linear models for the longitudinal part, which limit their biological interpretation. More biological models of tumor dynamics can be obtained by using nonlinear models, but they are limited by the fact that parameter identifiability require rich dataset. In that context Bayesian approaches are particularly suited to incorporate the biological knowledge and increase the information available, but they are limited by the high computing cost of Monte-Carlo by Markov Chains algorithms. Here, we aimed to assess the performances of the Hamiltonian Monte-Carlo (HMC) algorithm implemented in Stan for inference in a nonlinear joint model. The method was validated on simulated data where HMC provided proper posterior distributions and credibility intervals in a reasonable computational time. Then the association between tumor size dynamics and survival was assessed in patients with advanced or metastatic bladder cancer treated with atezolizumab, an immunotherapy agent. HMC confirmed limited sensitivity to prior distributions. A cross-validation approach was developed and identified the current slope of tumor size dynamics as the most relevant driver of survival. In summary, HMC is an efficient approach to perform nonlinear joint models in a Bayesian framework, and opens the way for the use of nonlinear models to characterize both the rapid dynamics and the intersubject variability observed during cancer immunotherapy treatment.