Based a Machine Learning Approach to Investigate the Factors Influencing Nirmatrelvir/ritonavir Exposure in Human Plasma: A Multicenter, Observational Study

1. Introduction

COVID-19, caused by severe acute respiratory syndrome coronavirus 2019 (SARS-CoV-19), poses a significant threat to global public health systems. Antiviral drugs, including remdesivir [1], lopinavir/ritonavir [2], molnupiravir [3], and nirmatrelvir/ritonavir (N/R) [4], have demonstrated promising therapeutic effects against COVID-1. N/R, an oral antiviral drug manufactured by Pfizer in the United States, is considered an ideal treatment for mild to moderate COVID-19 and severe disease risk factors [5,6]. N/R consists of two active substances: nirmatrelvir and ritonavir. Nirmatrelvir is a potent inhibitor of SARS-CoV-2 protease, which hinders the processing of viral polyproteins and prevents viral replication. Ritonavir, on the other hand, primarily inhibits the cytochrome P450 3A4 enzyme (CYP3A4), preventing premature metabolic inactivation of nirmatrelvir and increasing its serum concentrations [7]. Thus, the combination of nirmatrelvir and ritonavir can enhance the effectiveness of nirmatrelvir in alleviating clinical symptoms in patients with COVID-19.

Despite the effectiveness of N/R in combating COVID-19, the inhibition of CYP3A metabolism by ritonavir can lead to a rapid increase in the concentration of drugs such as tacrolimus and cyclosporine, which are also metabolized by CYP3A4. This increase can result in adverse effects such as brain damage, epilepsy, kidney damage, and even death [8,9,10,11]. N/R is not recommended for patients with severe renal or hepatic impairment. Therefore, for patients taking N/R, especially those with special conditions, timely monitoring of drug concentration is crucial for effective treatment. Therapeutic drug monitoring (TDM) has become an important tool for optimizing the use of biopharmaceuticals. It is primarily used to monitor drugs with challenging target concentrations, significant pharmacokinetic variability, and adverse reactions [12,13]. Therefore, based on previous experience [14,15], TDM may be an optimized tool for drug therapy to manage drug-drug interactiong (DDI), specifically for patients using potential related medications.

To ensure the effectiveness of N/R and minimize its side effects, TDM was utilized to monitor the blood concentrations of nirmatrelvir and ritonavir during COVID-19 treatment with N/R. In this study, we analyzed patients who received N/R treatment and underwent TDM monitoring in two study center (Changxing People's Hospital and Second Affiliated Hospital of Xi’an Jiaotong University). We identified several factors that could potentially influence N/R blood levels, providing a foundation and reference for the clinical use of N/R.

2. Patients and Methods

2.1. Study Design

This retrospective multicenter study included all patients who received N/R treatment at Changxing People’s Hospital and Second Affiliated Hospital of Xi’an Jiaotong University from December 2022 to June 2023. Informed consent forms were signed by all patients. All patients were administered two tablets of nirmatrelvir and one tablet of ritonavir orally or nasally every 12 hours. Patients are advised to take nirmatrevir and ritonavir simultaneously. The maintenance dose was 300 mg administered for 5 days.

2.2. Patients Enrollment

Inclusion criteria: age≥18 years; patients adhered to the prescribed medication regimen. Blood samples were collected from all patients on the third day after drug administration, 30 minutes prior to drug administration. At least one sample from each patient was tested for TDM. If patients discontinued medication for reasons other than the study, the time of re-administration was considered as the initial treatment. The primary disease under investigation was neocoronitis. Exclusion criteria: incomplete clinical data recording, failure to meet blood sample collection and testing requirements.

2.3. Clinical Data and Definitions

At baseline, we collected and recorded the following information from patients: sex, age, department, diagnosis, smoking and alcohol abuse history, underlying diseases (hypertension), blood routine index, liver function indicators (alkaline phosphatase [ALP], glutamyl transpeptidase [GGT], aspartate transaminase [AST], and alanine transaminase [ALT]), renal function indicators (creatine kinase [CK], creatinine, and estimated glomerular filtration rate [eGFR]), clinical symptoms, and use of CYP inducers and inhibitors.

2.4. Quantification of Nirmatrelvir and Ritonavir Plasma Concentrations

2.4.1. Chromatographic and mass spectrometric conditions

A methodology for the analysis of N/R in the body has been established using UPLC-MS/MS [16]. Specifically, we utilized an ACQUITY UPLC BEH C18 (2.1 mm×50 mm, 1.7 μm) column, with 100% acetonitrile as the A-phase and water containing 0.1% formic acid as the B-phase. To enhance sample separation and improve separation efficiency, we performed gradient elution on the samples. The gradient elution consisted of the following conditions: 0-1.2 min, 22% A; 1.2-1.9 min, 100% A; 1.9-2 min, 100% A; 2-3 min, 22% A. The elution was carried out at a flow rate of 0.3 mL/min, with a column temperature of 45℃ and an injection volume of 2μL.

The electrospray ionization source (ESI)+ mode was utilized, with a capillary voltage of 3000 V, a desolvation temperature of 350℃, an N2 flow rate of 600 L/h, and a scanning mode of MRM. Positive and negative ions were alternated for scanning and analysis, and the parameters for mass spectrometry detection are provided in the table below (Table 1).

2.4.2. Preparation of Reserve Solution and Working Solution

The Nirmatrelvir and ritonavir standards were accurately weighed. Nirmatrelvir and its internal standard, nirmatrelvir-D9, were dissolved in acetonitrile, while ritonavir and its internal standard, ¹³C, ²H₃--ritonavir, were dissolved in acetonitrile: water (50:50, V/V). Nirmatrelvir and ritonavir were prepared as separate standard stock solutions with a concentration of 1.00 mg/mL each. Nirmatrelvir-D9 and ¹³C, ²H₃-ritonavir were dissolved separately, mixed, and then diluted with acetonitrile: methanol (3:1, V/V) containing 0.1% formic acid. This resulted in a mixed internal standard working solution with concentrations of 1.00 μg/mL for nirmatrelvir- D9 and 0.50 μg/mL for ¹³C, ²H₃-ritonavir.

Take the appropriate amount of nirmatrelvir and ritonavir standard stock solution. Dilute them step by step with acetonitrile: water (3:1, V/V) containing 0.1% formic acid to obtain a series of mixed standard working solutions. The concentrations of nirmatrelvir in the solutions are 300, 100, 50, 20, 10, 5, and 1 μg/mL, while the concentrations of ritonavir are 50, 40, 20, 10, 5, 1, and 0.5 μg/mL. Additionally, prepare low concentration quality control (QC) samples with nirmatrelvir at 5 μg/mL and ritonavir at 1 μg/mL, medium concentration QC samples with nirmatrelvir at 10 μg/mL and ritonavir at 10 μg/mL, and high concentration QC samples with nirmatrelvir at 80 μg/mL and ritonavir at 40 μg/mL.

2.4.3. Preparation of Standard Curve and Quality Control

The standard curve and QC were prepared by adding 10 μL of a mixed standard working solution containing nirmatrelvir and ritonavir to 90 μL of blank human plasma. The standard curve for nirmatrelvir ranged from 30, 10, 5, 2, 1, 0.5, 0.1 μg/mL, while the standard curve for ritonavir ranged from 5, 4, 2, 1, 0.5, 0.1, 0.05 μg/mL. The low concentration levels were 0.5 μg/mL for nirmatrelvir and 0.1 μg/mL for ritonavir. The medium concentration levels were 1 μg/mL for nirmatrelvir and 1 μg/mL for ritonavir. The high concentration levels were 8 μg/mL for nirmatrelvir and 4 μg/mL for ritonavir. The lower limit of quantification was 0.1 μg/mL for nirmatrelvir and 0.05 μg/mL for ritonavir.

2.4.4. Sample Pretreatment

Transfer 100 μL of plasma sample into a tube, then add 300 μL of a mixed internal standard working solution. Mix the contents thoroughly and centrifuge at 12000 rpm/min for 5 minutes. Finally, take 2 μL of the supernatant for analysis.

2.4.5. Statistical Methods

To ensure the stability of the subsequent model, we conducted an initial assessment of the sample (n = 139) for multicollinearity. Through Spearman correlation analysis, we identified highly correlated features and selectively eliminated statistically significant ones to mitigate the effects of multicollinearity.

To construct a logistic regression model, we treated prognostic outcome as the dependent variable and nirmatrelvir blood concentration as the independent variable. The purpose of this model was to analyze the correlation between the two variables and determine the optimal blood concentration thresholds for distinguishing prognostic outcomes. To assess the predictive performance of the model, we calculated the receiver operating characteristic curve (ROC) and the area under the curve (AUC)[17]. To tackle the imbalance in prognostic outcomes, we employed cost-sensitive learning by incorporating inverse frequency weights derived from the distribution of prognostic categories. This technique aimed to improve the model’s sensitivity to less frequent categories. Furthermore, we improved our model through an iterative threshold adjustment process. In each cycle, we systematically eliminated the sample with the highest nirmatrelvir blood concentration. Subsequently, we reevaluated the model to assess its impact on AUC and the threshold value. This iterative exclusion aimed to optimize the threshold for prognostic discrimination, ensuring resilience against outliers or extreme values.

We employed a random forest model [18] to determine the factors influencing blood drug concentration. All samples were divided into a training cohort (n = 112) and a validation cohort (n = 27) at a 4:1 ratio using R software (version 4.3.0). The random forest model consisted of 1,000 trees and underwent 10-fold cross-validation. Features that showed a positive percentage increase in mean squared error (IncMSE) during training were considered significant and chosen for the XGBoost prediction model. We evaluated the performance of the XGBoost model by generating learning curves through iteration over various hyperparameter combinations. The learning curves helped determine the optimal combination to fine-tune the model. Metrics such as Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), R-squared (Coefficient of Determination), and Mean Squared Error (MSE) were used to assess the accuracy and calibration of the model. Additionally, we confirmed the model’s internal validity through five-fold cross-validation.

3. Results

3.1 Patient Summary

This study included a total of 139 subjects, all of whom were initially diagnosed with SARS-CoV-2 infection. The age of the subjects ranged from 18 to 97 years, with a mean age of 73.42 ± 21.36 years. A total of 77 male and 62 female patients were included. Additional baseline laboratory findings are presented in Table 2.

3.2. Changes in Clinical Symptom Score

We assessed the severity of common symptoms in COVID-19 patients, including fever, cough, nasal congestion, headache, muscle ache, nausea, chills, vomiting, and diarrhea, using a quantitative scoring system (0-none, 1-mild, 2-moderate, 3-severe). As depicted in Figure 1, the clinical symptom scores of nearly all patients showed a gradual decrease over the course of N/R administration. The aforementioned symptoms exhibited a gradual decrease in patients who received continuous N/R administration. These findings suggest that N/R exhibits efficacy in the treatment of SARS-CoV-19, and the timing of administration may influence its effectiveness.

3.3. Feature Correlation Analysis

Before performing feature selection, we conducted a correlation analysis on all variables to identify and eliminatehighly correlated features. Highly correlated features were defined as those with an absolute correlation coefficient greater than 0.5 and a P-value less than 0.05. Figure 2A shows the Spearman correlation heatmap for all features, and Figure 2B displays the heatmap after removing features with correlations beyond the defined threshold. Stronger correlations are represented by darker shades. We observed a strong correlation between the concentrations of nirmatrelvir and ritonavir (Correlation: 0.6063156, P＜0.001). Therefore, we decided to analyze only the concentrations of nirmatrelvir in the subsequent analyses.

3.4. Nirmatrelvir Blood Concentration Predicts The Prognosis

We evaluated the predictive value of Nirmatrelvir blood concentration for COVID-19 patient outcomes using the receiver operating characteristic (ROC) curve. The area under the ROC curve (AUC) was 0.467, indicating that the model’s predictive ability is similar to random guessing. This finding suggests that Nirmatrelvir blood concentrations may not have a significant correlation with prognostic outcomes, indicating limited usefulness in predicting patient prognosis (Figure 3).

3.5. Feature Selection from Random Forest Model

The Random Forest model was used to identify the features that had the greatest impact on our analysis. This model ranked the features based on their importance. The most important feature identified was the eGFR, which is an indicator of renal function. It was followed by creatine kinase, AST, ALT, Lymph, PCT, the use of the cytochrome P450 enzyme inhibitor omeprazole (CYP-4), and the patient’s smoking and drinking history. These results suggest that factors such as renal and liver function, lymphocyte and platelet counts, and lifestyle choices like smoking and alcohol consumption may significantly affect the metabolism and elimination of N/R, as shown in Figure 4A.

Further analysis was performed to investigate the relationship between these features and blood drug concentration. Scatter plots were created to visually represent the distribution of each feature’s value in relation to blood drug concentration. Linear regression was then applied to these distributions, and the resulting regression lines were used to illustrate the linear relationships between each feature and blood drug concentration, as depicted in Figure 4B.

3.6. Learning Curve Analysis for Hyperparameter Optimization

After selecting important features using the random forest model, we trained the XGBoost model with these features to determine the optimal hyperparameters for the model’s generalization performance. We analyzed the learning curves to observe how the model’s error, measured by RMSE, changed with increasing boosting iterations under different settings of eta (learning rate) and subsample (random sampling rate of samples), as shown in Figure 5.

The results shown in Figure 5 indicate that the model, with an eta of 0.25 and a subsample of 1, consistently achieved a low and stable RMSE. This suggests that the model had an appropriate learning rate and effectively utilized the dataset. Additionally, a max tree depth of 6 was found to be effective in maintaining the model’s performance without overfitting, even as the number of iterations increased.

In conclusion, we chose a hyperparameter set that included eta = 0.25, subsample = 1, and max tree depth = 6. The training process was terminated after 500 iterations, as it was found to yield the best validation performance. This particular hyperparameter configuration allowed us to effectively utilize the selected features and capture the intricate relationships within the data, while avoiding overfitting.

3.7. XGBoost Model

After training the XGBoost model with the optimally selected hyperparameters, we assessed the importance of each feature in the model. The model output includes three metrics of feature importance: Gain, Cover, and Frequency. ‘Gain’ represents the average gain of splits that utilize the feature, ‘Cover’ corresponds to the average coverage of splits that utilize the feature, and ‘Frequency’ indicates the proportion of times a feature is used in splits across all trees.

The results, depicted in Figure 6, emphasize that eGFR has the highest gain among the features, implying its significant influence on model predictions. Additionally, liver enzymes ALT and AST, CK, Lymph, and PCT exhibit notable importance, signifying their relevance in the model’s decision-making process. Analyzing these metrics reveals the biological factors that are most predictive of blood drug concentrations based on the trained model (Figure 6).

We trained the XGBoost model with the selected features and hyperparameters, and then evaluated its predictive performance using the validation set. The model’s predictions were assessed using various statistical metrics to measure prediction error and goodness of fit. The mean squared error (MSE) was 1.328, indicating the average squared difference between the observed actual outcomes and the model’s predictions. The root mean squared error (RMSE), which measures error in the same units as the response variable, was 1.152. The mean absolute error (MAE), representing the average absolute difference between observed and predicted values, was 0.717. Additionally, the model achieved an R-squared (R²) value of 0.779, suggesting that approximately 77.9% of the variability in the response variable was explained by the model. As shown in Table 3, these metrics collectively demonstrate a strong predictive performance, with the model accurately forecasting blood drug concentrations.

4. Discussion

N/R has been crucial in treating the SARS-CoV-2019 pandemic. N/R has been administered to adult patients with COVID-19 at risk of severe disease, as well as pediatric patients with moderate to mild disease (≥12 years of age weighing at least 40 kg). Studies have demonstrated that N/R use significantly reduces hospitalization and mortality rates in patients with COVID-19 [5,19,20]. However, N/R is not without flaws, and several studies have revealed that its use may lead to varying degrees of side effects in patients [10,21]. With the advancement of individualized medicine, TDM has emerged as a prominent precision medicine tool. TDM assists clinicians in tailoring dosage adjustments for patients’ medication, maximizing drug efficacy, and minimizing adverse drug reactions by considering individual differences in pharmacokinetics [22,23]. While TDM has been used for clinical monitoring of various drugs, no study has investigated the potential factors influencing N/R blood concentration in different patients using TDM. Therefore, it is promising and highly meaningful for TDM to guide individualized dosing in the clinical study of N/R.

In our baseline data, we observed that patients with a prognosis of worsening or death were significantly older compared to those with an improved or cured prognosis, which aligns with previous findings. Age is a widely recognized risk factor for COVID-19 [24,25]. Older adults exhibit weakened immune systems, making them more vulnerable to viral infections and prone to immune system overactivation caused by chronic inflammation, potentially resulting in pneumonia. Furthermore, compromised organ function and pre-existing chronic conditions in older adults play crucial roles in pneumonia development. The slowing down of the body’s metabolism with increasing age suggests that elderly individuals need regular monitoring of blood levels and adjustment of the therapeutic regimen when using N/R. Additionally, a wide range of laboratory biochemical indicators are considered the main risk factors for the severity and mortality of COVID-19. Our study found significant differences in eGFR, CRP, PCT, Cr, WBC, NE, and Lymph in both groups, which is consistent with the findings of a previous study [24]. Differences in biochemical parameters may be associated with the prognosis of patients. On one hand, they can alter the metabolism and absorption of drugs, and on the other hand, they can influence the response of the immune system. During inflammation, the immune system stimulates a rapid increase in WBC and NE to combat pathogenic bacteria. Additionally, Lymph is massively translocated to the site of inflammation to recognize and destroy virus-infected cells. Elevated levels of CRP and PCT, which are markers of acute inflammation, indicated the presence of severe inflammatory infections in our patients in the group with worsening prognosis or death. While a small percentage of patients on N/R still failed to recover from COVID-19, the majority of patients had a favorable prognosis, indicating the beneficial effects of N/R against COVID-19.

In the study conducted by Cao et al [26], all patients in the nirmatrelvir-ritonavir group achieved sustained clinical recovery. Clinical recovery was defined as the remission of all target symptoms related to COVID-19, with a total score of 0 or 1 for each symptom. Consistent with their findings, our study involved patients who received a maintenance dose of N/R for 3 days. During this period, there was a gradual decrease in clinical symptom scores. On the fifth day, the majority of patients (111 cases) had a total score of less than 2, indicating significant relief from symptoms. We can hypothesize that the extended duration of N/R administration allowed nirmatrelvir to maintain high drug concentrations and sustained efficacy in patients. Despite the initial indication of a relationship between nirmatrelvir and patient outcomes, the logistic regression modeling analysis unexpectedly showed weak correlation between prognostic outcomes and nirmatrelvir blood concentrations. This could be attributed to the older age and presence of severe underlying diseases, including cardiovascular disease (including hypertension), diabetes, and leukemia, among the patients. The poorer prognostic outcome of the patients was likely a result of the underlying disease, causing organ dysfunction, failure, or excessive immunocompromise, unrelated to SARS-CoV-2 infection.

N/R blocks the activity of the SARS-COV-2-3Cl protease, which is necessary for coronavirus replication, to exert its antiviral effect [27,28,29]. After excluding highly correlated individual features, Spearman correlation analyses revealed a high correlation between nirmatrelvir and ritonavir concentrations, indicating that low-dose ritonavir had an effect. Nimatrevir alone may be eliminated by the body due to lack of recognition. The combination of nimaltegravir and ritonavir slows down the metabolism and breakdown of PF-07321332, thereby prolonging its high activity and exerting antiviral effects. The random forest model screened eight important features that affect blood drug concentrations by comparing the training cohort (n = 112) with the validation cohort (n = 27). Interestingly, the most important feature was eGFR, followed by Creatine Kinase, PCT, AST, ALT, Lymph, CYP4, and smoke-drink. Previous studies reported significantly higher Creatine kinase concentrations in deceased patients compared to recovered patients [30]. Creatine kinase is a risk factor for patient severity and mortality [31,32]. In this study, we discovered for the first time that both Creatine kinase and PCT have significant effects on N/R blood concentrations, which is a meaningful finding. During inflammation, cells are destroyed and release creatine kinase into the bloodstream, leading to elevated concentrations. Elevated levels of creatine kinase in the blood indicate substantial damage to existing cells, which affects the absorption of N/R. PCT, which is normally expressed at low levels, is released in large quantities into the bloodstream during infection and is considered an important indicator for assessing the severity of severe infections. In our study, we found that PCT had a significant effect on N/R concentration. The PCT results in the good prognosis group showed mild infection in patients (0.5-2 ng/mL), while in the poor prognosis outcome group, PCT reached 3.12 ± 8.01 ng/mL, indicating severe infection in patients. N/R inhibits the replication of the SARS-CoV-2 virus, but it does not affect established viral loads. In severe cases, where the viral load is high, the inhibition of viral replication alone is limited. Nirmatrelvir is metabolized and eliminated from the body when it fails to exert its antiviral effect. Therefore, the range of PCT significantly affects the concentration of nirmatrelvir and ritonavir. The liver plays a crucial role in drug metabolism, and AST and ALT are commonly used to assess liver damage. Elevated AST and ALT levels indicate impaired liver function, leading to reduced metabolism of N/R and affecting its efficacy. Additionally, omeprazole (CYP-4) may also play a significant role in influencing the concentration of N/R. Omeprazole down-regulates immunoglobulin E (IgE) receptor sites on plasma cell-like dendritic cells, leading to increased production of interferon alpha (INF-alpha) and enhanced antiviral immunity. We hypothesize that this could be attributed to the alteration of immune status by omeprazole, which may affect the cellular metabolism of N/R through immunokines or immune cells. The primary function of lymph is to maintain immune system homeostasis and to identify and eliminate abnormal cells. The quantity of lymph is indicative of the body’s immune function, which can potentially impact cellular behavior and consequently alter N/R blood levels.

This study represents the first application of the XGBoost modeling algorithm to identify the features that influence the concentration of N/R in the blood. Through careful selection of hyperparameters, we optimized the model’s performance in validation, striking a balance between complexity and generalization to mitigate the risk of overfitting. Learning curves played a crucial role in this process, allowing us to iteratively fine-tune the model until we achieved the optimal configuration. The XGBoost model identified six features, namely eGFR, CK, AST, ALT, Lymph, and PCT, as significant predictors of N/R concentration levels in COVID-19 patients. These features were consistent with those identified by the Random Forest model, thereby strengthening the validity of our feature selection process. The identification of these key features highlights the potential physiological and metabolic factors that may impact the response of patients to N/R. Integrating these findings into clinical practice could improve the therapeutic monitoring of N/R, potentially leading to more personalized and effective treatment strategies for COVID-19 patients.

It is important to acknowledge the limitations of the current study. Firstly, we did not continuously monitor N/R blood concentrations using TDM, which may not provide an optimal approach for observing continuous changes in drug concentrations. Secondly, the sample size we utilized was relatively small, and we plan to conduct a larger study to further validate our findings. Additional research is necessary to gain a better understanding of the factors contributing to N/R in COVID-19 patients and to accurately predict patient outcomes.

5. Conclusions

We identified a set of variables through therapeutic drug monitoring (TDM) to influenced N/R exposure levels and developed a new machine-learning model to identify them. These variables offer valuable predictive information about patients’ N/R exposure levels and prognosis. Furthermore, we provide approximate ranges of nirmatrelvir and ritonavir concentrations, which can guide the establishment of standardized ranges for steady-state concentrations. Clinicians can refer to our findings to customize management strategies for patients, while researchers can utilize our findings to develop multivariate prognostic models for enhancing patient prognosis.