Chromatographic separations were performed in an LC-20AD analytical pump (Shimadzu) using a SIL-20A HT autosampler (Shimadzu). A Luna 5-μm C18 100A 150 × 4.6 mm column (Phenomenex) was employed as stationary phase, and acetonitrile (B) and water with 0.1% formic acid (A) were employed as mobile phase. The mobile phase flow rate comprised of 65% of B was set at 1.7 ml min⁻¹, resulting in a total run of 2.8 min. Autosampler temperature was set at 22°C, while the injection volume was 20 µl.

Mass spectrometry analyses were performed using a Quattro Micro (Micromass) mass spectrometer equipped with an electrospray ionization (ESI) source. Analyses were performed in positive ionization mode using nitrogen as desolvation gas. ESI source parameters were set as follows: source temperature of 105°C, desolvation temperature of 480°C, desolvation flow of 800 l h⁻¹, and capillary voltage of 3 kV. Rosuvastatin was detected by a multiple reaction monitoring (MRM) transition of 482.18 > 258.14 using a cone voltage of 35 V, and diazepam was detected by a MRM transition of 285.20 > 193.30 using a cone voltage of 30 V. Data were acquired using MassLynx 4.1.

A stock solution of rosuvastatin was prepared by accurately weighing and dissolving the standard in acetonitrile/water (80:20 v/v) to a concentration of 100 μg ml⁻¹. Then, the solution was sequentially diluted to working solutions of 1, 0.152, and 0.0024 μg ml⁻¹. IS diazepam was dissolved in acetonitrile/water (80:20 v/v) to obtain a 100-μg ml⁻¹ stock solution and then diluted to obtain a 0.015-μg ml⁻¹ working solution. Method validation, calibration curve, and quality control analysis were carried out using blank plasma calibration standards including 0.2, 1, 5, 10, 20, 30, 40, and 50 ng ml⁻¹ blanks and quality control (QC) of 0.6, 25.0, and 38 ng ml⁻¹ samples prepared by spiking the above-described working solutions of rosuvastatin and IS in a blank biological matrix. Stock and working solutions were stored in a refrigerator and allowed to equilibrate (30 min) at room temperature before use.

Three hundred microliters of plasma were aliquoted into a microcentrifuge tube. Then, 25 µl of IS solution were added to the plasma. Next, 25 µl of HCl 1 M as well as 1,000 µl of a diethyl ether/dichloromethane (70:30 v/v) were added to the same tube. The content was shaken for 5 min and then centrifuged for 10 min, at 132,000, at 4°C. After that, the supernatant was transferred to another tube and dried under nitrogen stream. An acetonitrile/water solution (80:20 v/v) was used for resuspension of the dried content, which was shaken and finally transferred to an insert before injection.

The analytical method was validated for the following parameters: specificity, carryover, matrix effect, calibration curve, accuracy, precision, and stability.

Specificity was ascertained by analyzing blank human plasma samples from six individuals and comparing the chromatograms with chromatograms from blank human plasma spiked with rosuvastatin in the lower limit of quantification (LLOQ) concentration and IS.

Carryover was ascertained by analyzing three injections of the same blank sample, one before and two after an injection of a sample in the upper limit of quantification (ULOQ) concentration. Results were compared with LLOQ chromatograms.

Matrix effect was ascertained by spiking eight different extracted blank human plasma samples with rosuvastatin at QC concentrations and the IS. Peak areas of extracted spiked samples were compared to those of standard solutions.

The calibration curve was prepared using six different blank human plasma samples (three for weighting tests and three for linearity tests) spiked with concentrations at eight levels ranging from 0.2 to 50 ng ml⁻¹. Analyte concentrations in samples were calculated by linear regression equation, typically described by equation y = ax + b where y corresponds to the analyte/IS peak-area ratio and x corresponds to the ratio of rosuvastatin to IS concentration. Due to the range of the calibration curve and lower value of the sum of the relative errors of the nominal values of the calibration versus its values obtained by the curve equation, the weighting factor of reciprocal concentration squared (1/x ²) was applied.

Intra- and inter-batch accuracy and precision were evaluated at three different levels (0.6, 25, and 38 ng ml⁻¹) of QC samples in quintuplicate in three different batches. The accuracy of the method was expressed as relative error (RE), whereas precision was obtained by calculating the within- and between-run coefficient of variation (CV). The acceptance criteria for each quality control were that CV must not exceed 15% for QC and 20% for LLOQ.

Freezing and thawing stabilities for rosuvastatin in human plasma samples were ascertained after four cycles and the analytical process conducted at low and high QC concentrations. Samples were frozen at −70°C in four cycles of 12, 24, 36, and 48 h. Autosampler stability was investigated over a 24-h storage period in the autosampler tray with low and high quality control concentrations. Long-term stability was also evaluated over a 160-day storage period at −70°C. The acceptance criteria for each quality control were that CV and accuracy must not exceed 15%.

The validated LC-MS/MS method was applied to measure rosuvastatin concentrations to obtain additional pharmacokinetic parameters.

Pharmacokinetic parameters were calculated from plasma levels employing a non-compartmental statistic using WinNon-Lin 8.3 software (Pharsight, United States). Following Food and Drug Administration (FDA) guidelines, blood samples were drawn up to a period of three to five times the terminal elimination half-life (t _1/2) and the mean AUC_0–t/AUC_0–∞ ratio was required to be higher than 80%. The area under the concentration–time curve (AUC_0–t) was calculated by the trapezoidal method. The total area under the curve (AUC_0–∞) was obtained up to the last measurable concentration, and extrapolations were performed using the last measurable concentration and the terminal elimination rate constant (K _e). The terminal elimination rate constant, K _e, was estimated from the slope of the terminal log₁₀ transformed exponentially and multiplied by −2.303 phase of the plasma of rosuvastatin concentration–time curve (by means of the linear regression method) adjusted in the three last values. The terminal elimination half-life, t _1/2, was then obtained as 0.693/K _e. The C _max and the time to reach maximum plasma concentration (T _max) values were determined by visual inspection of the plasma rosuvastatin concentration–time profiles. Results are presented as mean ± standard deviation (SD).

Chromatographic analyses were performed in an ACQUITY FTN H Class (Waters) liquid chromatograph. Separations were carried out using an ACQUITY CSH C18 column (Waters) with dimensions 2.1 × 100 mm × 1.7 µm as stationary phase using acetonitrile (B) and water with 0.1% formic acid (A) as mobile phase. Flow rate was set at 0.4 ml min⁻¹. Segmented gradient was applied as follows: 0–2 min, 10% of B; 7 min, 90% of B kept for 2 min; and 9–11 min returning to initial conditions and kept for 2 min for column re-equilibration resulting in a 13-min run. Autosampler temperature was set at 20°C, whereas the injection volume was 0.5 µl for positive mode and 5 µl for negative mode.

High-resolution mass spectrometry analyses were performed using a XEVO-G2XSQTOF (Waters) equipped with an ESI source. Analyses were performed in both positive and negative ionization mode using nitrogen as the desolvation gas. ESI source parameters were set as follows: source temperature of 140°C, desolvation temperature of 550°C, desolvation flow of 900 l h⁻¹, capillary voltage of 3 kV, sampling cone of 30 kV, and cone gas flow of 10 l h⁻¹ for positive mode and capillary voltage of 2.5 kV, sampling cone of 40 kV, cone gas flow of 50 l h⁻¹ for negative mode.

Spectra were acquired within an acquisition mass range of 50–1,700 Da using a data-independent MS^E approach. MS^E fragmentation settings were applied using 6 V of collision energy for low-energy scan and a 15–30-V ramp for high-energy scan. Leucine encephalin (Tyr-Gly-Gly-Phe-Leu, formula C₂₈H₃₇N₅O₇, molecular weight = 555.2693 g mol⁻¹ at 200 pg μl⁻¹ in acetonitrile/water 1:1 v/v) was used as lockmass to ensure exact mass measurements during data acquisitions. A 0.5-mM sodium formate solution was used for instrument calibration. MassLynx 4.1 was used for mass spectrometer control and data acquisition management. Samples were randomly analyzed. Pooled quality control samples were analyzed at the beginning of the batch in order to equilibrate the system as well as sample injections at regular intervals during the batch to monitor instrumental drifts that may occur. Additionally, at the end of the batch, serially diluted QC samples were analyzed in order to verify features that follow dilution trends and filtered.

A 50-µl aliquot of each sample was taken, and 100 µl of methanol containing the IS p-bromo-phenylalanine (200 µM) were added, both in the same centrifuge tube. Then, samples were shaken for 10 min for protein precipitation and centrifuged at 14,000 rpm for 1 min at 8°C. Finally, 100 µl of supernatant was transferred to a vial for instrumental analysis. Samples were randomly prepared. A pooled QC sample was prepared by mixing equal small aliquots of each study sample, distributed equally in regular intervals across the sample preparation batch, and extracted as such. In addition, five-level serially diluted QC samples were prepared.

Raw data acquired by UPLC-QTOF-MS^E were directly imported to Progenesis QI (Nonlinear Dynamics) for peak detection and deconvolution, alignment, retention time correction, data filtering, and MS^E-based ion annotation and identification to generate a suitable data matrix table for statistical analysis. For the feature-identification step, MassBank of North America (MoNA)-containing metabolomics libraries such as HMDB, LipidMaps, LipidBlast, and MassBank were used.

QC sample acquisition allowed for quality data checking and cleaning, and molecular features that showed relative standard deviation (RSD) > 25% after data cleaning and normalization were stripped out from the multivariate data modeling.

Machine learning modeling was performed by means of Elastic Net (sklearn.linear_model.ElasticNet), with the linear regression with combined L1 and L2 penalty priors as a regularizer, which penalizes models based on their complexity favoring simpler and more generalizing models. This algorithm was applied as estimator for AUC and C _max predictions using Python 3.7.10 with scikit-learn (Pedregosa et al., 2012) version 0.22.2. Other Python programing standard packages for data manipulation, such as Pandas (McKinney, 2010) and Numpy (Kristensen and Vinter, 2010), and data visualization, such as Matplotlib (Hunter, 2007), Seaborn (Waskom, 2021) and ScyPy (Virtanen et al., 2020), were also employed.

Results were mainly assessed in terms of explained variance R ² (sklearn.metrics.explained_variance_score), main absolute percentage error (MAPE) (sklearn.metrics.mean_absolute_percentage_error), and root mean squared error (RMSE) (sklearn.metrics.mean_squared_error, squared = False). For overfitting, checking, and model validation, the same metrics were employed obtained by leave-one-out cross-validation (LOOCv) predicted data for each validation sample [sklearn.model_selection.cross_val_predict, cv = LeaveOneOut()]. Models were refined based on selection of features with non-zero coefficients where, among all initially selected features, different number of features starting from the most significant, and by steps of one, each feature was evaluated so that in the end all significant features were evaluated. For each iteration, the corresponding quality metrics were collected for further evaluation. The optimal numbers of significant features were further used in the process of model hyperparameter tuning (sklearn.model_selection.GridSearchCV). The code is available in the following Github repository: https://github.com/GustavoHBDuarte/Metabolomics_PK_proj.

Specificity assessment results show that the developed method for the analysis of rosuvastatin using the selected mass transitions for MRM function was selective enough for quantitative purposes, with no apparent interferences either from endogenous compounds or from matrix effects around the analyte and IS retention times. Figure 1 portrays the chromatograms obtained from blank plasma samples and blank plasma samples spiked with rosuvastatin (0.2 ng ml⁻¹) and IS.

To evaluate the presence of carryover effects, three injections of blank samples were analyzed, one before and two after ULOQ. No signal of analyte was detected in any of the blank injections.

For matrix effect evaluation, the mean values obtained for QCs and IS were 108.4% (SD = 10.9%) and 104.4% (SD = 8.6%), respectively, demonstrating that no apparent matrix effect affects the rosuvastatin determination in human plasma samples employing the developed method.

The calibration curves were plotted as the peak area (rosuvastatin/IS)/drug concentration ratio. Linearity was achieved over the concentration range of 0.2–50 ng ml⁻¹ with the correlation coefficient (R ²) greater than 0.99 for all the curves. The RE of the mean concentration measured ranged from −6.17% to 9.50%, and the RSD ranged from 3.02% to 14.46%. Table 1 portrays the overall results for linearity evaluation.

Accuracy and precision were calculated in terms of intra- and inter-batch variation at QC levels, and the corresponding results are shown in Table 2. Intra- and inter-batch precision (RSD) ranged from 5.87 to 14.16 and 5.84 to 10.59, respectively, whereas intra- and inter-batch accuracy ranged from −0.31% to 5.83% and −3.04% to 1.88%, respectively, indicating that the method is reliable and reproducible within its analytical range.

The four-cycle freezing and thawing stability (−70°C to room temperature) evaluations have shown that the analyte is stable in human plasma at QC levels for all cycles. Additionally, no analyte degradation was detected in an autosampler tray stability assessment over a 24-h storage period, with the measured analyte values ranging from 83.17% to 111.58% of their corresponding nominal values. Moreover, for long-term evaluation, results have shown that rosuvastatin was stable over a 160-day storage period at −70°C, with measured analyte values ranging from 90.42% to 97.17% of the nominal values.

The validated analytical method was then applied to measure rosuvastatin concentration in healthy volunteers in order to obtain the main pharmacokinetic parameters. This study was conducted with 40 volunteers after a single oral dose (20 mg) of the drug. Figure 2 portrays the typical plasma concentration versus time profiles. Plasmatic concentrations of rosuvastatin ranged within the standard curve and remained above the LLOQ of 0.2 ng ml⁻¹ for the entire sampling period.

The observed AUC value from time 0 to the last sampling time (AUC_0–t) was 154.13 ± 63.06 ng h/ml, whereas the AUC from 0 to ∞ (AUC_0–∞) was 162.85 ± 64.96 ng h/ml. The elimination half-life (t _1/2) was 15.10 ± 7.91 h.

The observed C _max was 19.53 ± 8.64 ng ml⁻¹, whereas the T _max was 4.18 ± 0.72 h. Moreover, the AUC_0–t/AUC_0–∞ mean ratio was higher than 80%.

The refined Elastic Net model was then applied to predict AUC values of study volunteers from predose endogenous metabolite profiling acquisition. Figure 3 portrays the predicted-against-experimental AUC values for ESI (+) and ESI (−) data where both training and LOOCv testing datasets are evaluated simultaneously. Ideally, in a good prediction model, all the samples should fall as close as possible to the regression line. The R ² coefficient scales the extension of variance in AUC values that can be explained by the model. For ESI (+) data, R ² achieved 1.00, whereas for ESI (−) data, R ² achieved 0.98 for the training set. Although training R ² results were close to 1.0, it is crucial to consider the results of the LOOCv external set in order to evaluate the model’s ability to predict unseen data so that it can be properly validated. The LOOCv R ² Cv was 0.93 for ESI (+) and 0.86 for ESI (−). The closeness between training and Cv for R ² indicates both models’ good ability to predict AUC values for external sets of samples.

Furthermore, root mean square error of cross-validations (RMSE Cv) and mean absolute percentage error of cross-validations (MAPE Cv) were also calculated. These parameters specify error measurements for fitted model observations. For ESI (+), RMSE Cv reached 19.03 ng h/ml and MAPE Cv reached 10.12%, while for ESI (−), RMSE Cv reached 27.92 ng h/ml and MAPE reached Cv 14.28%.

From predicted and true AUC values, it is possible to obtain a residual plot. In the dispersion diagram, it is possible to note both training and LOOCv testing set samples similarly distributed around the ordinal least squares (OLS) R ² regression line. Moreover, samples with higher AUC values had higher residuals. Considering that the frequency of observations with higher AUC values was substantially lower than the frequency of observations with lower AUC values, the selected machine learning algorithm possibly turned out to be trained differently across this AUC range, and apparently, this difference reflected on the calculated residuals. The OLS R ² regression line value for LOOCv residuals was 0.25 for ESI (+), meaning that 25% of the variation in AUC-predicted values is explained by variation in residuals shown in the dispersion diagram, whereas for ESI (−), the OLS R ² achieved 0.19. Furthermore, shared Y-axis reveals higher residuals for the ESI (−) model compared to the ESI (+) model. Figure 4 portrays the residual evaluation in the form of a dispersion diagram.

Similarly, the information provided by residuals was used to generate both histograms and probability plots to evaluate the distribution of frequencies for residuals and obtain a more in-depth visualization. Figure 5 portrays the results. A histogram inspection suggests residuals around zero as being the most frequent values for both the ESI (+) and ESI (−) models, whereas the probability plot reveals that, although not falling straight to the normality Gaussian fit line, the distribution of residuals against theoretical quantiles follows a similar trend for both ESI (+) and ESI (−), showing the consistency of the obtained data.

Regarding machine learning modeling, applying exclusively clinical features, unlike the performance of Elastic Net model using metabolomics data as predictors, it was not possible to achieve similar results using volunteers’ clinical features as pharmacokinetic parameter predictors for AUC. RMSE Cv reached 63.92 ng h/ml, whereas MAPE Cv reached 33.00%. The LOOCv R ² coefficient was 0.24. Predicted AUC values from the ESI (+), ESI (−), and volunteers’ clinical feature modeling were also compared to each other from the perspective of the frequencies of their corresponding predicted values. Figure 6 demonstrates a clear representation of the extension of the lesser predictability of the clinical feature model compared to the ESI (+) and ESI (−) models by using frequency plots. Results demonstrate that the clinical feature model generates an evidently higher difference in the distance of both training and LOOCv testing frequencies of predicted values from true values compared to the same measure of both the ESI (+) and ESI (−) models, showing its considerably lower predictive ability compared to information provided by metabolomics.

The final evaluation of AUC modeling included an integrated model applying both metabolomics [ESI (+) and ESI (−)] and clinical features. The information provided by the combination of both sources did not substantially improve the performance of predictions for ESI (+). Instead, for ESI (−), a slight decrease was observed in the quality metrics of the model. LOOCv R ² reached 0.94 and RMSE Cv and MAPE Cv reached 18.32 ng h/ml and 10.21%, respectively, for ESI (+), whereas LOOCv R ² reached 0.85 and RMSE Cv and MAPE Cv reached 28.21 ng h/ml and 14.39%, respectively, for ESI (−).

Refined ESI (+) and ESI (−) datasets were also employed to predict C _max values of the same study volunteers from predose endogenous metabolite profiling experiments, and their corresponding performance was assessed using the same metrics used for AUC predictions. The same trend observed for AUC predictions, whereby the ESI (+) model performed better than the ESI (−) model, was also observed for C _max models. Figure 7 portrays predicted-against-experimental C _max values for ESI (+) and ESI (−) for each observation. For ESI (+), training R ² achieved 0.99 and LOOCv testing R ² was 0.94, whereas for ESI (−), training R ² achieved 0.9 and LOOCv testing R ² achieved 0.79. RMSE Cv and MAPE Cv for ESI (+) were 2.54 ng ml⁻¹ and 11.56%, respectively, whereas ESI (−) RMSE Cv and MAPE Cv were 4.55 ng ml⁻¹ and 23.24%, respectively.

Figure 8 displays the residual analysis for C _max predictions as a dispersion diagram. The same trend seen for AUC predictions is also seen for C _max predictions in the corresponding plot where residuals for the ESI (−) model were higher compared to the ESI (+) ones. Although it can be observed that unlike the AUC predictions, there was a sample with a higher residual value [close to −10 in the ESI (−) model]. The sample in question was not the sample with the highest C _max value. For ESI (+), the samples with higher residuals were not necessarily the samples with higher C _max values either. The observations with higher C _max values were also less frequent than with the AUC model. The fact that, for the C _max model, the range of numerical values to predict is lower and samples are closer to each other could allow the C _max model to train better than the AUC model for predicting the numerical value of samples and spreading the residuals more evenly throughout the samples. Besides, the OLS R ² regression line value for LOOCv residuals shown in the dispersion diagram was 0.02 for ESI (+) and 0.01 for ESI (−). Furthermore, for the distribution analysis of residuals, shown in Figure 9, the residuals with values close to zero were mainly the most frequent values. The shared X-axis of histograms also allowed for the visualization of the lesser spreading of residuals of the ESI (+) model compared to the ESI (−). Similarly, the shared X-axis of the probability plots shows, as in the AUC models, the similar distribution of data regarding residuals against theoretical quantiles, regardless of ionization mode, confirming the consistency of obtained data as the data points apparently seem to fall closer to the Gaussian regression line compared to the AUC model and, therefore, a more normal-like distribution of residuals around zero.

Regarding C _max predictions using volunteer’s clinical features, a similar trend observed for AUC models was also seen for C _max modeling with poorer prediction ability when comparing to the model applying only metabolomics data. The LOOCv R ² coefficient reached 0.08, whereas RMSE Cv and MAPE Cv reached 9.77 ng ml⁻¹ and 34.77%, respectively, showing substantially lower correspondence between clinical features and C _max values. Figure 10 shows a similar trend for AUC predictions using volunteer’s clinical features, where predictions for both training and LOOCv testing demonstrated a considerably lower ability to predict pharmacokinetic values with a clear distortion in the frequencies of predictions for both training and LOOCv testing models. A more detailed visualization in the distribution plot of predictions for the C _max ESI (+) model reveals a closer overlapping between the frequencies of training and LOOCv testing predictions compared to the frequencies of the C _max ESI (−) and even closer if compared to the AUC ESI (+) and AUC ESI (−) training and LOOCv predictions. Though there is no apparent overfitting for any of these models, this described observation is consistent with the lower distance between R ² and LOOCv R ² for the C _max ESI (+) model compared to the same measurement of the C _max ESI (−), AUC ESI (+), and AUC ESI (−) models.

In addition, regarding the integrated model, the same behavior observed for AUC predictions was observed for C _max predictions with the application of both clinical and metabolomic information in the same dataset not providing significant improvements over the model using only metabolomic information. Supplementary Table S2 in the Supplementary Material data file provides an overall visualization of the performance metrics for all the evaluations carried out in this work.