Exploring Tafamidis Effects Through PBPK-QSP Modelling

1. Introduction

Transthyretin amyloidosis (ATTR) is a progressive, life-threatening disorder caused by destabilization of transthyretin (TTR), a homotetrameric protein primarily synthesized in the liver [1,2]. Tetramer dissociation into monomers is the rate-limiting step of amyloid formation. Liberated monomers can misfold, aggregate, and deposit in tissues such as the myocardium and peripheral nerves, giving rise to ATTR cardiomyopathy (ATTR-CM) or polyneuropathy. The disease occurs in hereditary (ATTRv) or wild-type (ATTRwt) forms and is associated with substantial mortality. The introduction of kinetic stabilizers that bind to the thyroxine-binding sites of TTR and inhibit tetramer dissociation, such as tafamidis, has significantly improved outcomes in ATTR cardiomyopathy [3,4,5].

A consistent pharmacodynamic (PD) signature of tafamidis therapy is an increase in circulating TTR concentrations, typically around 33% above baseline in both ATTRv and ATTRwt. Although this rise is clinically reproducible, its mechanistic basis remains unresolved [6,7]. Stabilization reduces tetramer dissociation, but it is unclear whether, under physiological, non–mass-conserved conditions, this effect can fully account for the magnitude of the observed increase [8,9,10]. Therefore, mechanistic understanding of the origin of the TTR concentration rise requires models that incorporate both drug–protein binding and TTR turnover physiology.

Several TTR-focused modeling approaches have been proposed over the past decade. Subunit-exchange kinetic models accurately quantify tetramer dissociation under physiological conditions using in vitro or ex vivo plasma experiments [8]. Minimal turnover models couple synthesis, dissociation, reassociation, and degradation of TTR and show that slowed dissociation alone cannot reproduce the >30% clinical rise in circulating TTR without invoking additional mechanisms, such as altered clearance or changes in monomer disposition [9,10]. Population PK–PD models, most notably Tess et al. (2023), integrate in vitro binding, TTR concentrations, and clinical endpoints from the Transthyretin Amyloidosis Cardiomyopathy Clinical Trial (ATTR-ACT) to link tafamidis exposure with TTR stabilization and disease progression [5,11]. These analyses confirm near-complete stabilization at the approved 80 mg tafamidis meglumine dose and provide insight into the relationship between binding-site occupancy and clinical biomarkers.

Physiologically based pharmacokinetic (PBPK) modeling provides a mechanistic representation of drug disposition across organs and tissues, accounting for absorption, distribution, metabolism, and excretion based on anatomical and physiological principles [12]. Quantitative systems pharmacology (QSP) extends this to molecular and cellular processes underlying drug action (here, ligand–protein binding, tetramer stabilization, monomer generation and loss, and TTR turnover) [13]. An integrated PBPK–QSP model is therefore uniquely suited to translate tafamidis pharmacokinetics into mechanistic predictions of TTR homeostasis, offering the ability to formulate and test biological hypotheses and identify parameters that drive uncertainty.

As mentioned above, the study by Tess et al. established an important link between tafamidis binding and disease progression by modeling TTR concentrations and clinical endpoints as a function of steady-state tafamidis exposure [11]. While valuable, this approach relies on population PK rather than multiorgan PBPK, and thus does not capture tissue-level distribution, dynamic binding competition, or physiological constraints on tafamidis disposition. Additionally, the model proposed by Tess et al. explicitly assumes that monomer degradation greatly exceeds reassociation and fits several kinetic parameters simultaneously to reproduce the ~30% TTR increase observed after treatment initiation.

In contrast, the PBPK–QSP framework proposed in this manuscript incorporates experimentally supported parameters wherever possible, including human in vivo TTR turnover data and in vitro kinetic constraints from subunit exchange assays [8,14,15]. This reflects a key conceptual stance: mechanistic models must respect the origin and quality of parameter estimates, rather than expanding the parameter set to force agreement with a single clinical time point. Any sufficiently flexible model can be tuned to match sparse observations; explanatory value arises predominantly when a model reproduces data based on a physiological description, and not because its degrees of freedom are abundant. Our integration therefore makes explicit which physiological measurements, such as organ-level clearance, monomer internalization, and elimination pathways, remain critical gaps for better mechanistic understanding of tafamidis-modified TTR dynamics.

The objective of this study was to develop and validate an integrated PBPK–QSP model of tafamidis kinetics and TTR turnover. The developed model: (1) captures tafamidis pharmacokinetics across organs and tissues; (2) quantitatively describes drug binding, stabilization, and TTR turnover under physiological, non–mass-conserved conditions; (3) enables mechanistic hypothesis testing to evaluate which biological processes can plausibly account for the clinically observed increase in circulating TTR; and (4) identifies physiological and kinetic parameters that constrain our mechanistic understanding of TTR homeostasis under stabilizer therapy, thereby mapping explicit blind spots, i.e., parameters and processes we know we do not yet know well enough, to guide future experimental and clinical data collection.

2. Materials and Methods

2.1. Clinical Trials Data Used for Model Development and Validation

Clinical pharmacokinetic data from selected tafamidis studies were used to develop and evaluate model performance against observed human exposure across a range of dosing conditions. Study characteristics and data sources are summarized in Table 1.

Model accuracy was assessed using fold error, defined as the ratio of the predicted to the observed pharmacokinetic summary statistic ($Predicted/Observed$), and reported for the corresponding mean, median, or geometric mean values, as applicable. To provide a global assessment of model bias and precision, Average Fold Error (AFE) and Absolute Average Fold Error (AAFE) were additionally calculated using fold errors derived from exposure metrics (AUC and $C_{max}$) only, consistent with standard PBPK model evaluation practice. AFE was defined as the geometric mean of the Predicted/Observed ratios, calculated as ${10}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\left(\mathrm{l}\mathrm{o}\mathrm{g}\right(\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}/\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d}\left)\right)}$, while AAFE was defined as the geometric mean of the absolute fold deviations, calculated as ${10}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}(\mid \mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}/\mathrm{O}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d})\mid )}$. Model performance was interpreted against commonly applied PBPK acceptance criteria, with fold error values within 0.5–2.0 considered acceptable and values within 0.8–1.25 considered indicative of strong predictive performance.

2.2. PBPK Model Development

The general structure of the PBPK model was previously proposed and implemented in Simcyp Simulator (version 23) [20]. The model follows a whole-body, perfusion-limited PBPK framework in which all anatomical compartments are represented as well-stirred, kinetically homogeneous units, assuming instantaneous distribution within each tissue relative to inter-compartmental blood flow. Drug transport between compartments is governed by organ-specific blood flows and tissue-to-plasma partition coefficients.

The PBPK model-simulated plasma concentration of tafamidis is then used as an external, dynamic effect driving variable for the downstream QSP model. The simulated plasma concentration–time profile of tafamidis serves as a dynamic forcing function for the downstream QSP model, with no feedback from the QSP system to the PBPK component.

2.1.1. Drug Physicochemical Properties and Blood Proteins Binding

Tafamidis is a monoprotic acid with a pKa of 3.73, molecular weight of 308.12 g/mol, with the molecular weight of tafamidis meglumine salt being 503.33 g/mol [21]. A logP value (3.91), calculated using the ALOGPS model, was used to parametrize the PBPK model [22]. The plasma unbound fraction of the drug has a value of 0.008 (0.8%). Apart from its affinity for TTR, which is a key element of its mechanism of action, tafamidis also exhibits high affinity for plasma albumin, consistent with interactions between a weakly acidic small molecule and proteins with basic binding domains [23]. Binding to plasma proteins is part of the developed QSP model, as described further in the QSP model development and Results sections.

2.1.2. Dissolution and Absorption

Tafamidis has been reported as having low solubility and variable permeability; however, in vitro-in vivo extrapolation of Caco-2 permeability results suggests that it is actually a well-permeable drug, since its apparent permeability ($P_{app}$) equals $30.90\pm 3.28\times {10}^{-6} {\displaystyle \frac{m}{s}}$ from the apical to basolateral side and $31.90\pm 2.82\times {10}^{-6} {\displaystyle \frac{m}{s}}$from basolateral to apical side [24]. The above-mentioned $P_{app}$ value was used to estimate effective permeability ($P_{eff}$) using linear scaling, as described by Sun et al. [25].

According to regulatory data, tafamidis meglumine solubility at pH 6.8 buffer is 3.121 ${\displaystyle \frac{mg}{mL}}$. It is also known that in an alkaline environment (0.1 M NaOH), solubility is higher than 4.19 ${\displaystyle \frac{mg}{mL}}$, while solubility in water is higher than 4.63 ${\displaystyle \frac{mg}{mL}}$ [26]. The presence of meglumine, an alkalising agent, explains the higher solubility in pure water compared to the pH 6.8 buffer. In the developed PBPK model, tafamidis meglumine solubility was described as a pH-dependent parameter, incorporating the adjusted Henderson-Hasselbalch equation for prediction, based on the measured value at pH 6.8 [27]. Predicted intrinsic solubility, S_o = 0.0025 mg/mL, seems reasonable considering low solubility in low pH. The predicted pH dependent solubility is presented in Figure S1. The commercially available tafamidis meglumine is manufactured as gelatine-encapsulated powder; thus, it is assumed that it is released rapidly once it enters the duodenum, allowing for the immediate start of the dissolution process and absorption. For description of the absorption process, the ADAM (Advanced Dissolution, Absorption, and Metabolism) model was used, with the stomach serving as the dose-accepting compartment and upstream driver of intestinal absorption. Fed and fasted conditions were handled implicitly within the Simcyp ADAM framework by applying condition-specific gastrointestinal physiology (including gastric emptying time, luminal fluid volumes, bile salt concentrations, and intestinal transit), as implemented in the Simcyp Healthy Volunteer population. The above-reported Caco-2 $P_{app}$ data were used to estimate $P_{eff}$ as $6.97\times {10}^{-4}{\displaystyle \frac{cm}{s}}$[28]. It was assumed that permeability is the same for each of the seven intestinal subcompartments, whereas local pH, fluid dynamics, GIT transit time follow inter-individual variability as implemented in the Simcyp Healthy Volunteer population [20]. Additionally, the gastric emptying time varies between virtual individuals [29].

2.1.3. Volume of Distribution and Tissue-to-Plasma Partitioning

The apparent volume of distribution after oral administration of tafamidis meglumine to healthy volunteers was reported as 16 L (0.2 L/kg for an 80 kg patient) [30]. For a full PBPK distribution model, the volume of distribution at steady state ($V_{ss}$), accounting for distribution in the following tissues: adipose, bone, brain, gut, pancreas, heart, kidney, liver, lung, muscle, skin, spleen, was calculated following the Sawada equation [31].

Partition coefficients between each tissue and plasma ($K_p$) – except for bone, pancreas, muscle, and skin – were obtained from rat PK study on tafamidis distribution, following the assumption that drugs are distributed in rats similarly to humans [24,32]. For prediction of $K_p$ for bone, pancreas, muscle, and skin, the Rodgers and Rowland model was incorporated as described previously [33]. The model estimated $V_{ss}$as 13.53 L for an 80 kg patient (0.17 L/kg).

2.1.4. Metabolism and Excretion

Tafamidis is metabolized primarily by UGT-mediated glucuronidation, predominantly by UGT1A1, UGT1A3, and UGT1A9 [34]. Unfortunately, no intrinsic clearance values or contributions of each of the enzymes can be found in the publicly available literature; thus, metabolism was described with a holistic parameter of oral clearance – 0.263 L/h – as stated in drug prescribing information [35]. Therefore, the model is not suitable for populations with altered UGT abundance (i.e. certain diseases including liver cirrhosis) or for assessing drug–drug interactions, as it lacks the mechanistic scalability offered by intrinsic clearance data. Based on documentation submitted to the FDA, a large fraction of pharmacologically inactive tafamidis metabolite is expected to be excreted with bile; the proposed model does not account for this route of elimination as there is not enough data to mechanistically describe such a process.

2.1.5. Enzyme- and Transporter-Mediated Drug–Drug Interactions

In vitro data presented in the regulatory documents do not exclude inhibition of intestinal UGT1A1 and induction of CYP2B6 and CYP3A4. It has been suggested that tafamidis might also inhibit OAT1, OAT3, and BCRP transporters [34]. Nevertheless, interaction with enzymes and transporters is not relevant for this model as drug-drug interactions (DDI) are not in scope of this study.

2.2. QSP Model Development

Tafamidis drug action was described using a basic ligand–protein binding model assuming two binding sites on the TTR molecule. The model assumes that after tafamidis binding, TTR undergoes conformational changes as previously observed [36,37]. Since structural modifications might alter molecular properties, the bound protein is hypothesized to exhibit distinct attributes compared to its unbound form, making it less susceptible to clearance or more susceptible to undergo tissue internalization [6,37,38]. In our analysis, we assumed that a kinetically stabilised TTR molecule does not dissociate to monomers at all. It is known that the dissociation rate of total TTR in plasma decreases as tafamidis concentration increases; however, it is not possible to specify quantitatively how various bound states of TTR contribute to this change [8]. The total TTR concentration remains at a dynamic steady-state level governed by its synthesis and degradation processes. These processes are quantified by rates and can be described using a zero-order synthesis and first-order degradation turnover framework as previously reported [9]. Under this formulation, total TTR dynamics can be expressed as:

$$\begin{array}{c}{\displaystyle \frac{dTTR}{dt}}\left[\mu M\right]=k_{synthesis}\left[{\displaystyle \frac{\mu M}{day}}\right]-k_{degradation }\left[{\displaystyle \frac{1}{day}}\right]\times TTR \left[\mu M\right],\end{array}$$

(1)

yielding a steady-state concentration of:

$$\begin{array}{c}TT{R}_{baseline}\left[\mu M\right]={\displaystyle \frac{k_{synthesis}\left[\frac{\mu M}{day}\right]}{k_{degradation}\left[\frac{1}{day}\right]}} .\end{array}$$

(2)

Consequently, any change in the steady-state concentration of circulating TTR must arise from a change in the ratio between its synthesis and elimination rates. Let ${\mathrm{T}\mathrm{T}\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}}$ denote the baseline steady-state TTR concentration prior to tafamidis administration, and ${\mathrm{T}\mathrm{T}\mathrm{R}}_{\mathrm{s}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{w}}$ the steady-state concentration after drug saturation. An observed 33% increase corresponds to:

$$\begin{array}{c}{\mathrm{T}\mathrm{T}\mathrm{R}}_{\mathrm{s}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{w}}\left[\mu M\right]=1.33\times {\mathrm{T}\mathrm{T}\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}}\left[\mu M\right].\end{array}$$

(3)

When monomer elimination is negligible, the steady-state tetramer concentration is governed by the balance between synthesis and degradation. In this framework, proportional scaling of either the synthesis rate or the degradation rate leads to mathematically identical steady-state solutions. Consequently, steady-state data alone cannot distinguish whether the observed increase in circulating TTR arises from enhanced synthesis or reduced elimination; these mechanisms bring structural non-identifiability of the model. This structural ambiguity is not specific to the present PBPK-QSP framework but reflects a general limitation of turnover models constrained solely by steady-state observations, as discussed previously in a phenomenological context [10].

To reproduce the observed population-average increase while keeping synthesis unchanged, we introduced a scaling factor $s$ on the effective degradation rate such that:

$$\begin{array}{c}{\displaystyle \frac{k_{\mathrm{s}\mathrm{y}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}} \left[\frac{\mu M}{day}\right]}{s\times k_{\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}} \left[\frac{1}{day}\right]}}=1.33\times {\displaystyle \frac{k_{\mathrm{s}\mathrm{y}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}} \left[\frac{\mu M}{day}\right]}{k_{\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}} \left[\frac{1}{day}\right]}},\end{array}$$

(4)

yielding the solution where $s=1/1.33\approx 0.75$.

Tafamidis binds directly to circulating TTR tetramers and induces stabilizing conformational changes, which in the model are assumed to alter internalization or clearance pathways of the protein; this assumption is supported by the high fraction of circulating TTR present in the drug-bound state.. In contrast, increased hepatic TTR synthesis is unlikely to underlie the observed rise in circulating TTR, as in vitro evidence supports post-translational modulation of secretion at the site of synthesis, requiring intracellular (endoplasmic reticulum) access; however, these findings do not necessarily translate quantitatively to the in vivo setting [39]. Accordingly, the QSP model is calibrated to a known population-level steady-state effect and is not intended to provide predictive discrimination between synthesis- and clearance-based mechanisms.

The main differences between the current and previously presented minimalistic QSP model include the application of two stable bound states of transthyretin and the introduction of a drug concentration variable as an input from the above-described PBPK model [9]. The previously presented minimalistic model was developed to provide a conceptual and mathematical framework for understanding TTR dynamics in vivo, with particular emphasis on disentangling the roles of tetramer stabilization, synthesis, and elimination. The model describes TTR as an open system with production, reversible tetramer–monomer interconversion, and irreversible elimination from the circulation. Under its core assumptions—rapid monomer–tetramer reassociation relative to monomer degradation—the model demonstrates analytically that changes in tetramer dissociation rate alone do not alter steady-state TTR concentration. Instead, clinically observed increases in circulating TTR can only be reproduced by concurrent changes in TTR production and/or elimination rates. Thus, the purpose of the minimal model development was mechanistic understanding to rationalize clinical observations and identify which physiological processes must be affected by tetramer stabilizers to explain increased serum TTR levels. The current work builds upon this framework by extending it to a mechanistic, state-resolved QSP model that explicitly represents multiple drug-bound TTR states, albumin binding, and dynamic tafamidis exposure provided by a PBPK model, thereby enabling quantitative simulation and prediction of clinically observed TTR responses. The general structure of the proposed model can be seen in Figure 1.

The analysis of drug effect relies on precise description of total TTR concentration (${TTR}_{total}$) in plasma, as it is the only surrogate of drug action that can be quantitatively tracked. To characterize the sensitivity of model-predicted TTR responses to tafamidis dose, a dose-sensitivity analysis was performed under steady-state conditions by systematically varying the administered tafamidis meglumine dose around the recommended clinical regimen (4 × 20 mg, corresponding to 48.8 mg tafamidis). Simulations were conducted for fractional and multiple doses equal to twofold, three-quarters, one-half, and one-quarter of the standard dose, selected to span a clinically relevant and practically divisible range. For each dose level, steady-state values of total, bound, free, and monomeric TTR species were extracted, with particular emphasis on the unbound-to-total TTR ratio as an integrative, model-derived quantity summarizing dose-dependent changes in TTR binding at the population level.

We suggested a mechanism combining TTR degradation rate modification with molecular-level TTR binding. However, it should be noted that in some cases, steady-state TTR levels may increase through stabilization alone, without affecting degradation, although such increases are typically modest, not exceeding, on average, 15% [10].

The full model code, in the form of a notebook, together with parameter descriptions, can be found in the Supplementary Materials. We provide results for models with two different association and dissociation constants, coming from separate experiments, described in two articles; the code can be easily adjusted to switch from one set of parameters to another [37,40]. Results of the simulations are available as Excel files attached as the Supplementary Materials.

2.2.1. QSP Model Parametrisation

While a static model would be sufficient for monitoring steady-state concentrations, incorporation of dynamics enhances the model. It allows tracking delays in effect and accounting for time-dependent variations in individual drug concentrations. A total of six ordinary differential equations (ODEs) were used in the model (see equations below). The QSP model parameters are categorised into two groups: drug-specific parameters and system (patient)-specific parameters. The latter, namely system parameters, were provided as individual values for virtual individuals as detailed in a subsequent section. After analysis of the binding affinities, two sets of association and dissociation rate constants were applied: one derived from the equilibrium binding constant reported in Nelson et al. and the other one directly adopted from Corazza et al. [37,40]. To account for the observed TTR dynamics, a fixed scaling factor was introduced to describe changes in the clearance of the drug–protein complex. Specifically, the model was calibrated using scaling factors of 0.75 for TTR1 and TTR2 clearances, which allows the model to predict an average increase of 33% in TTR plasma concentration in healthy male volunteers, as previously reported in the literature [6].

Identifiability of the synthesis and clearance scaling factors was evaluated using a local Fisher Information Matrix (FIM)–based analysis with a steady-state pseudo-observation reflecting the reported 33% increase in total TTR concentration. Results of the analysis demonstrated that synthesis and clearance scaling factors cannot be estimated independently using available data, exhibiting strong correlation and a rank-deficient FIM. Consequently, the attribution of the observed steady-state increase to reduced effective clearance represents a mechanistically motivated modeling assumption rather than a statistically identifiable fact. Full details of the identifiability analysis are provided in the Appendix (R scripts).

Generative AI support (GPT-5.2 large language model, OpenAI; accessed January 2026) was used in a limited technical-assistance role during development of the Fisher Information Matrix identifiability workflow. Specifically, the model was used to help for- malize the computational procedure for local FIM construction, including symbolic-to- numeric mapping of model parameters, structuring of finite-difference sensitivity calcu- lations, and implementation logic for the steady-state pseudo-observation approach in R. The AI tool did not generate primary scientific results, select model assumptions, or inter- pret identifiability outcomes; all equations, parameterizations, scripts, and conclusions were defined, verified, and validated by the authors.

All parameters used for building the QSP model are presented in the Supplementary Material (see Table A2).

2.2.2. Equations Describing Tafamidis-TTR Kinetics

$${\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}1}}{\mathrm{d}\mathrm{t}}}=-{\mathrm{T}}_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d},\mathrm{d}\mathrm{e}\mathrm{g}}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}1}+{\mathrm{k}}_{\mathrm{o}\mathrm{n}1}\times {\mathrm{C}}_{\mathrm{T}\mathrm{a}\mathrm{f}}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}-{\mathrm{k}}_{\mathrm{o}\mathrm{f}\mathrm{f}1}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}1}+{\mathrm{k}}_{\mathrm{o}\mathrm{f}\mathrm{f}2}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}2}-{\mathrm{k}}_{\mathrm{o}\mathrm{n}2}\times {\mathrm{C}}_{\mathrm{T}\mathrm{a}\mathrm{f}}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}1}$$

(5)

The concentration of single tafamidis-TTR complex increases when free transthyretin (${TTR}_{free}$) binds to a tafamidis molecule, initiating stabilization of TTR. The single tafamidis-TTR complex also rises when the double tafamidis-TTR complex dissociates, releasing one tafamidis molecule and reverting to the single-bound form. Conversely, the single tafamidis-TTR complex concentration decreases through three mechanisms: irreversible degradation of the complex, dissociation of tafamidis from TTR into their free forms, and the binding of an additional tafamidis molecule to the single complex, forming the double tafamidis-TTR complex.

In Equation (5), k_on1 represents the primary association rate of tafamidis with TTR ($\mu M^{-1}{h}^{-1}$), while $k_{off1}$ denotes the primary dissociation rate of tafamidis from TTR ($1/h$). The formation of the double complex is governed by $k_{on2}$, the secondary association rate of tafamidis with the single-bound complex ($\mu M^{-1}{h}^{-1}$), and its dissociation is described by $k_{off2}$ ($1/h$). The irreversible elimination is the same for single- and double-bound complexes, assigned as $T_{bound,deg}$ ($1/h$), which is scalar adjusted via $T_{deg}$.

$$\begin{array}{c}{\displaystyle \frac{d{C}_{TTR2}}{dt}}=-T_{bound,deg}\times C_{TTR2}+k_{on2}\times C_{Taf}\times C_{TTR1}-k_{off2}\times C_{TTR2}\end{array}$$

(6)

The concentration of the double tafamidis-TTR complex increases when a second tafamidis molecule binds to the single tafamidis-TTR complex, however without accounting for enhanced resistance to dissociation. Conversely, this concentration decreases through two primary mechanisms: irreversible degradation of the double complex and dissociation of one tafamidis molecule, which converts the double complex back into the single-bound form.

In Equation (6), $k_{on2}$ represents the secondary association rate of tafamidis with the single-bound TTR complex ($\mu M^{-1}{h}^{-1}$), while $k_{off2}$ denotes the secondary dissociation rate of Tafamidis from the double complex ($1/h$).

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}{\mathrm{d}\mathrm{t}}}=\mathrm{r}+{\mathrm{k}}_{\mathrm{a}}\times {{\mathrm{C}}_{\mathrm{M}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{r}}}^4-\left({\mathrm{k}}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}+{\mathrm{T}}_{\mathrm{d}\mathrm{e}\mathrm{g}}\right)\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}-{\mathrm{k}}_{\mathrm{o}\mathrm{n}1}\times {\mathrm{C}}_{\mathrm{T}\mathrm{a}\mathrm{f}}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}+{\mathrm{k}}_{\mathrm{o}\mathrm{f}\mathrm{f}1}\times {\mathrm{C}}_{\mathrm{T}\mathrm{T}\mathrm{R}1}\end{array}$$

(7)

The concentration of free TTR increases through its natural synthesis in the body and the dissociation of tafamidis from the single tafamidis-TTR complex, releasing TTR back into its free form, as well as due to reassociation of monomers back to tetramers with rate k_a ($\mu M^{-3}{h}^{-1}$). Conversely, free TTR is reduced by three key processes: natural degradation of TTR, dissociation into monomers, which reflects destabilization of its native tetrameric structure and increases the risk of aggregation, and binding with tafamidis to form the single tafamidis-TTR complex.

In Equation (7), $\mathit{r}$ represents the rate of TTR synthesis ($\mu M/h$), while ${\mathrm{T}}_{\mathrm{d}\mathrm{e}\mathrm{g}}$ denotes the natural degradation rate of TTR ($1/h$). The dissociation of TTR into monomers is described by ${\mathrm{k}}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}$ ($1/h$), and the binding of free TTR with tafamidis to form the single complex is governed by ${\mathrm{k}}_{\mathrm{o}\mathrm{n}1}$ ($\mu M^{-1}{h}^{-1}$). The reverse process, where tafamidis dissociates from the single complex, is described by ${\mathrm{k}}_{\mathrm{o}\mathrm{f}\mathrm{f}1}$ ($1/h$).

$$\begin{array}{c}{\displaystyle \frac{d{C}_{Monomer}}{dt}}=4k_{diss}\times C_{TTR}-4\times k_a\times {C_{Monomer}}^4-M_{deg}\times C_{Monomer}\end{array}$$

(8)

The concentration of monomers increases when free TTR dissociates from its tetrameric form, reflecting destabilization of the native structure. This dissociation is a critical step in the pathogenesis of amyloid diseases, as monomers are prone to misfolding and aggregation into amyloid fibrils. Conversely, the monomer concentration decreases when monomers re-associate to form stable tetramers, restoring TTR’s functional structure, undergo irreversible degradation, or clump to aggregates. For the proposed model the last two are indistinguishable, thus they are merged within single parameter.

In Equation (8), $k_{diss}$ represents the dissociation rate of tetrameric TTR into monomers ($1/h$), while $k_a$ denotes the association rate of unbound TTR monomers to reform tetramers ($\mu M^{-3}{h}^{-1}$). The irreversible degradation of monomeric TTR is defined by $M_{deg}$ ($1/h$). Since no data exist on the concentration of dissociated TTR monomers in plasma, we assumed it results from the steady-state balance of TTR and its dissociation and association constants as presented in the Appendix A (Table A2, Equation (A1)).

$$\begin{array}{c}{\displaystyle \frac{d{C}_{Alb1}}{dt}}=k_{on,albumin}\times C_{Taf}\times C_{Alb,free}-k_{off,albumin}\times C_{Alb1}\end{array}$$

(9)

Tafamidis can bind to albumin in the plasma, forming a tafamidis-albumin complex, thus reducing the availability of free tafamidis for transthyretin stabilisation. This interaction serves as a reservoir, modulating the pharmacokinetics of tafamidis, mainly by keeping the fraction unbound at a stable level. Conversely, the number of tafamidis-albumin complexes decreases when tafamidis dissociates from albumin, releasing free tafamidis back into the blood.

In Equation (9), $k_{on,albumin}$ represents the association rate of tafamidis with albumin ($\mu M^{-1}{h}^{-1}$), while $k_{off,albumin}$ denotes the dissociation rate of tafamidis from the tafamidis-albumin complex ($1/h$).

$$\begin{array}{c}{\displaystyle \frac{d{C}_{A\mathrm{l}\mathrm{b},\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}{dt}}=-k_{on,albumin}\times C_{Taf}\times C_{Albumin}+k_{off,albumin}\times C_{Alb1}\end{array}$$

(10)

Albumin is present in abundance; however, it is necessary to track its concentration to capture all drug interactions. We assume that there is only one binding site of tafamidis to albumin and that only unbound albumin described with Equation 10 is available for binding [11].

$$C_{Taf}=C_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-\left(C_{TTR1}+2\times C_{TTR2}+C_{Alb1}\right)$$

(11)

The concentration of free tafamidis in plasma is determined by the total administered dose, excluding the amount bound to protein complexes. We assume that only the unbound drug is available to interact with TTR. Total concentration ($C_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$) mentioned in Equation 11 represents the total tafamidis concentration treated as the merging variable for the PBPK and QSP model. All concentrations are expressed in $\mathrm{\mu }\mathrm{M}$.

2.3. Population-Level Sampling of TTR Concentrations

To conduct simulations in a virtual population, baseline TTR concentrations were generated for thirty virtual healthy male subjects aged 30–40 years, using summary statistics reported by Ingenbleek et al. for age- and sex-stratified groups [41]. This age range was selected deliberately for methodological consistency. First, the only publicly available multiple-dose tafamidis pharmacokinetic data approaching steady state were generated in healthy male volunteers within this age range and were used to evaluate the PBPK component of the model [19]. Second, the most frequently cited human estimates of transthyretin elimination kinetics were derived from in vivo radiolabeled protein kinetic studies conducted in healthy males aged 30–40 years, as reported by Oppenheimer et al. and Socolow et al. [14,15]. These studies reported first-order degradation rate constants of approximately 0.363 day⁻¹ (Oppenheimer et al.) and 0.268 day⁻¹ (Socolow et al.). Ingenbleek et al. reported that baseline TTR concentrations within each age–sex stratum are approximately normally distributed [41]. Individual baseline TTR values were sampled randomly from a normal$(\mu ,\sigma )$ distribution parameterized using the reported mean and standard deviation for this group. The number of virtual subjects (n = 30) was selected to match the size of the multiple-dose clinical pharmacokinetic study used to evaluate the PBPK component, thereby aligning simulated between-subject variability with the available observed study cohort [19]. One baseline value was assigned per virtual subject and held constant across simulations, thereby preserving the reported between-subject variability in baseline TTR concentrations.

Female subjects were not simulated in the present work. Although baseline TTR distributions by sex are available, the data required to parameterize and validate sex-specific tafamidis exposure under a multiple dosing regime and to quantitatively support transferability of the observed TTR increase across sexes are limited [41]. To avoid introducing additional unverifiable assumptions, the present analysis was restricted to males, and generalization of the results to female populations is acknowledged as a limitation.

Initial values for optimisation were derived from Oppenheimer et al., as this study reports turnover estimates in a demographic group matching the selected virtual population [14]. Parameter fitting and population sampling were implemented using the R programming language. Full optimisation results, sampling scripts, and associated code are provided in the Supplementary Material.

3. Results

The developed PBPK-QSP model effectively describes the overall drug effect on TTR and successfully replicates tafamidis pharmacokinetic outcomes observed in multiple-dose studies conducted in a healthy population. The comparison of observed and predicted PK parameters (AUC, $C_{max}$, $T_{max}$) is shown in Table 2.

The pharmacokinetic model’s performance in predicting tafamidis exposures after multiple-dose therapy in a 30–40-year-old male population is presented in Figure 2D. Moreover, simulated tafamidis plasma concentration–time profiles following a single supratherapeutic 400 mg dose were consistent with observed clinical data in a cohort of 42 healthy individuals (Figure 2A). The model also successfully captured the pharmacokinetics of single-dose administration in 16 healthy volunteers enrolled in trial B3461054, demonstrating robust predictive performance under both fasted and fed conditions (Figure 2B and 2C) [17]. Across all evaluated dosing regimens and nutritional conditions, predicted exposure metrics (AUC and $C_{max}$) showed good agreement with observations, with individual fold errors ranging from 0.74 to 1.12. Based on geometric mean AUC and Cmax values, the overall AFE was 0.93, indicating a slight underprediction, while the AAFE was 1.13, corresponding to a typical prediction error of approximately 13%.

Figure 3 illustrates the distribution of bound tafamidis across the single- and double-bound states of TTR, as well as its interaction with albumin, which acts as a buffering pool for circulating tafamidis and modulates its availability for TTR binding. Comparison of results obtained using the Nelson-based and Corazza-based binding parameter sets shows that the former yields a higher proportion of single-bound TTR–tafamidis complexes, whereas the latter shifts the distribution toward a higher abundance of double-bound complexes (Figure 3 and Supplementary Figure A2).

The steady-state unbound-to-total TTR ratio exhibited a clear dependence on tafamidis dose across the evaluated range. Relative to the recommended tafamidis meglumine regimen (4 × 20 mg; 48.8 mg tafamidis), stepwise dose reductions to three-quarters, one-half, and one-quarter of the standard dose were associated with progressively higher unbound-to-total TTR ratios, whereas a twofold dose increase resulted in a further reduction of the unbound fraction. The predicted dose–response relationship for the unbound-to-total TTR ratio is shown in Figure 4. Dose-dependent changes in total, bound, free, and monomeric TTR species are presented in the Supplementary Materials (Figure A3, Figure A4, Figure A5 and Figure A6).

4. Discussion

Tafamidis has been approved for the treatment of transthyretin amyloid cardiomyopathy based on demonstrated clinical benefit, including improved survival and reduced cardiovascular-related hospitalizations, together with strong experimental evidence showing its ability to kinetically stabilize transthyretin tetramers in vitro and ex vivo. The present modeling work does not aim to reproduce or predict clinical outcomes. Instead, it builds on the established mechanism of action of tafamidis—namely, stabilization of the transthyretin structure.

The presented PBPK-QSP model successfully replicates tafamidis pharmacokinetics in a healthy population, supporting its applicability for dose selection and clinical trial optimization. Tafamidis is a long-term therapy and accordingly the integrated PBPK–QSP model simulation was run over a sufficiently long time horizon to reach pharmacodynamic steady state. It should be noted that pivotal clinical evidence for tafamidis efficacy derives predominantly from elderly ATTR patients, often with multiple comorbidities. Age- and disease-related physiological changes may affect both tafamidis pharmacokinetics and TTR turnover; therefore, extrapolation of simulation results from young healthy males to the clinical ATTR population should be done cautiously. In the absence of age-matched longitudinal data suitable for parameterization and validation within this framework, the reported average ~33% increase in total TTR concentration was used as a calibration target, and this assumption is explicitly treated as a limitation of the model. Moreover, the model accurately predicts tafamidis plasma concentrations following a single supratherapeutic 400 mg dose (Figure 2A), supporting its applicability for dose selection.

On the PD side, the model captures trends in average TTR concentration changes; however, its predictive performance cannot currently be evaluated at the individual level. This is because the model was developed based on the data derived from population-level analysis. If the individual TTR trajectories were available, the model’s predictive accuracy for patient-specific responses could be assessed and further refined. The variability in TTR concentration appears to depend on the baseline protein levels and overall patient health, but the precise contributing factors remain unidentified. Addressing this limitation could potentially enhance the model’s predictive accuracy and increase its usefulness towards personalized medicine. Our analysis indicates that patients with low baseline TTR levels may achieve adequate stabilization at lower doses. For instance, a low-TTR patient (3.85 $\mu M$) dosed with 60 mg tafamidis meglumine reached steady-state fractional stabilization of 2.55%, comparable to the population average of 2.19%. In contrast, a high-TTR patient (7.51 $\mu M$) required 80 mg daily to reach 3.86% stabilization, a value close to that achieved by the average patient on 60 mg (4.59%). The dose levels in this analysis were chosen deliberately: while the standard regimen involves four 20 mg tablets, reducing the dose to three tablets in patients with low TTR could lower treatment costs without compromising therapeutic effect. Conversely, patients with high baseline TTR may require higher doses to achieve optimal stabilization. These simulations suggest that baseline TTR concentration may influence the exposure–stabilization relationship; however, any dose individualization based on baseline TTR would require prospective clinical validation and should be regarded as hypothesis-generating.

Dose sensitivity analysis has shown that there is no clear difference in the binding effect across the tested dose range (Figures S3–S6). However, the simulations suggest that baseline TTR may influence the exposure–stabilization relationship in extreme cases; this observation is hypothesis-generating and would require prospective clinical validation before informing dosing decisions. Patients with lower TTR levels could potentially decrease the dose or extend dosing gaps. A critical aspect of the model’s implementation is the necessity of applying a fixed scaling factor to describe changes in drug-protein complex clearance. Due to the lack of sufficient experimental data, this heuristic adjustment was required to align predictions with observed clinical data. Similarly, the rate of monomer degradation was set to a constant value in the absence of experimental data on its actual magnitude. This choice is justified by the fact that, under physiological conditions, the product of $k\_association \times M^3$ is expected to be substantially lower than $k\_deg,monomer \times M$, given that monomer concentrations are below 1% of tetramer concentrations (as reported by Sekijima et al. 2001) [42]. Consequently, variations in this parameter are unlikely to dominate system dynamics; however, this assumption should be revisited if future data indicate otherwise.

A key methodological point is that, in an open turnover system such as circulating transthyretin, a clinically observed steady-state increase in total TTR can generally be reproduced by multiple mechanistic hypotheses that act on different parts of the turnover balance. In particular, reduced tetramer dissociation, altered effective clearance/internalisation of stabilized species, and modulation of net production can lead to similar steady-state outcomes, rendering the underlying mechanism structurally ambiguous when only steady-state endpoints are available [10]. In this light, the present PBPK–QSP framework is intended not as a definitive mechanistic proof, but as a quantitative scaffold that preserves mechanistic structure where parameters are experimentally constrained (tafamidis exposure, binding stoichiometry and kinetics), while representing unresolved biology through a minimal number of calibrated terms. This modeling stance allows the integration of pharmacokinetic information with TTR turnover physiology and provides a principled way to constrain the space of plausible mechanisms consistent with observed PK and population-average TTR responses, thereby guiding which future measurements would be most informative for discriminating between competing hypotheses.

The proposed model does not uniquely identify the biological mechanism underlying the observed increase in TTR concentrations, but it constrains which turnover and clearance hypotheses are consistent with the available steady-state and pharmacokinetic data. Future research should focus on elucidating whether this phenomenon is driven by modified clearance, internalization of bound complexes, or potential effects of tafamidis on TTR synthesis rates. Additional studies are also needed to clarify albumin’s competitive role in tafamidis binding and its impact on pharmacokinetics. While the PBPK-QSP model demonstrates strong predictive capabilities, further refinements and validation studies are necessary to fully capture the complexities of tafamidis-TTR interplay. We acknowledge that the model relies on several assumptions without direct experimental validation, and further research is required to confirm its mechanistic foundations.