Sample assay
An improved HPLC-MS/MS method based previous literature [15] was used to measure plasma concentrations of PCr and Cr, using salidroside and creatine-(methly-d3) (Cr-d3) as internal standards, respectively. 50 μL of plasma and 10 μL of internal standard solution (20 μg/mL salidroside and 100 μg/mL Cr-d3 for PCr and Cr determination, respectively) were prepared by protein precipitation with acetonitrile/water (1000 μL, 1:1, v/v). After vortex and centrifugation, 100 μL of the supernatant fluid was diluted with 500 μL of solvent, and 10 μL was injected into the HPLC-MS/MS system for the determination of PCr. While, for the Cr determination, 20 μL of the supernatant fluid was diluted with 1000 μL of solvent, and 5 μL was injected into the HPLC-MS/MS system.
A Hypersil Gold C18 column (150×2.1 mm, 5 μm, Thermo Scientific, USA) was used for separation of PCr. The mobile phase was solution A: 2 mM ammonium acetate in water (pH 10, adjusted with ammonia), and solution B: methanol. Chromatography separation was obtained with gradient elution solvent, the gradient elution program was: 0-1.5 min, 98-20% A; 1.5-4 min, 20% A; 4-4.01 min, 20-98% A; 4.01-10 min, 98% A. Monitored ion pairs were m/z 210 → 79 and m/z 299 → 119 for PCr and salidroside in negative ionization mode.
For the separation of Cr, a Hypersil Gold C18 column (150×4.6 mm, 5 μm, Thermo Scientific, USA) was applied. The mobile phase was solution A: 2 mM ammonium acetate in water (0.4‰ formic acid), and solution B: acetonitrile. Isocratic elution with 30% A was obtained. Monitored ion pairs were m/z 132 → 90 and m/z 135 → 93 for Cr and Cr-d3 in positive ionization mode.
This method was accurate with the intra-day and inter-day precision less than 6%. The lower limit of quantification (LLOQ) for PCr and Cr were 0.5 μg/mL and 4 μg/mL, respectively. Measured concentrations below LLOQ were reported in the dataset.
Population pharmacokinetic model development
A sequential two-step analysis approach to model building was implemented [16]. First, a population pharmacokinetic model of PCr was developed and then parent parameters were fixed to develop the population pharmacokinetic model for Cr. The nonlinear mixed-effects modeling method was applied to establish the population pharmacokinetic model. The PCr and Cr concentrations were fitted using the Phoenix NLME software (version 8.2, Certara, St. Louis, MO) with the first-order conditional estimation-extended least squares (FOCE-ELS) method throughout the model building process. Pharmacokinetic data which were below LLOQ were handled using the M3 method [17]. Model selection criteria were based on goodness-of-fit plots, objective function value (OFV, equal to -2 log likelihood), akaike information criteria (AIC) and precision of parameter estimates. The complex parent-metabolite joint population pharmacokinetic model was established in a stepwise fashion. First, the base model including PCr and Cr was developed with residual error model and interindividual variability model. After that, covariate effects on base model were investigated to construct the covariate model.
For PCr, structural model was tested either one-compartment or two-compartment PK models. The one-compartment model parameters included volume of distribution (V) and central clearance (CL). For the two-compartment structural model, inter-compartmental clearance (Q) and volume of the second compartment (V2) were included. The structural model for Cr was also tested. The fraction (Fm) of PCr metabolized to Cr was fixed so as to avoid identifiability problem and the volumes of the compartment for Cr was estimated. A published literature reported that a large amount of PCr was metabolized as Cr in animals after intravenous injection of exogenous PCr, and the conversion rate was about three-quarters [10]. In this study, according to the references and the characteristics of drug concentration-time curves of PCr and Cr, the Fm of PCr converted to Cr was assumed to be 0.75. Because of the existence of endogenous Cr which was detected at each of the pre-dose sample, a parameter of the baseline of endogenous Cr levels was added. Allometric scaling based on bodyweight (BW) was applied to the PK parameters. An allometric power model was used with power exponents fixed at 0.75 for clearances and 1.0 for volumes of distribution, as described in the following equations [18, 19]:
\(\text{CL}_{i}\ =\ \theta_{\text{CL}}\ \times\left(\text{BW}_{i}/20\right)^{0.75}\)(1)
\(V_{i}\ =\ \theta_{V}\ \times\left(\text{BW}_{i}/20\right)^{1}\)(2)
In this expression, CLi and Vi are the typical clearance and central volume of distribution for an individual i with bodyweight BWi while θCL and θV are the respective parameter values for a subject with a bodyweight of 20 kilogram.
The inter-individual variability (η) of pharmacokinetic parameters was assumed to follow a log-normal distribution with a mean of 0 and a variance of ω2. Exponential formula was applied to account for IIV (Equation 3).
\({P_{\text{ij}}\ =\ P_{j}\ \times e}^{\eta_{\text{ij}}}\) (3)
where Pj and Pij represent the typical value of j th population parameter and i th individual’j th parameter.
Proportional error model (Equation 4), additive error model (Equation 5), additive and proportional error model (Equation 6), and power error model (Equation 7) were evaluated to describe the residual error.
\(C_{\text{ij}}\ =\ \text{IPRED}_{\text{ij}}\ \times(1+\ \varepsilon_{\text{ij}})\)(4)
\(C_{\text{ij}}\ =\ \text{IPRED}_{\text{ij}}\ +\ \varepsilon_{\text{ij}}\)(5)
\(C_{\text{ij}}\ =\ \text{IPRED}_{\text{ij}}\ \times\left(1+\ \varepsilon_{ij,1}\right)+\ \varepsilon_{ij,2}\)(6)
\(C_{\text{ij}}\ =\ \text{IPRED}_{\text{ij}}\ +\ {{\text{IPRED}_{\text{ij}}}^{0.5}\times\varepsilon}_{\text{ij}}\)(7)
where Cij was the observation concentration ofi th individual at the j th sampling point and IPREDij was the individual prediction value. The residual error (ε) is normally distributed with a mean of 0 and a variance of σ2.
The impact of covariates on pharmacokinetic parameters was evaluated. The relationship between covariates and IIV for pharmacokinetic parameters were investigated graphically ahead of covariate assessment. In this dataset, both categorical variables (gender) and continuous variables (height, BW, age, disease duration, WBC, RBC, hemoglobin, PLT, NEUT, LYMPH, MONO, EO, BASO, TP, ALB, ALP, AST, ALT, TBIL, DBIL, urea, creatinine, UA, CK, CK-MB, LDH, HBDH, HSTN-I, GFR) were tested. Categorical covariates were included in the population model with the use of indicator variables and the impact of the categorical covariates on each parameter was tested using power function. Continuous covariates were centered at their median values, and the impacts of each covariate on parameters were evaluated using power functions. The stepwise forward addition and backward elimination were applied to develop the covariate model. In the forward addition, a covariate was retained if it resulting in a reduction the OFV > 6.64 (P < 0.01, df = 1). In the backward elimination, each covariate was left out of the full model built by forward addition one at a time. A covariate was retained if the elimination of the covariate lead to an increase of OFV > 10.83 (P < 0.001, df = 1).