This was a single-center, prospective, open-label, observational study. The entire study design is shown in Figure 1. Patients were aged ≥18 years and underwent their first heart transplantation surgery at Guangdong Provincial People’s Hospital and received triple maintenance immunosuppressive therapy comprising MMF, TAC, and corticosteroids. Patients undergoing combined organ transplantation were excluded. The study protocols were approved by the Research Ethics Committee of Guangdong Provincial People’s Hospital [2018478H (R1)], and all participants provided written informed consent before inclusion. The study protocol was registered in the Chinese Clinical Trial Registry (No. ChiCTR2000030903). Volunteers were the sole source of transplants in the current study. There were no executed donors in the current study. Generally, in China, doctors select heart donors under the age of 50. The mean age of the donors in the present study was 35.1 ± 9.4 years, and six donors were females. The donor cause of death was trauma, cerebrovascular accident, brain tumor, or hypoxic brain injury.

MMF dispersible tablets (Cycopin®, manufactured by Huadong Medicine Co., Ltd.) were orally administered at a maintenance dose of 250–750 mg twice daily (bid). TAC (Prograf®) was dosed at 0.5–4 mg once or twice a day after transplantation and therapeutic drug monitoring was conducted to achieve the target concentration. Methylprednisolone (20 mg) was administered during the operation and was switched to 8 mg in the following days. Sample collection for the measurement of MPA concentrations started at least 7 days post-operation. Full pharmacokinetic profiles were obtained by intensive sampling, drawing 10 blood samples: pre-dose, 0.5, 1, 1.5, 2, 3, 4, 6, 8, and 12 h after the morning dose. For sparse sample collection, blood samples were drawn at 0.5, 1.5, 4, and 9 h post-dose after the morning dose. More than one pharmacokinetic sampling cycle might be collected in one patient.

Patient information was collected during sample collection from the medical electronic records, including demographic information, biochemical indices, concomitant medications, TAC concentrations and daily dose, and postoperative time (POT). Proton pump inhibitors (PPIs) included omeprazole injections or gastro-resistant tablets and pantoprazole injections or enteric-coated tablets.

Blood samples were collected in EDTA-coated tubes and centrifuged at room temperature (2,000 × g for 10 min) to separate plasma, which was analyzed within 24 h after collection. MPA concentrations were detected by the EMIT using MPA kits on the Viva-E system (Siemens Healthcare Diagnostics, Newark, Del, United States) according to the manufacturer’s guidelines. The linear range of MPA concentrations was 0.1–15 mg/L. The lower limit of quantification (LLOQ) was 0.1 mg/L. At the control levels of 1, 7.5, and 12 mg/L, the precisions were 5–6%, 3.1–4.5%, and 4.3–6.5%, respectively.

The intensive and sparse concentration data were pooled and simultaneously fitted using a nonlinear mixed-effects modeling software program (NONMEM, version 7.3.0, ICON Development Solutions, MD, United States). The units of MMF doses and MPA concentrations were unified as moles by multiplying the molecular weights. The first-order conditional estimation with interaction (FOCEI) method was used throughout the entire model-building process. One- or two-compartment models with time-lagged first-order absorption and first-order elimination or a two-compartment model with EHC were tried as the structure model. Oral bioavailability (F) was fixed as 0.95 as in Bullingham et al. (1996) and Armstrong et al. (2005). Inter-individual variability (IIV) of each PK parameter was modeled using the following exponential error model: P i = P t v · e η i , (1) where P _i represents the PK parameter of the ith individual, P _tv represents the typical population value, and η _i is the inter-individual random effect with mean 0 and variance ω². The covariance between IIV values was estimated using a variance–covariance matrix.

The influence of covariates, including sex, age, body weight (WT), height (HT), body mass index (BMI), hemoglobin (HGB), hematocrit (HCT), serum creatinine (Scr), serum albumin (ALB), total protein (TP), blood urea nitrogen (BUN), uric acid (URIC), total bilirubin (TBIL), direct bilirubin (DBIL), diuretic co-medication (1 for yes/0 for no), PPI co-medication (1 for yes/0 for no), and TAC daily dose and trough concentrations on the main pharmacokinetic parameters was evaluated. The estimated glomerular filtration rate (eGFR) was calculated by the Modification of Diet in Renal Disease (MDRD) equation (modified for Chinese): eGFR = 175 × (Scr)^−1.234 × (Age)^−0.179 × (0.79 if female) (Ma et al., 2006). Discrete covariates (such as sex and co-medication of PPIs) were modeled as follows: P i = P t v · θ C O V i · e η i . (2) Here, COV _i represents the discrete covariate value of the ith individual as a binary variable (1 or 0) and Θ represents the factor of the covariate.

For continuous covariates (such as WT and ALB), an exponential or linear model with average covariate values and an adjusting factor were tried as follows:     P i = P t v · ( C O V i C O V a v e ) θ · e η i , (3) P i = P t v · ( 1 + ( C O V i − C O V a v e ) · θ ) · e η i . (4) Here, COV _i represents the continuous covariate value of the ith individual, COV _ave represents the average values of the covariates, and Θ represents the factor of the covariate.

The objective function value (OFV) and goodness-of-fit (GOF) plots were used as the criteria for model selection. For forward inclusion, it was considered significant if the OFV decreased more than 6.63 (chi-squared distribution, d = 1, p < 0.01). Backward exclusion had a stricter criterion of increase in the OFV, that is, OFV >10.83 (chi-squared distribution, d = 1, p < 0.001).

The adequacy of the final model was assessed using GOF plots, the normalized prediction distribution error (NPDE) and prediction-corrected visual predictive check (pcVPC). 1,000 bootstraps for the 95% confidence intervals (CIs) of the parameters (n = 1,000 times) were conducted in Perl-speaks-NONMEM (PsN) version 4.9.0 (psn.sourceforge.net, Uppsala University, Sweden) (Lindbom and Jonsson, 2004). The GOF plots were generated using R software (version 3.6.3) (R Core Team, 2019). The NPDE analysis and plotting were implemented using the NPDE add-on package (version 2.0) (Scrucca et al., 2016). PcVPC was performed by simulations (n = 1,000) based on the final PPK model using PsN and plotted by R.

The predictive performance of AUC_0–12h by the final PPK model was evaluated using intensive data from 14 patients. The AUC_0–12h value of MPA based on the intensive concentrations versus the time profile (marked as AUC_obs) was calculated directly by the linear trapezoidal rule using the R add-on package PKNCA (version 0.9.2) (Denney et al., 2015). The individual-predicted AUC_0–12h (AUC_ipred) was calculated by the simulated pharmacokinetic profiles with a virtual sampling increment of 0.5 h, using the individual PK parameters estimated by the final PPK model. The bias was expressed as the mean error in percentage (%ME), and precision was evaluated using percent root mean squared error (%RMSE) as in the following equations. % M E = 1 N ∑ i = 1 i = N A U C p r e d , i − A U C r e f , i A U C r e f , i · 100 , (5) % R M S E = 1 N ∑ i = 1 i = N ( A U C p r e d ,   i − A U C r e f , i A U C r e f , i ) 2 · 100. (6) Here, AUC _{ref, i} is the AUC_obs of the i ^th patient and AUC _{pred, i} is the AUC_ipred of the i ^th patient. Passing–Bablok regression analysis and Bland–Alman plots were used to estimate the relationship between predicted and observed AUC_0–12h values using an R add-on package (MethComp 1.30.0) (Carstensen et al., 2020).

Simulations were performed on the final PPK model to investigate the impact of the significant covariates on MPA exposure. Full pharmacokinetic profiles at steady-state, with time point increments of 0.5 h during a 12-h dosing interval, were simulated at a virtual dose of 500 mg MMF in 24 scenarios, including eGFR values of 30, 60, 90, and 130 ml/min/1.73 m², ALB values of 30, 40, and 60 g/L, and with or without PPIs. The eGFR values were used to simulate patients with severe or mild renal impairment, normal renal function, and augmented renal clearance (ARC). The selected ALB value of 30 g/L was used to simulate patients with hypoproteinemia. The AUC_0–12h values were calculated by the linear trapezoidal rule based on population-predicted (PRED) values using R.

For the limited sampling strategy (LLS), we used two methods: 1) optimal sampling times for Bayesian estimation by ED-optimal sampling and 2) a fixed time-point linear equation using MLR.

For the accurate individual pharmacokinetic estimates, an ED-optimal sampling strategy was applied to optimize the sampling time points using the R package PopED (version 0.6.0) (Nyberg et al., 2012). For the current study, the sampling timepoints were constrained to four within 6 h after drug intake, including pre-dose (0 h). To evaluate the accuracy of AUC_0–12h estimation based on the optimal sampling times, we performed simulations to generate 100 virtual patients receiving MMF at an oral dose of 500 mg. The covariate values (i.e., eGFR, ALB, and co-medication of PPIs) were generated using the mean and standard deviation of the original patients. Intensive concentrations of the virtual patients during a dosing interval were generated (as a reference dataset) based on the final PPK model. For the optimal sampling dataset, we assumed that the plasma concentration samples were only collected at the optimal sampling times from the reference dataset. AUC_0–12h was calculated with the linear trapezoidal rule using R. AUC_ref was calculated based on the reference dataset. AUC_opt was calculated based on the simulated concentrations by the individual PK parameters estimated by maximum a posteriori Bayesian (MAPB) estimation. We also selected the observed concentrations at the optimal sampling timepointss from the intensive dataset (n = 14) and estimated the individual PK parameters by MAPB. The AUC_0–12h value was marked as AUC_opt14. Predictive performance was evaluated by comparing the individual reference values with those estimated based on the optimal sampling dataset, respectively, and the %ME and %RMSE calculated by Eqs 5, 6. Here, AUC _{ref, i} is AUC_ref of the ith virtual patient or the AUC_obs of the ith intensive sampling patient and AUC _{pred, i} is AUC_opt of the ith virtual patient or the AUC_opt14 of the ith intensive sampling patient. Passing–Bablok regression analysis and Bland–Alman plots were used to estimate the correlation using R.

For clinical convenience, a fixed time-point linear equation was generated using simulated data based on the original dataset. Using the actual doses and the patients’ pharmacokinetic parameters obtained from the final PPK model, we simulated the complete pharmacokinetic profiles with a time increment of 0.5 h during the dosing interval. The simulated AUC_0–12h values were calculated by the linear trapezoidal rule as the reference (marked as AUC_ipred). With simulated data of each patient, we searched for predictive models of the AUC_ipred using multiple stepwise linear regression analysis with a limited sampling window ≤6 h post-dose and a maximum number of four samples. Data were analyzed using the software package SPSS 25.0 for Windows (SPSS Inc., Chicago, IL, United States). We considered linear equations with high coefficients (r ²) and clinical convenience as acceptable. The variance inflation factor (VIF) was calculated to check for collinearity of the variables, and VIF <10 was preferred. AUC_0–12h calculated by the final linear model was marked as AUC_MLR. The predictive accuracy was evaluated by comparing the AUC_ipred to AUC_MLR and the AUC_obs to AUC_MLR, respectively. The %ME and %RMSE were calculated by Eqs 5, 6. In the equations, AUC _{ref, i} is AUC_ipred of the ith patient and AUC _{pred, i} is AUC_MLR of the ith patient. To compare with the observed AUC_0–12h, AUC _{ref, i} is found to be AUC_obs of the ith patient from 14 intensive sampling patients and AUC _{pred, i} is AUC_MLR of the ith patient. Passing–Bablok regression analysis was performed, and Bland–Alman plots were generated to evaluate the accuracy of the linear equation.

In addition, optimal sampling based on the time-points in the final MLR equation was also evaluated, for there might be missing samples in clinical practice. If it happened, AUC_0–12h could also be estimated by MAPB estimation based on the linear sampling times.

In total, 508 plasma samples with 507 above LLOQ were obtained from 91 adult heart transplant recipients with 105 pharmacokinetic sampling cycles. Seventy-seven patients had one PK profile, 13 patients had two, and one patient had three PK profiles. The sampling time after operation (POT) was extensive, from 7 days to nearly 3 years. The intensive sampling was collected from 14 patients, and the other data were sparse. The patient characteristics are shown in Table 1. The recipients (seven females) had a median age of 50 (range, 21–74) years and a median bodyweight of 60.0 kg (range, 33.4–95.0 kg). Furthermore, the doses of MMF were in a range of 250–750 mg, co-medicated with TAC dosed from 1 to 8 mg per day.

MPA plasma concentrations were adequately fitted using the two-compartment model with lag time (T_Lag) and three covariates in the final model. The combined residual error model resulted in a proportional error of 26.1% and an additive residual error of 0.144 mg/L. The PPK parameter estimates with 95% CIs are shown in Table 2. The impacts of co-medication of PPIs on the bioavailability (F) of eGFR on the clearance of the central compartment (CL/F) and ALB on the volume of the central compartment (V₂/F) were significant covariates incorporated in the final PPK model. On the use of PPIs, the F of MPA decreased to 72.4%. The typical value of CL/F was 7.36 L/h. The factor of the eGFR covariate model (θ in Eq. 4) was 0.00791. The inclusion of eGFR explained 11.6% of the IIV of CL/F. With an incremental decrease of 30 ml/min/1.73 m² of eGFR, CL/F decreased by 23.7%. The impact factor of ALB on V₂/F was −7.31 using the exponential covariate model (Eq. 3), which resulted in a 30.8% reduction of IIV of V₂/F. The Q/F was estimated as 17.0 L/h, and the V₃/F was as large as 560 L with very high IIV (189.5%). The model-building process is shown in Supplementary Table S1.

The GOF plots showed an adequate fitness of the final model (Figure 2). There was a good agreement between population-predicted concentrations (PREDs) and the observed concentrations, except for underestimating the higher observations (Figure 2A). It improved in the individual-predicted concentrations (IPREDs) vs. the observed concentrations. However, a similar underestimation also occurred in the higher concentration range (>15 mg/L). With an acceptable range of the conditional weighted residuals (CWRES), CWRES indicated no obvious trends throughout the time and the predicted concentrations (Figures 2C,D). The NPDE plots are shown in Figure 3. The quantile–quantile (QQ) plot and the distribution histogram of the NPDE showed a mean of 0.0613 and a variance of 0.9812 (p = 0.164) (Figures 3A,B). There was no trend for the NPDE versus time (Figure 3C) or the NPDE versus the predicted concentrations (Figure 3D). These results indicated that the final PPK model of MPA was accurate and reliable. The pcVPC plot is shown in Supplementary Figure S1. The pcVPC plot shows that most of the observed concentrations were between the 5th and 95th percentiles of model predictions, except that the predictions were higher than the observations around the peak timepoints (1–2 h post-dose). Generally, the model adequately characterized MPA concentrations.

Compared with the observed AUC_0–12h values (AUC_obs) of the intensive data, the accuracy (%ME) and precision (%RMSE) of the AUC_0–12h estimation were 5.0 and 11.8%, respectively. High correlations between the AUC_obs and AUC_ipred (r = 0.977) and most of the differences within 1.96 SD limits are shown in Supplementary Figure S2. It demonstrated that the final PPK model predicted the AUC_0–12h values well.

The simulated AUC_0–12h values at the steady-state calculated by the PRED are listed in Table 3. The change of ALB values resulted in few changes in AUC_0–12h. However, PPI co-medication and eGFR change remarkably affected MPA exposure. With PPIs, AUC_0–12h was reduced by approximately 26%, which was uniform with the influence on F. The simulated AUC_0–12h values under severe renal impairment, mild renal impairment, and ARC were about 1.60-, 1.23-, and 0.78-fold that for normal renal function, respectively. Without PPIs, severe renal-impaired patients acquired the highest MPA exposure.

The ED-optimal sampling timepoints included three timepoints: pre-dose and 1 and 4 h post-dose. We found that, except for 1 or 4 h post-dose alone, taking any of the optimal sampling times could obtain accurate AUC_0–12h by MPAB estimation, with the value of correlation coefficients above 0.97 (Supplementary Table S2). We got the highest r ² (0.9981) with the %ME value of −0.1% and the %RMSE value of 2.9% taking all three sampling timepoints. It showed a high linear correlation between AUC_ref and AUC_opt (r = 0.998) (Figure 4A). The Bland–Altman plot showed that only eight AUC_opt values were outside the 1.96SD limit (Figure 4B). The linear coefficient between AUC_obs and AUC_opt14 was 0.843 (Figure 4C), with an acceptable bias (%ME = –3.3%), but with a larger precision (%RMSE = 32.4%). There was only one outlier in the Bland–Altman plot (Figure 4D).

The MLR analysis showed several one-to four-timepoint linear models with R ² > 0.9 (Supplementary Table S3). Considering estimation reliability, collinearity, and convenience of clinical application, a best-fit linear equation with four timepoints was finally chosen: AUC_0–12h = 3.539 × C₀ + 0.288 × C_0.5 + 1.349 × C₁ + 6.773 × C_4.5 (r ² = 0.999). No intercept was included in the linear equation because it was not significant (p = 0.471). High correlations between the PPK-predicted, linear-calculated, and observed AUC_0–12h values indicated the accuracy of the linear equation, with slopes of 1.0 (r = 1.00) and 0.95 (r = 0.976), respectively (Figures 5A,C). The %ME and %RMSE values between AUC_ipred and AUC_MLR were 0.33 and 1.75%, respectively. Compared with the AUC_obs, the %ME and %RMSE values were 6.1 and 11.0%, respectively. The Bland–Altman plot showed that only six AUC_MLR values (5.7%) were outside the 1.96SD limits with AUC_ipred as the reference (Figure 5B). With AUC_obs as a reference, there was only one outlier calculated by the linear equation (Figure 5D). The linear sampling times also showed adequate accuracy of AUC_0–12h estimation by MAPB estimation (Supplementary Table S2).