This single-centre, retrospective pharmacokinetic study of VCZ was conducted at Union Hospital in Wuhan, China, from Feb 2017 to July 2018. The inclusion criteria were as follows: 1) age>12 years old; 2) patients treating with oral VCZ for IFIs (possible, probable or proven). The exclusion criteria included: 1) intolerant to VCZ treatment; 2) critical data were missing; 3) participated in another clinical trial simultaneously.

This study was approved by the Ethics Committee of Tongji Medical College, Huazhong University of Science and Technology (No: IORG0003571) and written informed consent was obtained from all patients.

VCZ was administered orally at a dose of 200 mg twice daily without a loading dose. Whole blood samples (5 ml) were collected for routine therapeutic drug monitoring of VCZ trough concentration and CYP2C19 genotyping. Most samples were taken when steady state was attained after 5 days of dosing, while others were collected and tested when the physicians deemed it necessary. The samples were drawn into EDTA-K₂-containing tubes before the following dose administration. The plasma was obtained after 15 min centrifugation under 1500 g and stored at –80°C until assay.

Demographic and physiological information of all patients were extracted from the electronic medical records database, including gender, age, height, body weight (WT), total bilirubin (TBIL), direct bilirubin (DBIL), indirect bilirubin (IBIL), total bile acids (TBA), alanine aminotransferase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP), gamma-glutamyl transferase (GGT), total protein (TP), albumin (ALB), globulin (GLB), blood urea nitrogen (BUN), uric acid (UA), serum creatinine concentration (SCR). The body surface area (BSA) was calculated based on the Mosteller formula. In addition, concomitant medications (proton-pump inhibitors and glucocorticoids) were recorded.

A liquid-liquid extraction (LLE) method using methyl tert-butyl ether (MTBE) as the extraction solvent was applied for sample preparation. Patient plasma (95 μL) was spiked with 5 μL of 80% MeOH followed by the addition of 10 μL internal standard (IS, propranolol). Then, 1,000 μL MTBE was added into the above mixture and vortexed for 3 min at 4°C. After centrifugation at 12000 g for 10 min, 800 μL supernatant was aspirated and dried in a clean tube under vacuum condition and then reconstituted in 80 μL of acetonitrile. Finally, 10 μL of the processed sample was injected into the liquid chromatography-tandem mass spectroscopy (LC-MS/MS) system for analysis.

Liquid chromatography was performed in a Shimadzu Prominence UFLC system (Shimadzu Corporation, Kyoto, Japan) with an Ultimate UHPLC XB-C18 column (50 mm × 2.1 mm, 1.8 μm, Welch, China) at 35°C. The mobile phase, pumped at a flow rate of 0.35 ml/min, consisted of 2 mM ammonium formate (mobile phase A) and acetonitrile (mobile phase B). The gradient elution conditions were as follows: 0–0.1 min, 10% B; 0.1–0.5 min, 10–50% B; 0.5–1.0 min, 50% B; 1.0–1.1 min, 50–60% B; 1.1–3.0 min, 60% B; 3.0–3.1 min, 60–10% B; 3.1–4.0 min, 10% B. An API-4000 Q Trap triple quadrupole mass spectrometer (AB Sciex, Foster City, CA, United States) was used for detection of the analytes. After optimization, tandem mass spectrometric detections were performed under the following operational parameters: curtain gas, 10 psi; collision-activated dissociation gas, 4; ion spray voltage, 5,500 V; source temperature, 500°C; gas 1, 50 psi; gas 2, 50 psi; interface heater, on. The quantification was accomplished by electrospray ionization in positive ion mode with multiple reaction monitoring (MRM). The MRM transitions were m/z 350.1→280.8, m/z 366.0→142.9 and m/z 260.1→116.0 for VCZ, VNO and IS, respectively. The lower limit of quantification was 0.5 ng/ml for both analytes. The calibration curves were linear over a range of 0.5–100 ng/ml for both analytes. The intra- and inter-day precision determined by coefficients of variation were less than 10% for both analytes. The accuracies for VCZ and VNO were 89.2–109.8% and 88.1–94.3%, respectively.

Total Genomic DNA was isolated from whole blood samples with QIAamp DNA blood kits (Qiagen, Hilden, Germany) according to the instruction supplied by the manufacturer. Genomic polymorphisms of CYP2C19*2 (c.681G > A, rs 4244285), CYP2C19*3 (c.636G > A, rs4986893), and CYP2C19*17 (c.-806 C > T, rs12248560) were identified by polymerase chain reaction-restriction fragment length polymorphism (PCR-RFLP) method on an ABI 3730XL Genetic Analyzer (Applied Biosystems, Foster City, CA, United States) as previously described (Mafuru et al., 2019). When CYP2C19*2 and *3 alleles were not detected, the allele was identified to be CYP2C19 *1 (wild type). All variant alleles were in accordance with Hardy-Weinberg equilibrium. The CYP2C19 phenotype was determined according to the guideline by the Clinical Pharmacogenetics Implementation Consortium (CPIC) (Moriyama et al., 2017). The UM and RM were defined as *17/*17 homozygote and *1/*17 heterozygote, respectively. And the NM was assigned to patients with deficient allele heterozygote (e.g., *1/*1). The IM was defined as a patient carrying one loss-of-function allele in combination with one normal function allele (e.g., *1/*2, *1/*3, *2/*17). Besides, the PM was defined as a patient with deficient allele compound heterozygote or homozygote (e.g.,*2/*2, *2/*3 or *3/*3).

A PopPK analysis was conducted using the software Phoenix^® NLME (Version 8.2.0.4383, Pharsight Corporation, United States) to estimate the pharmacokinetic parameters of VCZ. The first-order conditional estimation-extended least squares (FOCE ELS) algorithm was applied throughout the modeling procedure. The R program (version 4.0.2, http://www.r-project.org/) was utilized for data visualization and model validation.

All VCZ concentrations in the present study were trough concentration data, which could not be used for characterization of the absorption phase. Thus, according to data from published literature, the absorption rate constant (k_a) and the bioavailability (F) were fixed at 1.1/h and 0.895, respectively (Pascual et al., 2012; Wang et al., 2014). Considering the limited sampling times, one-compartment models based on first-order absorption with either linear or non-linear (Michaelis–Menten) elimination were tested to fit the concentration profile of VCZ. In addition, compartmental and residual error models for VNO were identified to describe the metabolic pathway from VCZ to its N-oxide metabolite.

The inter-individual variability in pharmacokinetic parameters was estimated using an exponential function as follows: P i = θ×exp ( η i ) (1) Where P _i represents the estimated pharmacokinetic parameter for individual i, θ is the population typical value of that parameter, and η _i represents the random variable for individual i, which is defined as normally distributed with a mean of 0 and a variance of ω ².

The residual variability was evaluated by additive, proportional or combined additive-proportional residual error models. The equations were as follows: C i j = C p r e d , i j + ε i j (2) C i j = C p r e d , i j × ( 1 + ε i j ) (3) C i j = C p r e d , i j × ( 1 + ε i j 1 ) + ε i j 2 (4) Where C _ij is the jth observed value in individual I; C _pred,ij is the jth predicted value in individual i; ε _ij is the residual random error, which is assumed to be Gaussian distributed with a mean of 0 and a variance of σ ².The selection of the appropriate structural model was based on the visual inspection of diagnostic plots and improvements of the statistic parameters, including Akaike information criterion (AIC), Bayesian information criterion (BIC) and the objective function value (OFV).

After development of the base model, the potential covariates were tested using stepwise forward selection followed by backward elimination steps. Pairwise correlations between all variables were assessed prior to covariate analysis, and highly collinear variables (correlation coefficient >0.5) would not be simultaneously incorporated into the final model. Covariates associated with a significant decrease of OFV more than 3.84 units (Chi-squared distribution, df = 1, p ≤ 0.05) were added to the base model to establish a full model. Then the covariates were removed from the full model one by one. Covariates resulting in a significant increase of OFV by at least 6.63 units (Chi-squared distribution, df = 1, p ≤ 0.01) were retained in the final model. The possible influences of all demographic and physiological variables were explored. In the screening process, power models (Eq. (5)) were used for continuous covariates such as age, body weight, BSA, renal and liver function parameters. While exponential models (Eq. (6)) were used for categorical covariates such as sex, concomitant medications, andCYP2C19 metabolic phenotypes. θ i = θ × ( C o v j C o v m e d i a n ) θ c o v (5) θ i = θ × exp ( θ cov ) (6) where Cov_j is the jth covariate, Cov_median is the median value of covariate, θ_i is the population prediction of the pharmacokinetic parameter, θ is the population typical value of the parameter, and θ_covdescribes the fixed effect of the covariate on the parameter.

A post hoc empirical Bayesian method was employed to estimate individual exposure pharmacokinetic parameters at steady state for both VCZ and VNO, including trough concentration (C_min), peak concentration (C_max) and the area under drug plasma concentration-time curve over 24 h at day 7 (AUC, 144–168 h) with twice-daily oral administration of 200 mg VCZ. Additionally, metabolic ratio (MR) was calculated using Equation (7) (Yamada et al., 2015). SPSS software version 19.0 (SPSS Inc. Chicago, IL, United States) was used for exploratory data analysis. MR = AUC VNO AUC VCZ (7) where AUC_VNO and AUC_VCZ are the area under drug plasma concentration-time curve of VNO and VCZ, respectively.

The goodness of fit plots, nonparametric bootstrap method, and visual predictive check (VPC) were employed for model validation. Different diagnostic plots such as observed concentrations (DV) vs individual predictions (IPRED), DV vs population predictions (PRED), Conditional weighted residuals (CWRES) vs time, and CWRES vs PRED graphs were drawn to visually assess the accuracy of the PopPK model. A nonparametric bootstrap procedure was carried out using 1,000 randomly resampled datasets generated from the original data. To evaluate model stability, the 2.5th, 50th, and 97.5th percentile of the bootstrap estimates were computed and compared to the final model parameter estimates. The VPC method was performed for both VCZ and VNO by simulating 1,000 virtual subjects based on the final model estimates. The simulated concentration profiles were graphically compared with observed concentrations to evaluate the predictive performance of the final model.

The Monte Carlo simulations were performed based on the parameter estimates derived from the final model. Each simulation runs for 1,000 times to predict the VCZ concentration profiles following multiple oral doses. To evaluate and optimize the currently used dosing regimen, oral maintenance doses of 150–400 mg twice daily (BID) or three times daily (TID) were simulated in patients with various CYP2C19 phenotypes. Different concentration cut-off values for efficacy and toxicity (1, 2 and 5.5 mg/L) were employed to determine the probability of supratherapeutic or subtherapeutic levels with each dosing scenario.

A total of 78 patients participated in this study and no patient was excluded from PopPK analysis. These patients ranged in age from 14 to 70 years and in weight from 44 to 111 kg. Of these patients, only four were adolescents between the ages of 14 and 17, while the rest were adults. All participants provided 214 plasma concentrations for VCZ (range 0.01–7.34 μg/ml) and 213 plasma concentrations for VNO (range 0.04–7.89 μg/ml) measured, with a maximum of eight samples per patient. The sampling time after the first dose ranged from 23 to 4,223 h. The plasma concentration versus time profiles for both VCZ and VNO are shown in Figure 1. The genotyping results were obtained from 75 patients, the majority (n = 32) of them were categorized as IMs, 27 as NMs, and 16 as PMs. The CYP2C19*17 allele was not detected in the genotype analysis, which may be due to the limited frequency of the *17allele in Chinese population (Sim et al., 2006) or the sample size of our study. The remaining three patients were excluded from subgroup comparison analysis due to lack of genotyping information. The demographic and clinical characteristics of all patients are summarized in Table 1.

The concentration-time profile of VCZ was well fitted to a one-compartment model with first-order absorption and mixed linear and concentration-dependent nonlinear (Michaelis-Menten) elimination. Meanwhile, the disposition of VNO was adequately described by a one-compartment model, with an input rate the same as the conversion rate from VCZ to VNO (Figure 2). For both VCZ and VNO, the inter-individual variability and the residual variability could be best expressed using an exponential model and a proportional model, respectively. In this joint pharmacokinetic model, VCZ was represented by a central compartment (A₁), parameterized by the nonlinear clearance of VCZ to the metabolite VNO (CL_nonlin), the clearance of VCZ other than the metabolic pathway converting to VNO (CL₁), and the volume of distribution (V₁). VNO was described by a single metabolite compartment (A₂), in which the clearance (CL₂) and the distribution volume (V₂) of VNO were used. Hence, the equations that illustrated the final structural model were as follows: d A g u t d t = − k a × A g u t (8) d A 1 d t = F × A g u t × k a − C L 1 V 1 × A 1 − C L n o n l i n V 1 × A 1 (9) d A 2 d t = C L n o n l i n V 1 × A 1 × k n − C L 2 V 2 × A 2 (10) C L n o n l i n e = V max C 1 + k m × ( 1 − I m a x ⋅ C 2 I C 50 + C 2 ) (11) C 1 = A 1 V 1 (12) C 2 = A 2 V 2 (13) Where A_gut is the amount of VCZ in the absorption site, k_a is the absorption rate constant, F is the oral bioavailability, V_max is the maximum elimination rate, K_m is the Michaelis-Menten constant, C₁ and C₂ represent the plasma concentrations of VCZ and VNO, respectively. I_max is the maximal inhibitory effect, IC₅₀ is the concentration of VNO yielding 50% of maximum clearance inhibition. In this model, k_m, I_max, and IC₅₀ were fixed at 1.15 mg/L, 0.75 and 14.6 mg/L, respectively (Liu et al., 2014; Hohmann et al., 2017).

For basic model, the OFV, AIC and BIC were 1,320.48, 1,342.48 and 1,387.11, respectively. Covariate screening analysis indicated that the age, gender, WT and BSA of the patient as well as liver function parameters did not have any impact on the pharmacokinetic parameters to a statistically significant extent in the study population.Only CYP2C19 phenotype was considered to have a significant influence on V_max.The final model with this covariate decreased the OFV, AIC and BIC by 17.55, 13.55 and 5.44 units, respectively.When the base model was updated with the final model, the number of parameters was increased from 11 to 13. The detailed PopPK parameter estimates derived from the final model are given in Table 2.

The exposure parameters of VCZ and VNO estimated from the final joint model are presented in Supplementary Table S1. The one-way ANOVA and LSD test were used to compare the pharmacokinetic parameters in three groups classified according to CYP2C19 phenotypes. As shown in Figure 3, the result revealed remarkable differences in VCZ exposure (C_min-VCZ, C_max-VCZ, and AUC_VCZ) among the three groups (p < 0.05).

The goodness-of-fit plots presented in Figure 4 revealed acceptable agreement between the predicted and observed concentrations. The CWRES were generally distributed around zero with no trend, and the majority of them were in the -3 to +3 range. The results of bootstrap analysis are summarized in Table 2. All parameter estimates obtained from the final PopPK model laid within the 95% confidence intervals (CIs) calculated using the bootstrap method and were close to the median values with small bias (<10%), suggesting the robust stability of the final model. The VPC plots depicted in Figure 5 showed that the model performed well in predicting the plasma concentrations of both VCZ and VNO, and no evidence of model misspecification was found.

The predicted VCZ concentration-time profiles during 20 days after administration of standard doses (400 mg BID for two doses followed by 200 mg BID) were simulated according to the CYP2C19 phenotypes. As shown in Figure 6, the trough concentrations of VCZ at steady state attained the therapeutic range of 2–5.5 mg/L in PMs, whereas the trough levels in the subjects with IM and NM genotypes were below the therapeutic range. The assessment of the probability of achieving target through level of >2 mg/L with 200 mg twice daily oral dosing regimen indicated that up to 70.31% of NMs, 54.48% of IMs and 41.71% of PMs failed to reach the lower end of the therapeutic range. These results suggested that dosage adjustment based on CYP2C19 phenotype was needed to attain adequate exposure and to improve clinical response.

Table 3 shows the probability of VCZ target C_min attainment for different dosage regimens stratified by CYP2C19 phenotype. Considering the trade-offs between toxicity and efficacy, we concluded that the following regimens were appropriate for patients with different CYP2C19 phenotypes: 325 mg bid or 200 mg tid for NM patients, 275 mg bid or 175 mg tid for IM patients, and 225 mg bid or 150 mg tid for PM patients. The results also demonstrated that a lower C_min target of >1 mg/L would lead to higher PTA.