UDCA May Promote COVID-19 Recovery: A Cohort Study with AI-Aided Analysis

Real-world data and cohort characteristics

In this retrospective study, 120 patients were diagnosed with COVID-19 between December 10, 2022, and December 30, 2022, of which 5 patients were not hospitalized despite having a fever for four or more days. Following the exclusion criteria, a total of 115 patients were retained for the univariate recovery analysis with adequate individual clinical data and documented diagnosis of COVID-19, 65 of whom had a claim for ursodeoxycholic acid (UDCA), while 50 did not receive any UDCA treatment. Among the UDCA-treated group, 29 (44.6%) and 36 (55.4%) patients received daily doses of 150mg and 300mg or more, respectively. However, due to missing demographic data from 8 patients, the final multivariate adjusted analysis was performed on a subset of 107 patients (Extended Data Figure 1). The baseline characteristics of the different cohorts and subgroups were summarized in Extended Data Table 1-2.

• Download figure
• Open in new tab

Extended Data Figure 1.

Study participant flowchart

It is worth noting that the patients included in this study were young individuals without any comorbidities, and all had received at least two doses of the COVID-19 vaccine. Although UDCA has been proposed as a potential treatment for COVID-19¹¹, the available data is limited to retrospective analyses on patients with cholestasis, progressive liver disease, cirrhosis, or liver transplantation, which are confounded by numerous effects related to the complexity of the target FXR. Specifically, gene mutations and differences in bile acid profiles among patients with different types of cholestasis can affect FXR expression, making it more vulnerable ^16–18. Additionally, studies have shown that the bile acid composition after UDCA treatment differs between healthy individuals and those with primary biliary cholangitis (PBC) ^19,20. Therefore, it is crucial to evaluate the effectiveness of UDCA in healthy individuals while controlling for comorbidities and eliminating confounding factors.

UDCA may promote the body temperature recovery of COVID-19

Among the 115 patients who were included in the univariate analysis (65 treated with UDCA and 50 without UDCA), 61 patients (93.8%) who received UDCA and 38 patients (76.0%) who did not receive UDCA had their body temperature returned to 37.3^°C or lower before discharge. Moreover, the median recovery time for the UDCA group was 4 days (95% CI, 4-4), compared to 5 days (95% CI, 4-5) for the no-UDCA group. The difference in the distribution of recovery speed between the two groups was significant (P = 0.030), as shown in Figure 1a.

• Download figure
• Open in new tab

Figure 1: Univariate and multivariate analyses of UDCA and no-UDCA groups.

a, Univariate analysis of UDCA and no-UDCA groups. b, Multivariate analyses of UDCA and no-UDCA groups with major adjusted, unadjusted and different minor adjusted.

• Download figure
• Open in new tab

Figure 2: Univariate and multivariate analyses of subgroups.

a, Univariate analysis of subgroups. b, Multivariate analyses of subgroups with major adjusted, unadjusted and different minor adjusted.

After adjusting for covariates, including age, sex, BMI, concomitant NSAID use, days from fever onset to hospitalization, and temperature on admission, the Cox proportional hazards models showed that UDCA treatment was associated with a 62% higher benefit on the speed of Body Temperature Recovery (HR 1.62, 95% CI 0.992.60) compared to the no-UDCA group (Figure 1b). In the unadjusted analysis, the hazard ratio of the UDCA group referring to the no-UDCA group was 1.87 (95% CI, 1.06-3.32). The minor adjusted results of the Cox proportional hazards regression model were consistent with the unadjusted analysis, even after different degrees of adjustments for various covariates, indicating the potential utility of UDCA in improving recovery time with a significant level of confidence. Furthermore, to examine the effect of missing covariates on the results, multivariate analyses with imputation were performed and shown in (Extended Data Table 3).

High dose performed better on primary endpoint

To provide further evidence for the efficacy of UDCA in treating COVID-19, we conducted a subgroup analysis based on dosage. Patients who received a daily dose of 300mg or more of UDCA showed a significant improvement in recovery time compared to the no-UDCA group, regardless of multiplicity adjustment. In contrast, patients who received a daily dose of 150mg did not show a significant difference in recovery time (Figure 1a).

The multivariate analyses were further conducted using Cox proportional models in two subgroups. Results revealed that the major adjusted hazard ratio for the 300mg or more subgroup was significantly associated with a shorter time to Body Temperature Recovery ([HR], 1.82 [95% CI, 1.07-3.10]) compared to the no-UDCA group (Figure 1b), while the 150mg subgroup did not demonstrate such an association without multiplicity adjusting. Additionally, minor adjustment methods yielded similar results to the primary analysis with major adjustments. Moreover, the estimated hazard ratios of both subgroups in the unadjusted analysis did not show significant differences from the baseline no-UDCA group, with multiplicity adjusting. Multivariate analyses with imputation were shown in (Extended Data Table 4), indicating that the high-dose subgroup had more significant results when the sample size was sufficient.

Overview of VirtualBody

VirutalBody comprises three integral components: the Human feAtures conditioned Diffusion (HAD) Model, the Skeleton Model, and the Skin Model. These three components underwent pre-training and fine-tuning using their distinct datasets, and their pipelines are shown in Supplementary Materials¹. Upon the completion of the training pipeline, VirtualBody had the capability to generate individualized PK data and inferred the dissemination of UDCA across multiple organ systems, even with extremely limited clinical data in the general population. Firstly, our approach leveraged high-dimensional human body characteristics that were automatically extracted from HAD Model. These characteristics were inputted into the Skeleton Model, which subsequently rolled out the venous drug concentration and drug concentrations in other organs under a specific dosage protocol autoregressively. Following this step, the Skin Model ultimately utilized the circulation data generated by the Skeleton Model to match the transition distribution of PK data, thus enabling the simulation of various clinical results.

Technically, VirtualBody models the transition of PK as a Markov Decision Process (MDP), defined by a tuple M = (S, A, M ^∗, r, γ, ρ₀) ²¹, where S is the state space, A is the action space, and M ^∗ is the optimal (i.e., real-world) transition model. As for our MDP, it mainly consisted of three key elements:

1. States incorporated detailed information about the drug, including its chemical properties and patientspecific features derived from the HAD Model. Time information is also included to ensure the Markov property is satisfied. In addition to the time-independent features, the dynamic drug concentration transition data is an essential dimension for the states, as our model is optimized to generate drug concentration-time curve.

2. Actions at a given time step were two-dimensional, with the first dimension being the medicine dosage and the second dimension representing patients’ diets.

3. State transitions defined the probability distribution over the next states, i.e. the PK data at the next time point, given a state and an action.

The learning objective of VirtualBody was to predict the subsequent state for each time step within the trajectory in this MDP, based on the corresponding state-action pair at that particular step.

Overall, VirtualBody enabled the handling of missing data during inference, the accommodation of various dosing strategies, and the prediction of venous and organ concentrations at any given time step. It is skilled in adapting to diverse dosage protocols and variations among individuals during inference, thus providing information to change the dosage regimen for the better.

VirtualBody generates accurate PK data for UDCA

The performance of VirtualBody was compared with the traditional top-down approach PopPK and other supervised learning methods, using longitudinal time-concentration data collected from a bioequivalence trial. The trial involved 120 individuals in a randomized, open-label, two periods, two sequences, crossover, single dose, fasting, and fasted design. The data was collected over a 72-hour period at various time points. The evaluation of PK-PD analysis was performed based on fundamental evaluation indices including Goodness of Fit (GOF), Area Under the Curve (AUC), and Maximum Concentration (C_max). It was demonstrated that VirtualBody surpassed other baselines in the analysis. The evaluation results are presented in Table.1.

Goodness of Fit (GOF) is a statistical measure that point-wisely evaluates the accuracy of a model in predicting observed drug concentration data. Based on the evaluation results, it is apparent that VirtualBody gained superiority over all other baselines on all metrics, achieving exceptional performance in modeling drug behavior in vivo. Specifically, R-square is a widely used indicator of how well a model fits the data, with a value close to one indicating a good fit. UDCA, being an endogenous bile acid with many uncertain PK properties, generally exhibited a negative baseline R-square. This can be attributed to factors such as accumulation in the bile acid pool, a typical hepatointestinal circulation profile, and susceptibility to the effects of diet and gallbladder contraction. However, a positive R-square (0.123) was achieved by VirtualBody in the prediction of UDCA PK data, making it the only model to our knowledge with this ability, and demonstrating its accuracy in individual variation prediction. Moreover, the outstanding performance of VirtualBody was further evidenced by its low Mean Absolute Error (MAE) and high correlation with the observed data, indicating its impressive ability not only to predict drug concentrations more accurately at each timestep but also to successfully capture the underlying trend of the real data.

In regards to one of the most crucial PK parameters, AUC_0-72h, VirtualBody outperformed other models with an MAE of only 60% and a MAPE of 50% or even smaller than other models. AUC is a crucial quantitative parameter for describing the total drug exposure over time in vivo, taking into account various factors in the ADME process. This parameter aided in establishing a more precise exposure-efficacy model, usually with less variability in comparison to the dose-efficacy relationship. With respect to C_max, it is noteworthy that although the MAPE of VirtualBody was slightly larger than that of PopPK, VirtualBody still achieved the lowest MAE. Remarkably, VirtualBody had the minimum standard deviation in almost all indicators, indicating that VirtualBody was able to maintain high robustness with accurate prediction.

These phenomena indicated that VirtualBody not only delivered more precise predictions of blood drug concentration but also offered better stability in its predictions. As for the conventional methodologies for PK/PD analysis, PopPK suffered from time-consuming human adjustments and its embedded nonlinear mixed-effects model had limited expressive capacity. Supervised learning methods, on the other hand, were likely to be data-hungry if we aimed to achieve strong generalization. Additionally, it can be observed that the most accurate inter-subject variation was captured by VirtualBody, as exhibited by its lowest uncertainty.

PK-PD analysis with VirtualBody

VirtualBody served as a crucial intermediary between pharmacokinetics and pharmacodynamics by precisely forecasting individual pharmacokinetic profiles, guaranteeing a comprehensive understanding of implicit knowledge in clinical data, and facilitating the discovery of novel insights. Logistic regression analysis was conducted to evaluate the impact of UDCA dose, predicted AUC and AAIC as continuous variables on the likelihood of body temperature recovery. These analyses were based on a longitudinal expansion dataset from original clinical data combined with predicted results, as shown in the schematic diagram of Figure 3a. The results showed that cumulative dose, AUC, and all of AAIC had statistically significant effects on the event of body temperature recovery, with odds ratios greater than 0, after the adjustment for the days from fever onset. However, the odds ratio of AUC was higher than that of dose, indicating a stronger correlation between AUC and febrile events (Figure 3b). Moreover, by increasing the minimum effective concentration of AAIC, there was a corresponding increase in its odds ratio which reached its peak at 1.5μg/mL. This observation highlighted the importance of excluding low-dose AUC and provided practical evidence for UDCA dosage recommendations in clinical administration. The histogram vividly showed that the AUC of patients with body temperature recovery was higher than that of patients without body temperature recovery on the same day after fever, whether in the complete cohort or in UDCA treated group (Figure 3c).

• Download figure
• Open in new tab

Figure 3: Data expansion and further analysis.

a, Depiction for longitudinal expansion for each individual. b, Logistic regression in longitudinal expansion data. c, Histogram for exposure-efficacy relationship in two groups, the AUC bar of corresponding days was shown as Mean±SD. d, Population separation according to AUC value by decision tree model from a virtual population of 1000 individuals

To illustrate the model’s expansion and identify drug suitability for specific populations, a decision tree was employed to evaluate the impact of various features. Firstly, the impact of age on UDCA exposure was found to be the most significant factor according to the root node in our decision tree. The model detected a substantial difference in average AUC between the two nodes under the age branches. Given UDCA’s high solubility and low permeability, it was more readily absorbed in the elderly due to the slowing of gastrointestinal peristalsis. This finding validated the model’s ability to learn the effects of age on drug absorption and confirmed the pharmacokinetic characteristics of UDCA. Therefore, appropriate dose reduction was recommended for elderly patients receiving UDCA treatment. Population differences were reflected in certain liver function indicators, such as AST and albumin levels. As liver function declined, the AUC of UDCA increased due to its reliance on liver metabolism and transport. This highlighted the sensitivity of UDCA exposure to changes in liver function, with elevated AST and decreased albumin levels resulting in higher AUC values (Figure 3d).

Study design and population

We accessed and reviewed the de-identified data through the Electronic Data Capture system of Nanjing University Hospital in December 2022. The Nanjing University Hospital institutional review boards approved this study as minimal risk and waived informed consent requirements. We excluded patients who experienced fever due to COVID-19 for four or more days ²⁶.

Exposure and endpoints variables

All patients had no previous history of UDCA medication before hospitalization. Patients in the UDCA group were initiated on UDCA treatment on the day of admission, while patients in the no UDCA group did not receive UDCA treatment throughout their hospitalization. We organized patient data of the cohort encompassing demographics (sex, age, BMI), the clinical baseline on admission (days from fever onset to hospitalization and temperature on admission), concomitant drug use (non-steroidal anti-inflammatory drugs, other symptomatic drugs recommended by Chinese National Diagnosis and Treatment of Novel Coronavirus Pneumonia (Version 9) for alleviating symptoms other than fever, such as sore throat, cough, shortness of breath, vomit, and diarrhea), and clinical body temperature recovery. Patients in the UDCA group were assigned to subgroups based on their daily dose, namely 150 mg and 300 mg or more, as only two individuals received 450 mg besides 150 mg and 300 mg.

Our primary endpoint was the time to Body Temperature Recovery, which was defined as the duration it took for body temperature to return to 37.3^°C or lower after admission and censored on the day of discharge without recovery. For patients with recurrent fever, we recorded the last time of Body Temperature Recovery. The final body temperature for analysis was the measured axillary temperature + 0.5^°C.

Preliminary statistical analysis

Baseline characteristics of patients were presented as the median and interquartile range (IQR), unless noted. The differences between UDCA and no UDCA groups in categorical and continuous variables were statistically examined using chi-square tests and Wilcox tests, respectively. All statistical comparisons were two-sided, with significance assumed at α = 0.05. · Univariate Kaplan-Meier was used to show recovery curves between the two cohorts and subgroups categorized by the daily dose of UDCA, using the primary endpoint as the evaluation measure without adjusting for confounding factors. The median recovery time of the two cohorts and subgroups was estimated with the Brookmeyer-Crowley method, with a 95 % confidence interval. The Log-rank test was applied to compare recovery curves with and without UDCA cohorts. Furthermore, the same method was adapted to compare two subgroups within the UDCA cohort based on their daily dose with the no-UDCA cohort.

To assess the impact of UDCA versus no-UDCA on the time to Body Temperature Recovery, the primary analysis was conducted using Cox proportional hazard model with all covariates adjusted. The covariates available for adjustment include demographical covariates (sex, age, and BMI), concomitant drug use (whether to use NSAIDs), and clinical baseline on admission (days from fever onset to hospitalization and temperature on admission). To avoid bias caused by the selection of covariates, we conducted a sensitivity analysis not adjusting for all covariates, minor adjusting for demographical features, concomitant drug use, and one of the clinical baselines on admission (days from fever onset to hospitalization or temperature on admission). The continuous covariates were normalized to the 0-1 range. The same method was applied to the evaluation between the two dose subgroups and the no-UDCA cohort.

To handle missing covariates, we employed the multiple imputations by chained equations method with ten imputed datasets and ten iterations. The missing covariates were allowed a maximal fraction of missing information (FMI)/m ≤ 10% (Supplementary Table S2). For each covariate, a random forest model was utilized. This method was adopted under the assumption that the missing observations for covariates were missing at random. The validity of the imputations was assessed by comparing the distribution of recorded and imputed values for all measurements through visual inspection. To further evaluate the robustness of our results, we performed a sensitivity analysis on both samples with complete covariates and those with imputed data. The imputation was performed using the mice package. All statistical analyses were performed using R (version 4.2.1). The P values of subgroup analyses were presented both before and after Benjamini & Hochberg adjustment.

The VirtualBody workflow

The model training pipeline was a three-stage framework that comprised feature alignment, VirtualBody pretraining, and VirtualBody finetuning. Firstly, in the feature alignment stage, noticing that there existed a significant gap between the features from PBPK data and the clinical data, we employed a state-of-the-art generative model to align the PBPK features and authentic features, and details were provided in the

Supplementary Information²

Secondly, in the pre-training stage, we separately trained each organ neural network component with precollected data, which was encoded into the neural network model by adopting a data-driven approach. We then integrated all the components into the systemic circulation, based on prior physiological mechanisms to form the Skeleton Model. Finally, the fine-tuning stage was to bias the model towards the experimental outputs by utilizing the authentic data and still retaining the pre-trained information in the meanwhile.

VirtualBody training and validation

There is no denying that the pharmacokinetic process is complex, highly varied among individuals, and influenced by multiple features. During the pre-training state, our goal was for VirtualBody to capture all of the critical prior knowledge.

The authentic UDCA data used for fine-tuning consisted of longitudinal time-concentration data collected from the bioequivalence trial initiated by Azpharm group limited in Guangdong, China (Registration Number: CTR20192232, http://www.chinadrugtrials.org.cn/), and the results of this trial has been reviewed by The Center for Drug Evaluation, National Medical Products Administration, China. The half of 120 subjects received single 500mg ursodeoxycholic acid tablets orally on fasted state, and the other half were administrated the same regimen on the fed state.

A 5-fold cross-validation was conducted on the authentic UDCA dataset to assess the prediction accuracy of the VirtualBody and other baselines. Our evaluation primarily focused on the Goodness of Fit, Area Under the Curve, and Maximum Concentration, to exhibit the comparison of several models from two levels, the fitting of points, and PK parameters.

Goodness of Fit (GOF)

It was assumed that the total number of drug-time curves was N, with each curve consisting of H sampling points. The observed concentration and predicted concentration of the j-th sampling point of the i-th curve were represented by and . Additionally, and were defined to represent the mean value of the observed concentration and predicted concentration of the i-th curve, respectively. The evaluation indicators of GOF were the R square, mean absolute error, and correlation, represented by Eqs. (1), (2), and (3), respectively.

Area Under the Curve (AUC_0-72h)

Under similar assumptions in GOF, we defined ao_i and ap_i to represent the observed AUC_0-72h and predicted AUC_0-72h of the i-th curve respectively. The evaluation indicators of AUC_0-72h included mean absolute error Eq.(4) and mean absolute percentage error Eq.(5).

Maximum of Concentration (C_max)

Under similar assumptions in GOF, we defined co_i and cp_i to represent the observed maximum concentration and predicted maximum concentration of the i-th curve respectively. The evaluation indicators of C_max) included mean absolute error Eq.(6) and mean absolute percentage error Eq.(7).

PK-PD analysis with VirtualBody

Using VirtualBody, we predicted the concentration of UDCA for each sample and calculated the Area Under the Curve (AUC) and the Area Above the Inhibitory Concentration (AAIC) at eight different concentrations ranging from 0.25μg/mL to 2μg/mL. The resulting data were then integrated into the original real-world dataset. In PK-PD analysis, since the baseline of hospitalization in the real-world dataset was misaligned, the time-series data for each individual was divided into multiple samples based on their total duration of fever from the onset, according to the number of days with a fever. Furthermore, each sample was labeled with clinical outcomes as the clinical endpoint to conduct PD analysis, where only the outcome of the non-censored last sample corresponding to temperature recovery was recorded as 1. All other samples and the censored last sample were labeled as 0, indicating that the temperature had not recovered. Logistic regression was performed on the above dataset to investigate the impact of cumulative drug dose or predicted AUC or predicted AAIC on the clinical outcome of temperature recovery. Fever duration, the use of non-steroidal anti-inflammatory drugs, and demographic characteristics of the patients were included as covariates for adjustment.

To gain deeper insights into the capabilities of VirtualBody, we generated a virtual population of 1,000 individuals with a wide distribution of ages ranging from 18 to 60 years old and heights ranging from 1.6m to 1.8m. The subjects were administered a fixed daily dosage of 500mg for 5 days, and VirtualBody was used to predict and calculate the AUC for each individual at the end of the fifth day. Finally, we incorporated a variety of relevant features, such as height, weight, age, sex, and various laboratory test values as input for a decision tree algorithm to predict the AUC for each individual.