An Open and Reproducible Digital Twin for Personalised Anticoagulant Therapy: Modelling Apixaban Pharmacokinetics and Pharmacodynamics

1. Introduction

Thrombotic events, such as deep vein thrombosis and pulmonary embolism, are life-threatening conditions that substantially contribute to morbidity and mortality. Although comprehensive epidemiological data remain limited, available studies consistently show that these events become more frequent with increasing age [1]. With population ageing, the clinical need for long-term anticoagulant therapies with reliable pharmacokinetic (PK) and pharmacodynamic (PD) profiles and improved safety is increasing [2,3].

Apixaban, an oral direct factor Xa inhibitor, addresses several of these needs. It has a favourable safety profile and is associated with a lower risk of haemorrhage than traditional anticoagulants [4,5]. Its PK and PD profiles are generally stable and predictable. Together with its limited dependence on fasting state, these characteristics make apixaban attractive for long-term outpatient management [6]. However, apixaban exposure and response may still be altered by patient-specific factors, including body weight, co-administered medications, genetic polymorphisms, and renal and hepatic function [7,8,9,10,11]. Translating these complex, multifactorial influences into individualised, quantitative dose optimisation remains challenging. This can increase the risk of both overdosing and underdosing, with clinically relevant consequences: excessive dosing increases bleeding risk, whereas underdosing may fail to prevent thrombotic events [3]. Although careful monitoring can support dose adjustment in controlled clinical settings, such resources are often limited in outpatient care. Transparent and practical approaches are therefore needed to optimise apixaban dosing for individual patients in real-world ambulatory settings.

Physiologically based pharmacokinetic/pharmacodynamic (PBPK/PD) modelling provides a potential solution. PBPK/PD models describe the body as a system of interconnected compartments, each defined by organ-specific physiological and biochemical properties. By integrating relevant physiological, biochemical, and drug-specific processes, these models enable quantitative analysis of the relationship between drug exposure and therapeutic response, as well as the evaluation of different dosing regimens and routes of administration. Moreover, PBPK/PD models can be individualised using patient-specific parameters, such as age, body weight, sex, and organ function. They therefore provide a mechanistic basis for personalised dosing predictions and represent an important step towards digital twins: personalised computational models capable of simulating individual patient physiology and treatment response [12,13].

Despite their promise, several factors limit the predictive performance of PBPK/PD models and their translation into clinical practice. Models are often parameterised using data from a limited number of sources, which may not adequately represent the relevant target populations and can introduce undercoverage bias. Restricted access to the underlying datasets further limits independent validation, reuse, and model refinement. In addition, the use of commercial modelling environments with closed-source implementations reduces transparency, while proprietary executable formats hinder interoperability. Together, these factors can result in isolated models built on closed data, code and software, limiting both the FAIRness–Findability, Accessibility, Interoperability, and Reusability–and the reproducibility of PBPK/PD modelling research [14,15,16,17].

To our knowledge, published physiologically based models of apixaban show similar limitations, although not every model exhibits all of the issues described above. In each case, at least one limitation related to data coverage, transparency, interoperability, or reproducibility remains. Moreover, most models focus primarily on apixaban pharmacokinetics and do not explicitly represent apixaban metabolites, pharmacodynamic components, or both. When such components are included, their scope is often limited. A systematic overview of existing models in adult populations, including their scope, software implementation, availability, and reproducibility, is provided in Supplementary Materials Section S2 [18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39].

To address these limitations, this study: (i) performs a systematic literature review of available apixaban PK and PD studies and corresponding published computational models; (ii) builds a comprehensive, curated dataset of apixaban pharmacokinetic and pharmacodynamic data and makes it openly available through the PK-DB database (https://pk-db.com) [40]; (iii) develops an expanded PBPK/PD digital twin of apixaban, with separate representation of metabolites and enhanced pharmacodynamic outputs; (iv) implements the model in a modular structure using the open standard format SBML, with detailed metadata linked to established ontologies [41,42,43]; (v) parameterises and evaluates the model using the curated dataset; (vi) quantitatively investigates key determinants of apixaban PK and PD, including dose, renal and hepatic impairment, food intake, and body weight; and (vii) provides an open public repository containing the model code, simulation workflows, curated data, documentation, and main results.

2. Materials and Methods

2.1. Systematic Literature Research and Data Curation

A systematic literature search was conducted to identify studies reporting time-resolved pharmacokinetic and/or pharmacodynamic data for apixaban. PubMed was queried on 2025-10-24 using the search terms “pharmacokinetics AND apixaban”, and the PKPDAI database was screened in parallel for “apixaban” [44]. Eligible studies included adult human studies in healthy volunteers and patient populations, including food-effect and body-weight studies, drug–drug interaction studies when data for apixaban administered alone were available, co-administration studies with other direct oral anticoagulants, microdose cocktail studies, and studies in special populations with renal or hepatic impairment. Animal studies, paediatric clinical studies, and reports without time-resolved PK or PD data were excluded. In vitro studies were analysed separately to identify initial parameter values for subsequent optimisation. An overview of the study selection process is provided in Supplementary Figure S1.

Data from eligible studies were curated in the open pharmacokinetics database PK-DB [40], following established protocols. Extracted information included population characteristics, such as demographics, health status, fasting or fed state, and genetic information; clinical protocols, including formulation and dosing regimens; plasma concentration–time profiles of apixaban and its metabolites; recovery and cumulative amounts excreted in urine and faeces; and reported coagulation parameters, including prothrombin time (PT), modified prothrombin time (mPT), international normalised ratio (INR), activated partial thromboplastin time (aPTT), and anti-activated factor X activity (anti-Xa activity). Figure-based data were digitised using WebPlotDigitizer [45], whereas tabular and textual data were reformatted into standardised PK-DB formats. The complete curated dataset is publicly available through PK-DB and is also included with the model files; an overview of the included studies is provided in Table 1. To ensure curation quality, all digitised data were independently checked by at least two team members.

2.2. Computational Model

The complete model, including simulation scripts and documentation, is available in SBML format under a CC-BY 4.0 license via GitHub (https://github.com/matthiaskoenig/apixaban-model) and is archived on Zenodo (v0.6.0) [46]. Mathematical descriptions of the submodels and the corresponding ordinary differential equations are provided in the Supplementary Materials Section S4.

2.2.1. Model Structure

The PBPK/PD model was developed in the Systems Biology Markup Language (SBML) [42,43]. Programmatic model construction and visualisation were performed using the sbmlutils [47] and cy3sbml [48,49] libraries. Numerical solutions of the underlying ordinary differential equations (ODEs) were obtained with sbmlsim [50], using the high-performance SBML simulation engine libRoadRunner [51,52].

Fractional organ volumes and blood flows were taken from the literature [53]. Fractional compartment volumes were set to $F{V}_{\mathrm{gu}}=1.71\%$ for the gut, $F{V}_{\mathrm{ki}}=0.44\%$ for the kidneys, $F{V}_{\mathrm{li}}=2.10\%$ for the liver, and $F{V}_{\mathrm{lu}}=0.76\%$ for the lungs. Fractional blood flows were defined as $F{Q}_{\mathrm{gu}}=18.00\%$ for the gut, $F{Q}_{\mathrm{ki}}=19.00\%$ for the kidneys, $F{Q}_{\mathrm{h}}=21.50\%$ for hepatic venous outflow, and $F{Q}_{\mathrm{lu}}=100\%$ for the lungs. Absolute organ volumes and blood flows were calculated by scaling the corresponding fractional values by body weight.

Parameters describing renal and hepatic function, as well as apixaban absorption, were implemented as scaling factors because of their relevance for intra- and inter-individual variability in apixaban PK and PD:

• Renal impairment was modelled as a progressive reduction in renal function. The dimensionless factor $f_{\mathrm{renal}}$ was used to scale renal clearance of the corresponding substances. Scaling values were derived from the KDIGO Clinical Practice Guideline for the Evaluation and Management of Chronic Kidney Disease [54,55]. These values were based on the estimated glomerular filtration rate (eGFR, mL/min/ $1.73$ m2) defining each kidney-function category: normal kidney function (eGFR $\ge 90$, $f_{\mathrm{renal}}=1.00$); mild impairment ($60\le$ eGFR $\le 89$, $f_{\mathrm{renal}}=0.69$); moderate impairment ($30\le$ eGFR $\le 59$, $f_{\mathrm{renal}}=0.32$); severe impairment ($15\le$ eGFR $\le 29$, $f_{\mathrm{renal}}=0.24$); and end-stage kidney disease (eGFR $<15$, $f_{\mathrm{renal}}=0.10$). Dialysis was not included in the model.
• Hepatic impairment was modelled as a progressive reduction in functional liver tissue and shunting. The dimensionless factor $f_{\mathrm{cirrhosis}}$ represents the increase in non-functional liver parenchyma and the development of portosystemic collaterals. Values were assigned according to Child–Turcotte–Pugh (CTP) classes, which are used to predict mortality in patients with cirrhosis: class A (mild impairment, 5–6 points, $f_{\mathrm{cirrhosis}}=0.40$) and class B (moderate impairment, 7–9 points, $f_{\mathrm{cirrhosis}}=0.70$) [56,57,58,59].
• Food effect was modelled by scaling the fraction absorbed, thereby modulating intestinal apixaban absorption. The baseline fraction absorbed was set to $F_{\mathrm{absorption},\mathrm{api}}=0.66$ [60] and scaled by the dimensionless parameter $f_{\mathrm{absorption}}$. Fasted and not-reported conditions were assigned $f_{\mathrm{absorption}}=1.0$. The value of $f_{\mathrm{absorption}}$ under fed conditions was fitted using the corresponding datasets [61,62,63].

2.2.2. Apixaban Framework

In the model, apixaban can be administered orally as a tablet or solution, or intravenously as a solution (Figure 1B–E). After oral tablet administration, the tablet dissolves and is transported in the intestine, where apixaban is either absorbed into the portal vein or excreted in faeces. Intestinal absorption is represented by Michaelis–Menten kinetics with saturation at higher doses [64].

Apixaban disposition involves several enzymes and transport processes. In hepatocytes, apixaban is metabolised by CYP3A4/5, producing the metabolites M2, M4, and M7. M2 is subsequently conjugated by sulfotransferase SULT1A1, resulting in the formation of M1. Apixaban and the metabolites M1 and M7 are excreted into the urinary filtrate by the kidneys. Apixaban and its metabolites can also be excreted in faeces after biliary transport into the intestinal lumen.

The inhibition of activated factor X by apixaban is represented by corresponding changes in PT, mPT, INR, aPTT, and anti-Xa activity. Apixaban metabolites were assumed to have no pharmacological activity.

2.2.3. Model Assumptions

Additional details on the model equations and assumptions are provided in Supplementary Material S4. The key model assumptions and simplifications are summarised below.

• Apixaban absorption was modelled as a process following Michaelis-Menten kinetics without any apixaban efflux back to the intestinal lumen;
• Apixaban absorption was represented by Michaelis–Menten kinetics. Efflux of apixaban from the enterocytes back into the intestinal lumen was not included.
• Several apixaban metabolites have been reported: M1, M2, M4, M7, M10, M13 [6,65,66]. However, after oral administration of apixaban, only M1 and M7 were detected in substantial amounts in urine, and only M1, M2, and M7 were detected in faeces, with recovery $\ge 1\%$ in both studied groups [65]. Consequently, only M1, M2, and M7 were included in the model.
• Apixaban metabolism and metabolite conversion were modelled exclusively in the liver as a series of irreversible reactions following mass-action kinetics, assuming rapid equilibration for reversible transport processes. Conversion of apixaban to M2 and M7 was modelled following mass action kinetics, whereas conversion of M2 to M1 followed Michaelis–Menten kinetics. These metabolic reactions, as well as the export of apixaban and metabolites from hepatocytes, were represented explicitly, but without separate mechanistic descriptions of CYP3A4/5, SULT1A1, P-glycoprotein (ABCB1), or breast cancer resistance protein (ABCG2) activity. This simplification was used because, to our knowledge, there is no strong evidence that genetic polymorphisms in these genes have a clinically relevant effect on apixaban pharmacokinetics [67,68,69].
• Pharmacodynamic outputs were modelled through the direct dependency on apixaban concentration, following linear (anti-Xa activity) or Michaelis-Menten (other PD variables) kinetics. Study-specific initial values for PT, aPTT, and mPT were taken if reported. Otherwise, default values $12.5$ s, $28.4$ s, and $53.4$ s were used [60];
• The international normalised ratio (INR) was modelled independently from prothrombin time because the international sensitivity index varies between laboratory instruments and manufacturers.
• Anti-Xa activity in the studies is reported in two units: ng/mL and international units/mL. In the model, it is represented by two independent pharmacodynamic outputs, as translating international units to mass requires an apixaban-specific conversion factor, which is not standardised, depends on the assay parameters, and is not reported in the articles;
• Anti-Xa activity was reported in the curated studies using two different units: ng/mL and international units/mL. In the model, these were represented as two independent pharmacodynamic outputs because conversion from international units to mass concentration requires an apixaban-specific conversion factor. This factor is not standardised, depends on assay conditions, and was not reported in the articles.
• For mild and moderate hepatic impairment, no abnormalities in the coagulation pathway were assumed. Severe hepatic impairment was not included in the model.

2.3. Simulation Methodology

For each curated clinical study (Table 1), in silico experiments were performed to reproduce the reported dosing regimen, study conditions, and subject characteristics. Intravenous, oral tablet, and oral solution doses were specified using the parameters ${\mathrm{IVDOSE}}_{\mathrm{api}}$, ${\mathrm{PODOSE}}_{\mathrm{api}}$, and ${\mathrm{SOLDOSE}}_{\mathrm{api}}$, respectively. The simulation duration was selected to match the design of the corresponding clinical study. When reported, study- and subject-specific characteristics were adjusted, including:

• Body weight and height, which were used to scale organ and tissue volumes as well as blood flows;
• Glomerular filtration rate, which was used to account for renal function;
• Child–Turcotte–Pugh class, which was used to account for hepatic function;
• Baseline prothrombin time, modified prothrombin time, and activated partial thromboplastin time, which were used to represent the initial state of the coagulation pathway.

Multiple-dose regimens were implemented by stepwise numerical integration between dosing intervals, with dosing events applied according to the study-specific protocol.

2.4. Parameter Optimization

Parameter optimization minimized a cost function $F\left(\overrightarrow{p}\right)$ defined as the sum of squared weighted residuals $r_{i,k}$ across all timecourses k and data points i:

$$F\left(\overrightarrow{p}\right)=0.5\sum _{i,k}{\left(w_{i,k}·r_{i,k}\left(\overrightarrow{p}\right)\right)}^2,$$

where weights $w_{i,k}=n_k/{\sigma }_{i,k}$ combine weighting by study sample size $n_k$ and inverse measurement uncertainty ${\sigma }_{i,k}$, such that data points with lower standard errors from larger sample sizes contribute more to the objective function.

A subset of curated data from fasted healthy subjects following single- and multi-dose administration without co-administration was used to optimise parameters. This subset included concentration-time profiles of apixaban and its metabolites, their recovery in urine and faeces, and changes in time of pharmacodynamic outputs. Pharmacokinetic timecourses of apixaban measured in fed healthy subjects were used to optimise the parameter, scaling the apixaban fraction absorbed separately. Pharmacokinetic and pharmacodynamic parameters, such as maximal concentration, area under the curve (AUC), and maximal relative change from baseline, as well as digitised scatter plots data, special populations, and fed pharmacodynamic data, were not used in parameter optimisation.

Figure 1. Whole-body PBPK/PD model of apixaban and key processes influencing its exposure and effects. A) Whole-body model showing apixaban (API) administration, systemic circulation, and key organs (liver, kidney, GI tract) involved in absorption, metabolism, distribution, and excretion of API and its metabolites M1, M7. B) Intestinal model illustrating API tablet dissolution, absorption by enterocytes into the portal blood, and faecal excretion of API and its metabolites M1, M2, M7 that are excreted into the intestinal lumen with bile. C) Renal model showing uptake and urinary excretion of API and its metabolites M1 and M7. D) Hepatic model showing API uptake into hepatocytes, its intracellular conversion to M7 and M2 via CYP3A4/5 enzymes, subsequent conversion of M2 to M1 by SULTA1, release of M1 and M7 to the systemic circulation, and excretion of M1, M2, and M7 into the intestinal lumen in bile. E) Coagulation model illustrating API pharmacodynamics: inhibition of activated factor X (Xa) leading to inhibition of blood clotting. The extent of API-induced inhibition of blood clot formation is reflected by quantitative changes in specific coagulation parameters: Xa activity, anti-Xa activity, international normalised ratio (INR), prothrombin time (PT), modified prothrombin time (mPT), and activated partial thromboplastin time (aPTT).

Figure 1. Whole-body PBPK/PD model of apixaban and key processes influencing its exposure and effects. A) Whole-body model showing apixaban (API) administration, systemic circulation, and key organs (liver, kidney, GI tract) involved in absorption, metabolism, distribution, and excretion of API and its metabolites M1, M7. B) Intestinal model illustrating API tablet dissolution, absorption by enterocytes into the portal blood, and faecal excretion of API and its metabolites M1, M2, M7 that are excreted into the intestinal lumen with bile. C) Renal model showing uptake and urinary excretion of API and its metabolites M1 and M7. D) Hepatic model showing API uptake into hepatocytes, its intracellular conversion to M7 and M2 via CYP3A4/5 enzymes, subsequent conversion of M2 to M1 by SULTA1, release of M1 and M7 to the systemic circulation, and excretion of M1, M2, and M7 into the intestinal lumen in bile. E) Coagulation model illustrating API pharmacodynamics: inhibition of activated factor X (Xa) leading to inhibition of blood clotting. The extent of API-induced inhibition of blood clot formation is reflected by quantitative changes in specific coagulation parameters: Xa activity, anti-Xa activity, international normalised ratio (INR), prothrombin time (PT), modified prothrombin time (mPT), and activated partial thromboplastin time (aPTT).

Optimisation was performed using a local optimiser and followed a stepwise approach: at first, all parameters influencing apixaban were optimised on apixaban pharmacokinetic timecourses. The resulting absorption parameter values were fixed, and all other PK parameters were fitted. Parameter, scaling apixaban fraction absorbed, was optimised separately, followed by pharmacodynamic parameters. To improve robustness against local minima, multiple optimisation runs ($n=100$) were conducted using different initial parameter values. The global parameter set was applied across all in silico experiments without further study-specific parameter optimisation.

Parameter optimisation results, including fitting conditions, optimised pharmacokinetic and pharmacodynamic parameter sets, convergence behaviour plot, and goodness-of-fit assessments, are reported in Supplementary Materials Section S5.

2.5. Parameter Scans and Model Evaluation

To systematically investigate the influence of body weight, dose, kidney and liver function, and changes in absorption on the pharmacokinetics and pharmacodynamics of apixaban, parameter scans were conducted across physiologically and/or clinically relevant ranges of these variables: body weight from 38 kg to 200 kg (linear scale, 10 values), dose from 0.5 mg to 100 mg (linear scale, 10 values), $f_{\mathrm{renal}}$ from −1 to 1 (log scale, 10 values), $f_{\mathrm{cirrhosis}}$ from 0 to 0.7 (linear scale, 10 values), $f_{\mathrm{absorption}}$ from −1 to 1.5 (log scale, 10 values).

To calculate pharmacokinetic parameters from each scan across simulated concentration-time profiles of apixaban and its urinary recovery, the standard noncompartmental approach was used. Pharmacodynamic parameters (maximal change) were calculated using the same method, but across relevant simulations. The resulting set of parameters, time-course pharmacokinetics, and pharmacodynamics data for special populations, pharmacodynamics data on the fed population, and data digitised from scatter plots were compared against curated, corresponding clinical data.

First-order apixaban elimination kinetics and the linear dependence of AUC, maximum concentration, and maximal relative changes in PD parameters on dose across ranges studied in experiments enabled us to scale these parameters by the ratio of simulated to administered dose in the study, where needed, for proper comparison.

3. Results

3.1. Apixaban Database

Initially, 142 records were identified through PKPDAI and 743 records through PubMed. During the literature review, 12 additional relevant studies were identified manually and included in the screening. All records were assessed against the inclusion criteria described in the Materials and Methods section. The complete study selection process is shown in Supplementary Materials Figure S1.

In total, 35 studies met the predefined inclusion criteria. Pharmacokinetic and pharmacodynamic data from these studies were digitised, standardised, and independently checked by at least two team members. An overview of the curated study dataset is provided in Table 1, and all curated data are available through the open-access PK-DB database.

In addition, 22 published computational models describing apixaban pharmacokinetics in adults were identified. Among these, 4 included pharmacodynamic outputs and 10 used a PBPK modelling framework. A detailed comparison of these models with the model presented here is provided in Section S2 of the Supplementary Materials.

3.2. Computational Model

The apixaban PBPK/PD model was developed, optimised, and evaluated using the curated dataset. The model comprises a whole-body framework with dedicated submodels for the intestine, liver, and kidney (Figure 1A), which together describe apixaban absorption, distribution, metabolism, and excretion. A separate pharmacodynamic submodel represents the effect of apixaban on the coagulation pathway and the corresponding pharmacodynamic outputs, including PT, mPT, INR, aPTT, and anti-Xa activity. The model follows a hierarchical structure, in which the whole-body model connects the individual organ submodels through the systemic circulation, with the coagulation submodel embedded within this framework. Visualisations of the individual submodels are provided in Supplementary Materials Section S3 (Figures S2–S4), and the full set of model equations is reported in Section S4.

The final optimised parameter set is summarised in Supplementary Materials Table S5. Goodness-of-fit metrics and the distribution of the resulting weighted residuals are shown in Supplementary Materials Figures S5–S7. Comparisons of individual in silico simulation experiments with experimental pharmacokinetic and pharmacodynamic data not included in the main text are provided in Section S6.3 (Figures S10–S30).

The complete model, including simulation scripts, documentation, and metadata linking model variables to established ontologies, is available in SBML format under a CC-BY 4.0 license via GitHub (https://github.com/matthiaskoenig/apixaban-model) and archived on Zenodo (v0.6.0) [46].

3.3. Dose Dependency

The dose dependency of apixaban pharmacokinetic and pharmacodynamic parameters was evaluated across an oral tablet dose range of 0.5 to 100 mg. The corresponding parameter scans and time courses are shown in Figure 2, Figure 3, Figures S8, and S9A. All studies in the curated dataset that investigated apixaban dose dependency were included in the analysis (Abdollahizad2025 [70], Bashir2018 [71], Chang2016 [72], Cui2013 [73], Frost2013 [61], Frost2013a [74], Frost2014 [75], Frost2014a [76], Frost2015 [77], Frost2015a [78], Frost2015b [79], Frost2017 [80], Frost2018 [64], Frost2021 [60], Frost2021a [81], Garonzik2019 [82], Leong2024 [83], Raghavan2009 [65], Shaikh2021 [84], Song2015 [85], Song2016 [63], Tirona2018 [69], Upreti2013a [86], Vakkalagadda2016 [87], Wang2014 [88], Wang2016 [89]).

The parameter scan was consistent with the experimental data used for model evaluation, particularly at doses below 30 mg. In this range, there is a dose-dependent linear rise in drug exposure metrics such as area under the curve approximated to infinity ($AU{C}_{inf}$) and maximal concentration ($C_{\max }$), while time to reach the maximum $C_{\max }$ and half-life of apixaban remained unchanged. At higher doses, the model tended to overestimate both $AU{C}_{inf}$ and $C_{\max }$, indicating that the saturation effect was underestimated. However, the time-related parameters remained dose-independent. The corresponding in silico and clinical concentration–time profiles supported these trends. In addition, the predicted concentration–time profile of the apixaban metabolite M1 agreed with the corresponding experimental data, including the relative differences between the two investigated apixaban doses.

The maximum ratios of INR, PT, and aPTT relative to baseline showed clear non-linear, dose-dependent behaviour, which was also supported by the corresponding time-course data. Model simulations were consistent with the experimental data, particularly for INR and mPT. At higher doses, the model predicted higher mean aPTT values than those reported in clinical studies, but the predictions remained within the 2-fold range. Anti-Xa activity was represented with a linear dose dependence, which captured the available experimental data over the dose range of 2.5–5 mg.

For multiple-dose regimens, model simulations reproduced the slight accumulation of apixaban and the corresponding increase and prolongation of pharmacodynamic effects observed in clinical studies.

3.4. Renal Impairment

Renal impairment was represented by a reduction in kidney filtration capacity. The corresponding renal-function scaling parameter was varied from 0.1 (minimal filtration) to 10.0 (hyperfiltration). Simulation results and clinical data from studies investigating apixaban pharmacokinetics and pharmacodynamics in patients with different stages of renal impairment (Chang2016 [72], Metze2021 [90], Vakkalagadda2016 [87], Wang2016 [89]) are shown in Figure 4 and Figure S9B.

Renal impairment was represented by a reduction in kidney filtration capacity. The corresponding renal-function scaling parameter was varied from 0.1, representing minimal filtration, to 10.0, representing hyperfiltration. Simulation results and clinical data from studies investigating apixaban pharmacokinetics and pharmacodynamics in patients with different stages of renal impairment (Chang2016 [72], Metze2021 [90], Vakkalagadda2016 [87], Wang2016 [89]) are shown in Figure 4 and Figure S9B.

With increasing severity of renal impairment, apixaban renal clearance decreased substantially, resulting in a moderate increase in $AU{C}_{inf}$ and half-life. In contrast, $C_{\max }$ and $T_{\max }$ were less strongly affected. The model captured the observed changes in renal clearance, $AU{C}_{inf}$, $T_{\max }$, and half-life, although it tended to underestimate $C_{\max }$. Comparisons between simulated and clinical apixaban concentration–time profiles showed similar overall trends when the variability in the clinical data was considered.

The model also described the pharmacodynamic data in this population, including the mean maximum ratios of INR, PT, and aPTT relative to baseline. In addition, it reproduced the corresponding time courses, capturing both the temporal changes in these pharmacodynamic variables and the relative differences between groups with different stages of renal impairment.

3.5. Hepatic Impairment

Hepatic impairment was represented by a reduction in functional liver volume and shunting. The corresponding scaling parameter was varied from 0.0, representing a healthy liver, to 0.7, representing moderate hepatic impairment. Simulation results and clinical data from the study investigating apixaban pharmacokinetics and pharmacodynamics in subjects with different stages of hepatic impairment (Frost2021a [81]) are shown in Figure 5 and Figure S9C.

The experimental data showed no significant differences in apixaban pharmacokinetics or pharmacodynamics between healthy subjects and subjects with mild or moderate hepatic impairment, provided that no coagulation pathway dysfunction was present. Model simulations tended to overestimate apixaban concentrations at later time points in subjects with moderate hepatic impairment, resulting in a higher $AU{C}_{inf}$ and larger relative differences between the study populations than observed clinically. Simulated cumulative urinary apixaban excretion and anti-FXa activity showed the same trend. In contrast, simulated changes in $C_{\max }$, $T_{\max }$, and half-life, as well as relative changes in INR and aPTT, were consistent with the experimental data.

3.6. Food Effect

The food effect was modelled as a decrease in the fraction of apixaban absorbed. The parameter describing differences between fasted and fed conditions was optimised using apixaban pharmacokinetic data from fed subjects. Simulation results and comparisons with experimental data are shown in Figure 6 and Figure S9D.

The optimised model reproduced the experimental pharmacokinetic data from both fasted and fed subjects, with no substantial differences between the two populations. Model predictions were within the 2-fold range for pharmacokinetic parameters ($AU{C}_{inf}$, $C_{\max }$, $T_{\max }$, and half-life), apixaban concentration–time profiles, and pharmacodynamic time courses, including PT and aPTT ratios relative to baseline and anti-FXa activity.

3.7. Body Weight Dependency

The dependence of apixaban pharmacokinetic and pharmacodynamic parameters on subject body weight was evaluated across a range of 38–100 kg and is shown in Figure 7 and Figure S9E. Studies from the curated dataset that reported participant body weight were included in the analysis (Abdollahizad2025 [70], Chang2016 [72], Cui2013 [73], Frost2013 [61], Frost2014 [75], Frost2015 [77], Frost2018 [64], Frost2021a [81], Upreti2013 [91], Tirona2018 [69], Vakkalagadda2016 [87], Wang2016 [89]).

The model captured the non-linear decrease in $AU{C}_{inf}$, $C_{\max }$, and half-life with increasing body weight, whereas $T_{\max }$ remained largely unchanged. It also described the mean relative decrease in apixaban plasma concentration–time profiles across groups with low, reference, and high body weight, although it tended to underestimate the cumulative urinary amount of apixaban in subjects with low body weight.

The predicted maximum INR and aPTT ratios relative to baseline were within the 2-fold range of the available experimental data. The model predicted a slight decrease in these parameters and a more pronounced decrease in anti-factor Xa activity with increasing body weight. The model was not consistent with the anti-factor Xa activity reported by Frost2014 [75], which it underestimated. In contrast, comparison with individual data points from Upreti2013 [91] showed good agreement between predicted and observed anti-factor Xa activity across all body-weight groups (low, reference, and high) at the lower apixaban concentration, while the model still underestimated this parameter at higher concentrations.

Table 1. Summary of curated studies. Overview of studies, curated and used in the project: identifiers, PubMed IDs, PK-DB IDs, route of apixaban administrations, dosing regimen, co-administered substance, and subjects’ characteristics, including overall health status (H), renal impairment (RI), hepatic impairment (HI), and fasting status. Available data: P = apixaban plasma, MP = apixaban plasma metabolites, U = apixaban urine, MU = apixaban urine metabolites, F = apixaban feces, MF = apixaban feces metabolites, B = apixaban bile, MB = apixaban bile metabolites, INR = international normalized ratio, PT = prothrombin time, mPT = modified)prothrombin time, aPTT = activated partial thromboplastin, FXa = (anti-)FXa activity.

Study	PubMed	PK-DB	Route	Dosing	Dose [mg]	Co-adm.	H	RI	HI	Fast	Fed	P	MP	U	MU	F	MF	B	MB	INR	PT	mPT	aPTT	FXa
Abdollahizad2025 [70]	40718445	PKDB01140	oral	single	5		✓			✓		✓
Bashir2018 [71]	29972633	PKDB01165	oral	single	10	✓	✓			✓		✓
Chang2016 [72]	26358690	PKDB01148	oral	single	10		✓	✓		✓		✓		✓						✓		✓		✓
Cui2013 [73]	24353445	PKDB01142	oral	multi	10 (BID)		✓			✓		✓		✓										✓
Frost2013 [61]	22759198	PKDB01149	oral	single	sol: 0.5, 1, 2.5 tab: 5, 10, 25, 50		✓			✓	✓	✓								✓		✓	✓
Frost2013a [74]	23451769	PKDB01147	oral	multi	2.5, 5, 10, 25 (BID); 10, 25 (QD)		✓			✓		✓								✓		✓	✓
Frost2014 [75]	25419161	PKDB01110	oral	multi	2.5 (BID)		✓			✓		✓												✓
Frost2014a [76]	24697979	PKDB01167	oral	single	10	✓	✓					✓								✓				✓
Frost2015 [77]	25573421	PKDB01098	oral	single	20		✓			✓		✓		✓						✓		✓		✓
Frost2015a [78]	25377242	PKDB01166	oral	single	10	✓	✓			✓		✓
Frost2015b [79]	25501868	PKDB01150	oral	multi	10 (QD), 50 (QD)		✓					✓	✓
Frost2017 [80]	28260951	PKDB01158	oral	single	10	✓	✓			✓		✓		✓
Frost2018 [64]	30498375	PKDB01151	oral	single	2.5, 10, 25, 50		✓			✓		✓		✓						✓		✓	✓
Frost2021 [60]	34342172	PKDB01152	oral, IV	single	PO: 5; IV: 0.5, 1.25, 2.5, 3.75, 5.0		✓			✓		✓		✓						✓		✓		✓
Frost2021a [81]	34363188	PKDB01143	oral	single	5		✓		✓	✓		✓		✓
Garonzik2019 [82]	31030414	PKDB01168	oral	single	10	✓	✓			✓		✓
Jeong2019 [92]	32055579	PKDB01169	oral	single	10	✓	✓					✓
Kreutz2017 [62]	28805299	PKDB01097	oral	multi	5 (BID)		✓				✓	✓								✓	✓	✓	✓	✓
Lenard2024 [93]	36870039	PKDB01099	oral	single	0.025	✓	✓					✓								✓			✓
Lenard2025 [94]	40290437	PKDB01100	oral	single	0.025	✓	✓					✓
Leong2024 [83]	38685874	PKDB01154	oral	single	5		✓			✓		✓
Metze2021 [90]	34097808	PKDB01137	oral	single	5 (BID)			✓				✓		✓
Mikus2019 [95]	30828771	PKDB01138	oral	single	0.050	✓	✓					✓		✓
Raghavan2009 [65]	18832478	PKDB01146	oral	single	20		✓			✓		✓	✓	✓	✓	✓	✓	✓	✓
Shaikh2021 [84]	34333583	PKDB01170	oral	single	5		✓			✓		✓
Song2015 [85]	26188837	PKDB01171	oral	single	sol: 5, 10; tab: 10		✓			✓		✓
Song2016 [63]	27292282	PKDB01172	oral	single	5, 10		✓			✓	✓	✓
Rohr2024 [96]	37568371	PKDB01101	oral	single	0.025	✓	✓			✓		✓		✓
Tirona2018 [69]	29472495	PKDB01153	oral	single	2.5	✓	✓			✓		✓		✓
Upreti2013 [91]	23488672	PKDB01144	oral	single	10		✓			✓		✓		✓										✓
Upreti2013a [86]	23637566	PKDB01173	oral	single	10	✓	✓					✓
Vakkalagadda2016 [87]	26749408	PKDB01174	oral, IV	single	PO: 10, IV: 5	✓	✓			✓		✓		✓
VandenBosch2021 [97]	33351142	PKDB01145	oral	single	2.5, 5			✓
Wang2014 [88]	24277644	PKDB01175	oral	single	20	✓	✓			✓		✓
Wang2016 [89]	26331581	PKDB01141	oral	single	5		✓	✓				✓								✓	✓		✓	✓

Table 1. Summary of curated studies. Overview of studies, curated and used in the project: identifiers, PubMed IDs, PK-DB IDs, route of apixaban administrations, dosing regimen, co-administered substance, and subjects’ characteristics, including overall health status (H), renal impairment (RI), hepatic impairment (HI), and fasting status. Available data: P = apixaban plasma, MP = apixaban plasma metabolites, U = apixaban urine, MU = apixaban urine metabolites, F = apixaban feces, MF = apixaban feces metabolites, B = apixaban bile, MB = apixaban bile metabolites, INR = international normalized ratio, PT = prothrombin time, mPT = modified)prothrombin time, aPTT = activated partial thromboplastin, FXa = (anti-)FXa activity.

Study	PubMed	PK-DB	Route	Dosing	Dose [mg]	Co-adm.	H	RI	HI	Fast	Fed	P	MP	U	MU	F	MF	B	MB	INR	PT	mPT	aPTT	FXa
Abdollahizad2025 [70]	40718445	PKDB01140	oral	single	5		✓			✓		✓
Bashir2018 [71]	29972633	PKDB01165	oral	single	10	✓	✓			✓		✓
Chang2016 [72]	26358690	PKDB01148	oral	single	10		✓	✓		✓		✓		✓						✓		✓		✓
Cui2013 [73]	24353445	PKDB01142	oral	multi	10 (BID)		✓			✓		✓		✓										✓
Frost2013 [61]	22759198	PKDB01149	oral	single	sol: 0.5, 1, 2.5 tab: 5, 10, 25, 50		✓			✓	✓	✓								✓		✓	✓
Frost2013a [74]	23451769	PKDB01147	oral	multi	2.5, 5, 10, 25 (BID); 10, 25 (QD)		✓			✓		✓								✓		✓	✓
Frost2014 [75]	25419161	PKDB01110	oral	multi	2.5 (BID)		✓			✓		✓												✓
Frost2014a [76]	24697979	PKDB01167	oral	single	10	✓	✓					✓								✓				✓
Frost2015 [77]	25573421	PKDB01098	oral	single	20		✓			✓		✓		✓						✓		✓		✓
Frost2015a [78]	25377242	PKDB01166	oral	single	10	✓	✓			✓		✓
Frost2015b [79]	25501868	PKDB01150	oral	multi	10 (QD), 50 (QD)		✓					✓	✓
Frost2017 [80]	28260951	PKDB01158	oral	single	10	✓	✓			✓		✓		✓
Frost2018 [64]	30498375	PKDB01151	oral	single	2.5, 10, 25, 50		✓			✓		✓		✓						✓		✓	✓
Frost2021 [60]	34342172	PKDB01152	oral, IV	single	PO: 5; IV: 0.5, 1.25, 2.5, 3.75, 5.0		✓			✓		✓		✓						✓		✓		✓
Frost2021a [81]	34363188	PKDB01143	oral	single	5		✓		✓	✓		✓		✓
Garonzik2019 [82]	31030414	PKDB01168	oral	single	10	✓	✓			✓		✓
Jeong2019 [92]	32055579	PKDB01169	oral	single	10	✓	✓					✓
Kreutz2017 [62]	28805299	PKDB01097	oral	multi	5 (BID)		✓				✓	✓								✓	✓	✓	✓	✓
Lenard2024 [93]	36870039	PKDB01099	oral	single	0.025	✓	✓					✓								✓			✓
Lenard2025 [94]	40290437	PKDB01100	oral	single	0.025	✓	✓					✓
Leong2024 [83]	38685874	PKDB01154	oral	single	5		✓			✓		✓
Metze2021 [90]	34097808	PKDB01137	oral	single	5 (BID)			✓				✓		✓
Mikus2019 [95]	30828771	PKDB01138	oral	single	0.050	✓	✓					✓		✓
Raghavan2009 [65]	18832478	PKDB01146	oral	single	20		✓			✓		✓	✓	✓	✓	✓	✓	✓	✓
Shaikh2021 [84]	34333583	PKDB01170	oral	single	5		✓			✓		✓
Song2015 [85]	26188837	PKDB01171	oral	single	sol: 5, 10; tab: 10		✓			✓		✓
Song2016 [63]	27292282	PKDB01172	oral	single	5, 10		✓			✓	✓	✓
Rohr2024 [96]	37568371	PKDB01101	oral	single	0.025	✓	✓			✓		✓		✓
Tirona2018 [69]	29472495	PKDB01153	oral	single	2.5	✓	✓			✓		✓		✓
Upreti2013 [91]	23488672	PKDB01144	oral	single	10		✓			✓		✓		✓										✓
Upreti2013a [86]	23637566	PKDB01173	oral	single	10	✓	✓					✓
Vakkalagadda2016 [87]	26749408	PKDB01174	oral, IV	single	PO: 10, IV: 5	✓	✓			✓		✓		✓
VandenBosch2021 [97]	33351142	PKDB01145	oral	single	2.5, 5			✓
Wang2014 [88]	24277644	PKDB01175	oral	single	20	✓	✓			✓		✓
Wang2016 [89]	26331581	PKDB01141	oral	single	5		✓	✓				✓								✓	✓		✓	✓

Figure 2. Dose-dependency of apixaban and its metabolites (pharmacokinetics). A) Graphical representation of simulations examining the dose-dependent pharmacokinetics of apixaban and its metabolites. The oral dose of apixaban was increased stepwise, with distinct colours representing each dose level. B) Comparison of simulated and calculated apixaban pharmacokinetic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), and half-life on drug’s concentration. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, IV - intravenous administration, chronic - multiple administrations (more than one), x2 - two administrations in one day, and n - number of study participants. Data was used for parameter optimisation. Taken from Abdollahizad2025 [70], Bashir2018 [71], Chang2016 [72], Cui2013 [73], Frost2013 [61], Frost2013a [74], Frost2014a [76], Frost2015 [77], Frost2015a [78], Frost2015b [79], Frost2017 [80], Frost2018 [64], Frost2021 [60], Frost2021a [81], Garonzik2019 [82], Leong2024 [83], Raghavan2009 [65], Shaikh2021 [84], Song2015 [85], Song2016 [63], Tirona2018 [69], Upreti2013a [86], Vakkalagadda2016 [87], Wang2014 [88], Wang2016 [89].

Figure 2. Dose-dependency of apixaban and its metabolites (pharmacokinetics). A) Graphical representation of simulations examining the dose-dependent pharmacokinetics of apixaban and its metabolites. The oral dose of apixaban was increased stepwise, with distinct colours representing each dose level. B) Comparison of simulated and calculated apixaban pharmacokinetic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), and half-life on drug’s concentration. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, IV - intravenous administration, chronic - multiple administrations (more than one), x2 - two administrations in one day, and n - number of study participants. Data was used for parameter optimisation. Taken from Abdollahizad2025 [70], Bashir2018 [71], Chang2016 [72], Cui2013 [73], Frost2013 [61], Frost2013a [74], Frost2014a [76], Frost2015 [77], Frost2015a [78], Frost2015b [79], Frost2017 [80], Frost2018 [64], Frost2021 [60], Frost2021a [81], Garonzik2019 [82], Leong2024 [83], Raghavan2009 [65], Shaikh2021 [84], Song2015 [85], Song2016 [63], Tirona2018 [69], Upreti2013a [86], Vakkalagadda2016 [87], Wang2014 [88], Wang2016 [89].

Figure 3. Dose-dependency of apixaban (pharmacodynamics). A) Comparison of simulated and calculated apixaban pharmacodynamic parameters with experimental mean data: maximal ratio to baseline of international normalised ratio (INRmax ratio), prothrombin time (mPTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity* (anti-Xa activity max). Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. B) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, SOL- oral solution administration, IV - intravenous administration, x2 - two administrations in one day, individ - individual data points, and n - number of study participants. Timecourse data were used for parameter optimisation, and other data for evaluation. Taken from Frost2013 [61], Frost2013a [74], Frost2014 [75], Frost2018 [64], Frost2021 [60], Wang2014 [88], Wang2016 [89].

* Experimental data point (Frost2014 [75]): 2.5 mg of apixaban, 1.12 IU/ml anti-FXa activity max.

Figure 3. Dose-dependency of apixaban (pharmacodynamics). A) Comparison of simulated and calculated apixaban pharmacodynamic parameters with experimental mean data: maximal ratio to baseline of international normalised ratio (INRmax ratio), prothrombin time (mPTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity* (anti-Xa activity max). Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. B) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, SOL- oral solution administration, IV - intravenous administration, x2 - two administrations in one day, individ - individual data points, and n - number of study participants. Timecourse data were used for parameter optimisation, and other data for evaluation. Taken from Frost2013 [61], Frost2013a [74], Frost2014 [75], Frost2018 [64], Frost2021 [60], Wang2014 [88], Wang2016 [89].

* Experimental data point (Frost2014 [75]): 2.5 mg of apixaban, 1.12 IU/ml anti-FXa activity max.

Figure 4. Influence of renal impairment on apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the influence of renal impairment on the pharmacokinetics and pharmacodynamics of apixaban. Renal function changed stepwise (from minimal filtration ability to hyperfiltration), with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic and pharmacodynamic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, renal clearance, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio) on subject’s renal function. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, MiRI, MoRI, SRI, EsRI - mild, moderate, severe, and end-stage renal impairment, respectively, normal renal function if not stated, individ - individual data points, and n - number of study participants. Timecourse data from healthy subjects were used for parameter optimisation, and other data for evaluation. Taken from Chang2016 [72], Metze2021 [90], VandenBosch2021 [97], Wang2016 [89].

Figure 4. Influence of renal impairment on apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the influence of renal impairment on the pharmacokinetics and pharmacodynamics of apixaban. Renal function changed stepwise (from minimal filtration ability to hyperfiltration), with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic and pharmacodynamic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, renal clearance, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio) on subject’s renal function. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, MiRI, MoRI, SRI, EsRI - mild, moderate, severe, and end-stage renal impairment, respectively, normal renal function if not stated, individ - individual data points, and n - number of study participants. Timecourse data from healthy subjects were used for parameter optimisation, and other data for evaluation. Taken from Chang2016 [72], Metze2021 [90], VandenBosch2021 [97], Wang2016 [89].

Figure 5. Influence of hepatic impairment on apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the influence of hepatic impairment on the pharmacokinetics and pharmacodynamics of apixaban. Hepatic function changed stepwise (from normal to moderate impairment), with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity (anti-Xa activity max) on subject’s hepatic function. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, MiRI, MoRI - mild and moderate hepatic impairment, respectively, normal hepatic function if not stated, individ - individual data points, and n - number of study participants. Timecourse data from healthy subjects were used for parameter optimisation, and other data for evaluation. Taken from Frost2021a [81].

Figure 5. Influence of hepatic impairment on apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the influence of hepatic impairment on the pharmacokinetics and pharmacodynamics of apixaban. Hepatic function changed stepwise (from normal to moderate impairment), with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity (anti-Xa activity max) on subject’s hepatic function. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, MiRI, MoRI - mild and moderate hepatic impairment, respectively, normal hepatic function if not stated, individ - individual data points, and n - number of study participants. Timecourse data from healthy subjects were used for parameter optimisation, and other data for evaluation. Taken from Frost2021a [81].

Figure 6. Influence of food intake (absorption) on apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the influence of the apixaban’s absorbed fraction on its pharmacokinetics and pharmacodynamics. The fraction absorbed changed stepwise (from minimal to full absorption), with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity (anti-Xa activity max) on subject’s hepatic function. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration and n - number of study participants. Timecourse pharmacokinetic data from fasted and fed subjects were used for parameter optimisation, and other data for evaluation. Taken from Frost2021a [61], Song2016 [63], Kreutz2017 [62].

Figure 6. Influence of food intake (absorption) on apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the influence of the apixaban’s absorbed fraction on its pharmacokinetics and pharmacodynamics. The fraction absorbed changed stepwise (from minimal to full absorption), with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity (anti-Xa activity max) on subject’s hepatic function. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration and n - number of study participants. Timecourse pharmacokinetic data from fasted and fed subjects were used for parameter optimisation, and other data for evaluation. Taken from Frost2021a [61], Song2016 [63], Kreutz2017 [62].

Figure 7. Body weight-dependency of apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the body weight influence on pharmacokinetics and pharmacodynamics of apixaban. Body weight was increased stepwise, with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic and pharmacodynamic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity (anti-Xa activity max) on subject’s body weight. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, low, reference or if not stated, high - subjects with low, normal, and high body weight, respectively, individ - individual data points, and n - number of study participants. Timecourse data were used for parameter optimisation, and other data for evaluation. Taken from Abdollahizad2025 [70], Chang2016 [72], Cui2013 [73], Frost2013 [61], Frost2014 [75], Frost2015 [77], Frost2018 [64], Frost2021a [81], Upreti2013 [91], Tirona2018 [69], Vakkalagadda2016 [87], Wang2016 [89].

Figure 7. Body weight-dependency of apixaban (pharmacokinetics and pharmacodynamics). A) Graphical representation of simulations examining the body weight influence on pharmacokinetics and pharmacodynamics of apixaban. Body weight was increased stepwise, with distinct colours representing each level. B) Comparison of simulated and calculated apixaban pharmacokinetic and pharmacodynamic parameters with experimental mean data: dependency of area under the curve approximated to infinity (AUC), maximal concentration ($C_{max}$), time to reach maximal concentration (Max. time), half-life, maximal ratio to baseline of international normalized ratio (INRmax ratio), prothrombin time (PTmax ratio), activated partial thromboplastin time (aPTTmax ratio), maximal value of anti-Xa activity (anti-Xa activity max) on subject’s body weight. Standard deviations (solid) and minimal and maximal values (dashed) were plotted when available. Data was used for model evaluation. C) Comparison of model simulations (sim) with study data (exp), including standard deviations (SD) when available. PO - oral tablet administration, low, reference or if not stated, high - subjects with low, normal, and high body weight, respectively, individ - individual data points, and n - number of study participants. Timecourse data were used for parameter optimisation, and other data for evaluation. Taken from Abdollahizad2025 [70], Chang2016 [72], Cui2013 [73], Frost2013 [61], Frost2014 [75], Frost2015 [77], Frost2018 [64], Frost2021a [81], Upreti2013 [91], Tirona2018 [69], Vakkalagadda2016 [87], Wang2016 [89].

4. Discussion

In line with the aims of this study, we conducted a systematic literature review and identified 35 clinical studies reporting pharmacokinetic and pharmacodynamic data for apixaban. These studies included data from healthy adults and covered dose dependency, multiple-dose administration, food intake, body weight, renal impairment, and hepatic impairment. This work resulted in a comprehensive curated database, including data extraction, digitisation, standardisation, and independent verification of curation accuracy. Compared with previously published models, this approach enabled the use of a broader evidence base, including not only apixaban plasma concentration–time profiles, but also urinary and faecal excretion, metabolite data, and all available pharmacodynamic outputs. Bile metabolite concentrations, which were available from one study, were curated but not used for model optimisation.

Using the curated data, we developed a digital twin of apixaban pharmacokinetics and pharmacodynamics as a comprehensive PBPK/PD model with explicit metabolite representation and an expanded pharmacodynamic component. The model structure balances the available data with biological complexity, avoiding both under- and over-parametrisation. The framework describes dose-dependent apixaban absorption, tissue distribution, metabolism, and excretion. Because concentration–time profiles were available only for M1, whereas the remaining metabolite data consisted of urinary and faecal excretion of M1, M2, and M7, the metabolic pathway was simplified and limited to these metabolites. Conversion reactions were represented explicitly, but without separate enzyme-specific activities. This simplification is supported by studies reporting no strong evidence that genetic polymorphisms in genes involved in apixaban metabolism have clinically relevant effects on apixaban pharmacokinetics [67,68,69]. However, this currently limits the use of the model for mechanistic drug–drug interaction simulations. The corresponding drug–drug interaction data were nevertheless curated, published in the PK-DB database (https://pk-db.com), and visualised together with other simulations in Supplementary Materials Section S6.3. These data provide a basis for future model extensions targeting this application.

The coagulation submodel includes several pharmacodynamic outputs: international normalised ratio (INR), prothrombin time (PT), modified prothrombin time (mPT), activated partial thromboplastin time (aPTT), and anti-factor Xa activity reported in two units. The primary objective was to represent the pharmacodynamic effects of apixaban rather than to mechanistically describe the coagulation cascade itself. Therefore, these outputs were modelled as direct functions of apixaban concentration, without activated factor X as an intermediate variable, and were treated as independent of one another. This simplified the model and reduced uncertainty in parameterisation without compromising model performance.

Model performance was assessed by comparing simulations with clinical pharmacodynamic data, including both the subset used for fitting and independent evaluation data. In both cases, simulations agreed with the observed data and remained within the 2-fold range for most comparisons. In future work, INR and anti-factor Xa activity reported in mass concentration units could be modelled as functions of PT and anti-factor Xa activity reported in IU per volume, respectively, to account for assay calibration differences between laboratories and studies. In addition, no study was identified that systematically investigated the dose dependency of anti-factor Xa activity for apixaban. This limited the current model to a linear relationship between apixaban concentration and anti-factor Xa activity, although a Michaelis–Menten relationship may provide a more realistic description once suitable data become available.

All fitted parameters remained within their optimisation bounds and did not reach imposed limits. The optimisation bounds for the conversion rate constant from apixaban to M2 were set four-fold lower than those for the conversion rate constant from apixaban to M7, reflecting relative conversion rates reported in the literature [66]. Nevertheless, the optimised value for conversion to M2 was higher than that for conversion to M7. This discrepancy is likely explained by structural simplifications in the metabolic submodel and the absence of M7 concentration–time profiles. In addition, the conversion rate from apixaban to M2 was informed indirectly by the available plasma concentration–time data for M1, which is formed from M2. Expanding the database with additional metabolite concentration–time profiles would enable a more robust assessment of metabolite parameters and dynamics.

The in silico experiments reproduced the dose dependency of apixaban pharmacokinetics and pharmacodynamics, as well as the concentration–time profile of the metabolite M1. Dose-dependent absorption captured the approximately linear increase in apixaban plasma concentrations at lower doses and the modest saturation effect observed at higher doses. At high doses of 30–50 mg, simulations tended to underestimate this saturation effect and therefore overestimated apixaban plasma concentrations and cumulative urinary excretion. However, predictions remained within the 2-fold range. Moreover, these high doses are uncommon in clinical practice, supporting the prioritisation of model performance in the lower, clinically more relevant dose range [6]. For multiple-dose regimens, the model also captured the slight accumulation of apixaban and the corresponding prolongation of pharmacodynamic effects observed in clinical studies.

Because apixaban is partly cleared by the kidneys, renal impairment is expected to affect its pharmacokinetics and pharmacodynamics. The mechanistic representation of reduced renal excretion enabled the model to describe mean changes in apixaban plasma concentrations in patients with mild and moderate renal impairment. In contrast, the model underestimated apixaban plasma concentrations in patients with severe and end-stage renal impairment, likely reflecting the substantial inter-individual variability reported in the corresponding clinical studies. Nevertheless, the model reproduced pharmacodynamic changes across all stages of renal impairment, including the increase and prolongation of apixaban effects on INR, PT, and aPTT. These findings support renal function as a key factor to consider in apixaban therapy, particularly during repeated dosing, where tissue accumulation may further increase exposure and pharmacodynamic effects.

To model the effects of hepatic impairment on the pharmacokinetics and especially pharmacodynamics of apixaban is challenging to represent mechanistically. The liver not only contributes to apixaban metabolism and metabolite excretion, but also synthesises and clears coagulation factors, including factor X, whose activated form is the pharmacological target of apixaban. In cirrhosis, coagulation factors may be in an unstable balance that can shift towards either thrombosis or bleeding, substantially complicating anticoagulation therapy [98].

Only one clinical study investigating apixaban in subjects with hepatic impairment was identified, and this study excluded subjects with significant coagulopathy [81]. Therefore, severe hepatic impairment was not modelled, and model evaluation was performed under the same assumption of preserved coagulation function as in the clinical study [81]. Under these assumptions, in silico experiments tended to overestimate the relative changes in apixaban plasma concentrations and cumulative urinary excretion in subjects with moderate hepatic impairment. In contrast, mild hepatic impairment was well described, as were the relative increases in INR and aPTT and the absolute increase in anti-factor Xa activity across impairment stages. The experimental data were highly variable and included unexpected observations, such as lower cumulative urinary apixaban levels in subjects with moderate hepatic impairment than in those with mild hepatic impairment, likely reflecting the complexity of the underlying processes [81]. Overall, the simulations support the need for careful management of anticoagulation therapy in patients with hepatic impairment, while further evaluation using additional clinical studies and extension to severe hepatic impairment remain important future goals.

No clinically meaningful effect of food intake on apixaban pharmacokinetics or pharmacodynamics has been reported [61,62,63]. In the model, the parameter describing changes in apixaban absorption under fed conditions was optimised using apixaban pharmacokinetic data from fed subjects. This approach enabled reproduction of the experimental pharmacokinetic data and was further evaluated using pharmacodynamic outputs. The model reproduced the aPTT ratio and anti-factor Xa activity, while tending to overestimate PT changes. Overall, the agreement between model simulations and experimental data was good, particularly considering the high variability of PT data and the limited availability of corresponding fasted-state data used for parameter optimisation.

Although the current approach captures the mean experimental data and individual data, when data such as body weight, GFR, etc., are available, it cannot account for variability arising from inter- and intra-individual differences among subjects when this data is not reported. Future work, in addition to expanding the clinical database and other goals mentioned above, could incorporate parameter variability based on typical population distributions, moving beyond mean-based simulations and thereby substantially improving the model’s utility.

Although the current approach captures mean experimental data and can incorporate individual-level information when available, such as body weight or GFR, it cannot fully account for inter- and intra-individual variability when these data are not reported. Future work should therefore expand the clinical database and incorporate parameter variability based on representative population distributions. This would move the model beyond mean-based simulations and substantially improve its utility for virtual population studies and personalised predictions.

While several PBPK and PBPK/PD models for apixaban have been published, many are not easily reproducible, mainly because they were developed in commercial software with closed-source code and model equations, making verification challenging. As a result, it can be difficult to analyse significant differences between models that share the same objectives and framework but were created using different software [99]. Moreover, an executable model file is not shared, and calibration datasets are not reported [15,16]. Altogether, these limit the transparency, reproducibility, reusability, and cumulative development of models in the field, thereby reducing the feasibility of their translation into practice. The primary focus of this work was to address these key issues, treating them as fundamental principles throughout the project. The model is encoded in the SBMl and follows a hierarchical modular structure, allowing it to be reused not only as a whole, but also as independent submodels.

Several PBPK and PBPK/PD models of apixaban have been published, but many are difficult to reproduce because they were developed in commercial software environments with closed-source code or inaccessible model equations. This limits independent verification and makes it difficult to compare models that share similar objectives and structures but were implemented in different software platforms [99]. In addition, executable model files and calibration datasets are often not shared [15,16]. Together, these limitations restrict transparency, reproducibility, reusability, and cumulative model development, thereby reducing the feasibility of translation into practice. Addressing these issues was a central objective of this work. The model is encoded in SBML and follows a hierarchical modular structure, allowing reuse both as a complete PBPK/PD model and as individual submodels.

In summary, this digital twin of apixaban pharmacokinetics and pharmacodynamics integrates diverse clinical data into a comprehensive mechanistic PBPK/PD framework that supports mean and individual simulations across a broad range of dosing regimens and patient populations. The model provides quantitative insights into clinically relevant scenarios, including increased body weight, renal impairment, and hepatic impairment. The complete modelling framework, including curated clinical data, simulation scripts, model files, and documentation, is structured according to FAIR principles [14,17]. It is available under permissive MIT and CC-BY licenses, supporting transferability, reuse, and future model extensions with potential scientific, clinical, and industrial impact.