A Systems Pharmacology Approach to Rivaroxaban: Physiologically Based Modeling of Pharmacokinetics and Coagulation Dynamics

1. Introduction

Direct oral anticoagulants (DOACs) have transformed the prevention and treatment of thromboembolic disorders, offering predictable pharmacokinetics, and fewer drug–food interactions. One of the most widely prescribed DOACs is rivaroxaban (brand name: Xarelto), a selective, direct inhibitor of activated coagulation factor X (Factor Xa). By binding directly to factor Xa, rivaroxaban inhibits both its free and clot-bound form, reducing thrombin generation and fibrin clot formation [1,2]. Rivaroxaban is administered orally, which increases its efficacy and patient compliance, therefore eliminating the need for constant monitoring that is necessary with vitamin K antagonists (VKAs). These characteristics make rivaroxaban particularly effective for managing conditions that may require lifelong treatment, such as stroke consequences and systemic embolism in patients with non-valvular atrial fibrillation. It also prevents recurrent deep vein thrombosis or pulmonary embolism, especially in patients with permanent risk factors or unprovoked venous thromboembolism [3,4].

Rivaroxaban is rapidly absorbed after oral administration and normally reaches peak plasma concentrations within 2 to 4 hours. It exhibits high bioavailability, between 80% to 100% for a 10 mg dose, independent of food intake. Interestingly, the patient’s fasting state plays a more important role when increasing the dose. In particular, to achieve the same bioavailability of 80% to 100% for a 15 mg and 20 mg dose, it is necessary to take rivaroxaban with food [5,6,7]. On the contrary, the bioavailability for a dose of 15 mg to 20 mg can drop to around 66% if taken while fasting [8]. Rivaroxaban is highly protein-bound (fraction unbound is equal to 5–8%) and has a moderate volume of distribution (around 50 L). It is eliminated via two different pathways: approximately one-third is excreted unchanged by the kidneys, while the rest is metabolized in the liver via CYP3A4/A5 and CYP2J2, and the inactive metabolites are subsequently excreted through biliary and renal routes [9].

Two key laboratory assays used to evaluate coagulation are prothrombin time (PT) and activated partial thromboplastin time (aPTT), which assess the intrinsic and extrinsic coagulation pathways, respectively [10]. Both PT and aPTT are prolonged by rivaroxaban in a dose-dependent manner [11,12,13,14], reflecting its anticoagulant effect, which follows an E_max model of Factor Xa inhibition [15].

Despite its predictable pharmacologic profile, significant inter-individual variability in drug exposure remains a clinical concern. Factors such as renal and hepatic functional state can, for instance, influence rivaroxaban pharmacokinetics and pharmacodynamics. In particular, in patients with moderate hepatic impairment (Child-Pugh B), total body clearance of rivaroxaban is decreased, leading to an increase in systemic exposure. This is accompanied by enhanced pharmacodynamic effects, such as greater inhibition of Factor Xa, raising the risk of bleeding complications [16]. Renal impairment also affects rivaroxaban pharmacokinetics and pharmacodynamics as the drug is partially eliminated by the kidneys [17,18]. Therefore, reduced renal function leads to a higher plasma concentration and increased bleeding risk. It is recommended to adjust the dose in patients with moderate renal impairment (creatinine clearance 30–49 mL/min), and rivaroxaban should be used with caution in these individuals, especially if they are taking other medications that may further increase rivaroxaban levels (i.e. potent CYP3A4 and/or P-gp inhibitors, such as ketoconazole) [19,20].

Understanding the influence of renal and hepatic function, age, weight, dose, and food interaction on rivaroxaban PK/PD is essential for optimizing therapy, especially in vulnerable populations. To better account for the variability caused by these factors, physiologically based pharmacokinetic/pharmacodynamic (PBPK/PD) modeling has emerged as a powerful approach. PBPK/PD models combine drug characteristics with patient-specific physiological data and can be employed to simulate how rivaroxaban behaves under different clinical conditions, providing insights for supporting individualized dosing strategies.

The aim of this work was to study rivaroxaban pharmacokinetics and pharmacodynamics using PBPK/PD modeling. By addressing key factors influencing intra- and inter-individual variability we aim to better understand how altered metabolism and excretion impact drug behaviour. Additionally, we will investigate the role of dosing regimens and fasting states in modulating rivaroxaban’s action. The outcomes of this research will contribute to the development of predictive models and support precision medicine approaches to optimize rivaroxaban therapy across diverse patient populations.

2. Materials and Methods

2.1. Systematic Literature Research

A systematic literature review was conducted for clinical pharmacokinetic and pharmacodynamic data on rivaroxaban (Figure 1). The keywords rivaroxaban AND pharmacokinetics were used to search the PubMed and PKPDAI databases. This search, conducted on February 3rd, 2025, yielded a wide range of studies, including single and multiple dose pharmacokinetic profiles, pharmacodynamic data, as well as population-specific data (renal and hepatic impairment, and food effect). To determine study eligibility, the following criteria were applied: Clinical trials reporting plasma concentration-time data, dose levels, and elimination parameters in adult humans were included. Preclinical studies, pediatric trials, or studies lacking primary pharmacokinetic endpoints were excluded.

2.2. Data Curation

Relevant pharmacokinetic and pharmacodynamic data were extracted from selected studies and curated in the PK-DB database [21]. Information on demographics (age, sex, weight), clinical condition (renal/hepatic status), rivaroxaban dose, co-administration with other medications, and fasting state was recorded. Data from figures were digitized using WebPlotDigitizer [22], and all numerical entries were standardized according to PK-DB’s data schema [21]. Pharmacokinetic parameters such as peak plasma concentration (C_max), time to reach peak plasma concentration (T_max), half-life (T_1/2), and area under the curve ((AUC_0–∞)) were included.

2.3. Model

A whole-body PBPK/PD model of rivaroxaban was developed in SBML format using sbmlutils [23], with all simulations performed using sbmlsim [24] and libRoadRunner [25,26]. Visualization was performed with cy3sbml [27]. The model structure includes absorption (oral input), distribution across major organ compartments, hepatic metabolism, renal and hepatic excretion, and coagulation dynamics.

Hepatic Impairment. Hepatic clearance was scaled using the parameter f_cirrhosis, where 0 represents normal function. Simulated impairments followed Child-Pugh classifications where 0.40 corresponded to mild, 0.70 to moderate and 0.81 to severe degree of cirrhosis [28].

Renal Impairment. Renal elimination was adjusted using KI__f_renal_function, based on estimated glomerular filtration rate. The model simulated normal (1.0), mild (0.5), moderate (0.35), and severe (0.2) renal impairment, consistent with KDIGO staging [29].

Food Effect. Food effect was controlled through GU__F_riv_abs, which represents the fraction of drug absorbed in the gastrointestinal tract. This parameter was scaled to account for the impact of food on rivaroxaban’s bioavailability. The model assumes a 18% reduction in absorption under fasted conditions (0.82) compared to food (1.0), consistent with clinical observations of reduced bioavailability when the drug is taken without food [5].

Coagulation Dynamics. The pharmacodynamic model describes rivaroxaban’s anticoagulant effects by simulating its inhibition of Factor Xa. This interaction is represented using an E_max model, with rivaroxaban venous plasma concentration (C_{ve_riv}) as input and fractional Factor Xa activity as output. The resulting inhibition cascade affects downstream coagulation markers, including prothrombin time (PT) and activated partial thromboplastin time (aPTT). These are modeled using similar E_max relationships, capturing both the absolute clotting times and their fold-changes relative to baseline values ($P{T}_{\mathrm{ref}}$ and $aPT{T}_{\mathrm{ref}}$).

$$\begin{array}{cc}\hfill \mathrm{Xa}\_\mathrm{inhibition}& =E_{\max ,\mathrm{Xa}}·{\displaystyle \frac{C_{\mathrm{ve}\_\mathrm{riv}}}{C_{\mathrm{ve}\_\mathrm{riv}}+E{C}_{50,\mathrm{Xa}}}}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \mathrm{PT}& =P{T}_{\mathrm{ref}}\left(1+E_{\max ,\mathrm{PT}}·{\displaystyle \frac{C_{\mathrm{ve}\_\mathrm{riv}}}{C_{\mathrm{ve}\_\mathrm{riv}}+E{C}_{50,\mathrm{PT}}}}\right)\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \mathrm{aPTT}& =aPT{T}_{\mathrm{ref}}\left(1+E_{\max ,\mathrm{aPTT}}·{\displaystyle \frac{C_{\mathrm{ve}\_\mathrm{riv}}}{C_{\mathrm{ve}\_\mathrm{riv}}+E{C}_{50,\mathrm{aPTT}}}}\right)\hfill \end{array}$$

(3)

2.4. Parameter Optimization

Both PK and PD model parameters were optimized using a subset of the curated clinical data. The fitted parameters of the PK model consisted of gastrointestinal dissolution rate (GU__Ka_dis_riv), absorption rate (GU__RIVABS_k), hepatic metabolic rate (LI__RIV2RX_Vmax), and renal elimination rate (KI__RIVEX_k). For the PD model, E_max and EC₅₀ values for Factor Xa inhibition (Emax_Xa, EC50_riv_Xa), prothrombin time (Emax_PT, EC50_riv_PT), and activated partial thromboplastin time (Emax_aPTT, EC50_riv_aPTT) were optimized. Optimization was performed across multiple dose levels and under both fed and fasted conditions. Multiple local optimization runs ($n=100$) were conducted to first optimize the PK parameters and subsequently the PD parameters. Model performance was evaluated using independent datasets not used for calibration.

3. Results

3.1. Rivaroxaban Database

An open database of pharmacokinetic and pharmacodynamic data from 14 clinical studies was established, including diverse populations, dosing regimens and physiological conditions (Table 1). The data was used to develop the rivaroxaban PBPK/PD model.

3.2. Computational Model

A whole-body PBPK/PD model of rivaroxaban was created that includes key factors that determine intra- and inter-individual variability (Figure 2). The model includes the organs involved in the pharmacokinetics of rivaroxaban (gastrointestinal tract, liver, and kidney), connected through the systemic circulation, as well as the main pharmacodynamic outputs (activated Factor Xa inhibition, prothrombin time, and activated partial thromboplastin time). The main patient-specific factors that are known to affect the pharmacokinetics and pharmacodynamics of rivaroxaban are included in the model. Specifically, four key factors were included: dose dependency, liver impairment, renal impairment, and food effect.

3.3. Dose Dependency

The model confirmed the dose dependency of rivaroxaban pharmacokinetic and pharmacodynamic profiles in a wide range of oral doses (1-100 mg), as shown in Figure 3. Simulated plasma concentration-time profiles matched observed data from multiple studies, both with single and multiple doses, demonstrating increases in C_max and (AUC_0–∞) with increasing doses. More specifically, the (AUC_0–∞) increased proportionally from 1 to 20 μmolh/L as the dose increased from 1 to 100 mg. The sum of rivaroxaban metabolites also exhibited a proportional, if not as sharp, increase, going up to about 4 μmol/L when the dose reached 100 mg. The half-life remained constant across all doses. The simulations matched the data from Graff2007 [14], Jiang2010 [12], Kubitza2005a [11] and Zhao2009 [13]. The urinary data from Jiang2010 [12] and Zhao2009 [13] were also consistently in good agreement with the simulations, which exhibited higher levels of rivaroxaban in the urine at higher doses. Additionally, pharmacodynamic responses are also scaled with the dose, showing higher values of maximum prothrombin time (up to 34 $\mathrm{s}$ for a 100 mg dose), activated partial thromboplastin time (up to 43 $\mathrm{s}$), and activated factor X inhibition (up to 65%) for higher doses.

3.4. Food Effect

The model reproduced the known effects of food on rivaroxaban pharmacokinetics and pharmacodynamics (Figure 4), with enhanced oral bioavailability and increased coagulation effects when the drug is co-administered with food. Following administration with food, and the subsequent increase in the fraction of the drug that is absorbed (from about 0.82 to 1.0), the (AUC_0–∞)_0–∞ of both rivaroxaban and its metabolites increased by approximately 20% compared to the fasted state. Despite the changes in exposure, the half-life of rivaroxaban remained relatively unchanged, indicating that the observed increases were primarily driven by increased bioavailability, rather than alterations in clearance. Pharmacodynamically, food intake amplified rivaroxaban’s anticoagulant effects. The maximum prothrombin time increased from 15.4 to 15.9 seconds, and maximum aPTT rose from 38.2 to 38.9 seconds. The peak Factor Xa inhibition increased from 0.30 to 0.34, reflecting the higher systemic exposure to active drug under fed conditions. These findings align with the clinical data from Kubitza2006 [6] and Liu2024 [7], as well as with label recommendations to administer rivaroxaban with food, especially for doses greater than 15 mg.

3.5. Hepatic Impairment

Simulation of rivaroxaban pharmacokinetics and pharmacodynamics under increasing degrees of hepatic impairment showed a consistent increase in rivaroxaban exposure (Figure 5). (AUC_0–∞) values increased with worsening liver function from approximately 2 to 2.5 μmolh/L. In contrast, rivaroxaban metabolites displayed a decreasing trend, with (AUC_0–∞) falling from 0.5 to 0 μmolh/L, which is consistent with the reduced hepatic metabolic activity in the case of cirrhosis. The simulated half-life only increased minimally with the progression of liver cirrhosis.

On the pharmacodynamic side, prothrombin time (PT) increased from 15.5 to 18.5 $\mathrm{s}$, activated partial thromboplastin time (aPTT) from 39 to 44 $\mathrm{s}$, and Factor Xa inhibition from 0.32 to 0.44 across the hepatic impairment scale. These results align with the higher exposure that leads to increased anticoagulation effects.

3.6. Renal Impairment

Figure 6 shows the effect of renal function on the pharmacokinetics and pharmacodynamics of rivaroxaban. With decreasing renal function, systemic exposure to rivaroxaban and its metabolites increases. Specifically, the model predicts an increase in rivaroxaban (AUC_0–∞) from 1 to 3 μmolh/L, and in metabolite (AUC_0–∞) from 0 to 5 μmolh/L.

These trends reflect enhanced drug accumulation associated with worse renal clearance. The minor increase of half-life (from 4.35 to 4.36 $\mathrm{h}$ for rivaroxaban and from 4.35 to 4.42 $\mathrm{h}$ for metabolites) also supports a marginal decrease in the drug excretion. All this is aligned with the underlying physiology and with clinical findings from Greenblatt2018 [18] and Moore2014 [17], which show greater exposure with worsening renal function.

4. Discussion

In this study, we developed and evaluated a whole-body physiologically based pharmacokinetic/pharmacodynamic (PBPK/PD) model of rivaroxaban, a direct oral anticoagulant commonly prescribed for the prevention and treatment of thromboembolic events. The model was constructed using a comprehensive, literature-derived database of clinical pharmacokinetic and pharmacodynamic data, enabling simulation of plasma concentrations, drug exposure metrics ((AUC_0–∞), C_max, half-life), and PD characteristics (Factor Xa inhibition, PT, aPTT) across a range of physiological conditions. This work demonstrates the potential of PBPK/PD modeling to support precision dosing and address questions that could arise in real-world clinical scenarios, even in the absence of complete empirical datasets.

The curated database used in this project not only allowed robust model calibration and validation but also laid the foundation for future clinical decision support systems (CDSS) and digital twin applications. By aggregating and harmonizing clinical trial data, including various dosing regimens, patient populations, and intervention contexts, the database allows flexible reuse for both academic and clinical purposes. Importantly, it serves as a scalable resource that can be linked to other pharmacological and physiological models in systems pharmacology frameworks.

Our dose dependency simulations successfully reproduced the expected increases in rivaroxaban pharmacokinetic and pharmacodynamic response across a range of doses. These results are consistent with clinical findings and support the model’s capacity to simulate rivaroxaban pharmacology under varying therapeutic scenarios.

In terms of food effects, the model correctly captured the qualitative trend of increased rivaroxaban exposure when administered with food, in alignment with known clinical data. Further improvements could include making the food effect itself dose-dependent. Notably, the fact that food has a greater impact on rivaroxaban bioavailability at higher doses (particularly above 15 mg) [34] is an effect that was not fully captured. Incorporating this knowledge, e.g. by modifying intestinal permeability or gastric emptying rates in a dose-dependent manner, could enhance the model’s translational relevance, especially for clinical scenarios that require the administration of high doses of rivaroxaban.

As expected, simulations of hepatic impairment demonstrated increases in rivaroxaban exposure and pharmacodynamic effect with progression of liver dysfunction. This is consistent with the known role of hepatic metabolism and hepatic excretion in rivaroxaban clearance, and with clinical observations [16]. The model also predicted a progressive decline in the formation of rivaroxaban metabolites, which is biologically plausible given the reduced metabolic capacity in cirrhosis. Additionally, Factor Xa inhibition, PT, and aPTT all increased correspondingly, further validating the model’s ability to accurately reproduce liver impairment.

Renal impairment simulations likewise reflected physiologically consistent trends. Exposure metrics such as (AUC_0–∞) increased with worsening renal function, reflecting increased drug accumulation due to reduced renal clearance. While the model slightly underestimated the extent of drug retention seen in severe renal impairment in some studies [17], it was nonetheless able to replicate concentration–time profiles observed in plasma and urine. In particular, although the detailed pharmacodynamic timecourses under renal impairment conditions are lacking, the data reported by Moore2014 showed increases in (AUC_0–∞) for the PD curve and in E_max in individuals with impaired renal function. The model accurately represented this trends by predicting higher values in Factor Xa inhibition, prolonged PT, and elevated aPTT when the renal function is impaired, which is consistent with higher systemic exposure that would be expected from impaired clearance. This demonstrates one of the major strengths of PBPK/PD modeling: the ability to bridge data gaps and make physiologically grounded predictions even where detailed empirical data are lacking.

Looking forward, the model offers numerous opportunities for expansion. Demographic and physiological factors such as age, sex, and body weight have well-documented effects on rivaroxaban pharmacokinetics [30,35] and could be incorporated to support population-level simulations or individualized dosing. In particular, future model versions could differentiate between lean body mass and adipose tissue in order to more accurately capture distribution volumes and clearance dynamics in obese or underweight individuals. This distinction could be especially important given the moderate lipophilicity of rivaroxaban [36] and the increasing clinical emphasis on dosing optimization in diverse body types.

Additional refinements may include better representation of transporter proteins (e.g., P-gp) and enzyme polymorphisms (e.g., CYP3A4/5), as well as simulation of drug–drug interactions and adherence scenarios. Ultimately, the integration of this model into digital twin platforms could enable real-time therapeutic monitoring and decision-making support, particularly for vulnerable populations such as the elderly or patients with organ dysfunction.

5. Conclusions

In this study, we developed a whole-body physiologically based pharmacokinetic/pharmacodynamic (PBPK/PD) model of rivaroxaban to simulate its behaviour under varying physiological conditions and treatment scenarios. Using the established clinical database, we successfully reproduced key pharmacokinetic and pharmacodynamic outcomes across different doses, food intake states, and levels of hepatic and renal function. This work demonstrates the value of PBPK/PD modeling in understanding complex drug behaviour, supporting personalized dosing strategies, and informing clinical decision-making. The database and model together form a strong foundation for future integration into digital twin platforms and clinical decision support systems.