1 Introduction

Thyroxine (T4) and triiodothyronine (T3) are thyroid hormones (THs) that regulate several vital physiological processes: growth, differentiation, development of the nervous system, cardiovascular function, reproduction, and energy metabolism (1–5). The circulating TH levels are under robust homeostatic control to keep within narrow physiological ranges. Structured as a negative feedback loop, the hypothalamic-pituitary-thyroid (HPT) axis is the primary homeostatic mechanism. The thyrotropin-releasing hormone (TRH) produced by neurons in the hypothalamic paraventricular nucleus (PVN) stimulates the synthesis and release of the thyroid-stimulating hormone (TSH) in the anterior pituitary. TSH then stimulates the production and release of T4 and T3 by the thyroid follicles (6). Negative feedback from T4 and T3 occurs through the inhibition of both TRH and TSH, which keeps the circulating free TH levels close to the setpoints in the long term. In contrast, a separate mechanism is in place to ensure TH homeostasis against transient perturbations, which involves three major TH binding/distributor proteins (THBPs/THDPs) in the blood: thyroxine-binding globulin (TBG), transthyretin (TTR), and albumin (ALB). Overall, greater than 99.97% of plasma T4 is bound to THBPs, while about 99.7% of plasma T3 is bound to THBPs (7–9). The majority of the THBP molecules exist unoccupied by THs (9, 10). The unique concentration ratios of T4, T3, THBPs, and their complexes (TH-THBPs) generate a strong buffer system where transient changes of free T4 and T3 can be quickly compensated by the absorption or release of T4 and T3 by the THBPs. With high amounts of THs stored in THBPs, TH-THBPs also act as a reservoir in the event of short-term deficiency of THs. Moreover, THBPs are required to maintain a uniform delivery of circulating THs throughout the perfused tissues (11, 12).

A variety of environmental chemicals have structures similar to THs and can interfere with TH binding to THBPs. These endocrine-disrupting chemicals (EDCs) include polychlorinated biphenyls (PCBs), dibenzo-p-dioxins, dibenzofurans, and brominated flame retardants such as polybrominated diphenyl ethers (PBDEs), pentabromophenol (PBP), and tetrabromobisphenol A (TBBPA) (13–18). Some of them such as certain PBDEs even have higher binding affinities for TBG or TTR than the endogenous hormones (18, 19). It has been suspected that these EDCs can displace THs from THBPs, causing an increase in free plasma THs (14, 20). However, given the fact that the predominant fractions of THBPs are unoccupied by THs (thus available for EDC binding) and the cross-buffering effect between THBPs, it is unlikely that the free THs will be altered tangibly by THBP-binding EDCs. The effect size of these thyroid disruptors is not well characterized and a significant knowledge gap exists.

Mathematical modeling has been used to understand the various aspects of the thyroid system and its responses to dietary, therapeutic, and environmental perturbations (21–29). A physiologically based kinetic (PBK) model of THs containing sufficient details such as all the THBP binding events will help to provide critical insights to the EDC effects. Existing TH kinetic models are insufficient for our purpose here because none of them have all the necessary model components to address these questions (21–24, 26, 27, 29). Here we constructed a minimally sufficient human PBK model of THs and used the model to investigate the effects of THBP-binding EDCs. We started by formulating a multicompartment model to include TH production, distribution, and metabolism. The model contains the thyroid, liver, and rest-of-body compartments with explicit descriptions of the reversible binding kinetics between THs and THBPs in the blood compartments. The model was parameterized and validated extensively based on literature data. Here we reported several novel findings on TH kinetics and the effects of EDCs.

2 Methods2.1 Construction of the human PBK model2.1.1 Model structure

The PBK model contains four overall compartments: Body Blood, Thyroid, Liver, and Rest-of-Body (RB) (Figure 1A). The Liver, RB, and Thyroid compartments are further divided into tissue blood (vasculature) and tissue proper (extra-vasculature) subcompartments. Since we are primarily concerned with the steady-state TH behaviors in the present study, the negative feedback regulation of TSH by THs is not considered. The potential limitations in this regard are addressed in Discussion.

Figure 1 Schematic illustrations of the human PBK model for T4 and T3. (A) The overall structure of the PBK model. RB: rest-of-body. (B) Reversible binding of T4 and T3 with TBG, TTR and ALB (green-shaded area), and reversible binding of an EDC, X, with TBG or TTR (pink-shaded area). These binding events occur in the plasma of all blood compartments.

2.1.2 Kinetic processes

In each blood (sub)compartments, there are 11 molecular species (state variables): free T4 (fT4), T4-bound TBG (T4TBG), T4-bound TTR (T4TTR), T4-bound ALB (T4ALB), free T3 (fT3), T3-bound TBG (T3TBG), T3-bound TTR (T3TTR), T3-bound ALB (T3ALB), free TBG (TBG), free TTR (TTR), and free ALB (ALB) (Figure 1B, green-shaded area). Each blood (sub)compartment is assumed to be well-mixed and has its own set of ordinary differential equations (ODEs) tracking the (sub)compartment-specific rates of change of these variables. The interaction between each TH and each THBP species follows the law of mass action, which involves a reversible binding process: second-order association and first-order dissociation. For each specific interaction, its binding kinetic constants (i.e., the association rate constant and dissociation rate constant) are assumed the same across all blood (sub)compartments.

In the Body Blood compartment, in addition to the binding events between THs and THBPs, the venous blood inflows coming back from the Thyroid, Liver, and RB compartments and arterial blood outflow are also described for each of the 11 molecular species. Similarly, these species enter each tissue blood subcompartment through the corresponding arterial blood inflow and exit through the venous blood outflow. In Liver and RB, fT4 and fT3 are able to move bidirectionally as first-order processes between the tissue blood and tissue proper subcompartments (Figure 1A: k21, k23, k25, and k27 for influx; k28, k29, k30, and k31 for efflux). Using first-order processes to describe these exchanges is well justified because although they are membrane transporter-mediated, the binding affinities (Km) of T4 and T3 for these transporters are normally in the high nM to µM range (30, 31), which are hundreds of times higher than the free T4 and T3 concentrations in the plasma at least (32, 33). The Thyroid compartment is simplified by describing zero-order, unidirectional releases of fT4 and fT3 into Thyroid blood (Figure 1A: k20 and k22 respectively) while the concentrations of T4 and T3 in Thyroid tissue per se are not modeled.

Within the tissue proper of the Liver and RB compartments, T4 and T3 partition into free and bound forms which are assumed to be in fast equilibrium, such that their free fractions are described by constants (FuT4 and FuT3). Only free T4 and free T3 in the tissue proper are metabolized and cleared according to first-order processes. Conversion of T4 into T3 by DIOs is described by the k24 (RB) and k26 (Liver) processes (Figure 1A). The metabolism of T4 into other metabolites, including reverse T3, glucuronidated and sulfated T4, which are not explicitly tracked in the model, is described by the k32 (RB) and k34 (Liver) processes. T3 is cleared by the k33 (RB) and k35 (Liver) processes.

After the completion of the model construction, a mass balance check was performed for T4 and T3 by terminating their productions and tracking the amounts of T4 and T3 in various compartments and metabolized over time (Figure S1). The sum was then compared with the steady-state amounts before the productions were turned off. As indicated in Figure S1, both T4 and T3 in the model are mass balanced.

To investigate the phases of plasma T4 and T3 clearance and their half-lives, TH tracer experiments were simulated, for which T4 tracer or T3 tracer are coded as separate state variables that have identical interactions and properties as the respective endogenous counterparts.

2.2 Modeling the effect of EDCs

To examine the effects of THBP-binding EDCs, a chemical X is introduced in the PBK model. Reversible binding of X to either TBG or TTR is considered (Figure 1B, pink-shaded area). In each blood (sub)compartments, an additional ODE for free X and two additional ODEs for the X-bound TBG (XTBG) or X-bound TTR (XTTR) complexes are included. For simplicity, X is introduced to the model through a zero-order production process in the Body Blood compartment, resembling some environmental exposures not explicitly considered here, and X is cleared through a first-order process also in the Body Blood compartment.

2.3 Sensitivity analysis

Local sensitivity analysis was performed by increasing or decreasing the relevant parameters by 1% from the default values and allowing the system to reach steady state. The averaged ratio of the percentage change of T4 or T3 in Liver or RB to the percentage change in the parameter was calculated as the sensitivity coefficient.

2.4 Model parameters, equations, and simulation tools

An average human individual of 75 kg body weight is simulated. Model parameter values and details of the source references, justifications, and mathematical calculations from which the parameter values were derived are provided in Table S1-S3 and their footnotes. The ODEs and algebraic equations are provided in Table S4. The models were constructed in MATLAB R2019a (MathWorks Natick, Massachusetts, USA), and ode15s was used to numerically solve the ODEs. All MATLAB code is available at https://github.com/pulsatility/2023-TH-PBK-Model.git.

3 Results3.1 Basic steady-state characteristics of the PBK model3.1.1 TH concentrations and distribution in blood

Steady-state TH concentrations and percentage distribution among free and THBP-bound THs in the Body Blood compartment are provided in Table 1. The plasma concentration of fT4 is 15 pM, which accounts for 0.016% of Total T4 in this compartment. T4TBG is at 75%, accounting for the majority of Total T4, followed by T4TTR at 17% and T4ALB at 7.8%. The plasma concentration of fT3 is 5 pM, which accounts for 0.3% of Total T3 in the Body Blood compartment. T3TBG is at 75%, accounting for the majority of Total T3, followed by T3ALB at nearly 20% and T3TTR at nearly 5%. The TH concentrations in the Liver blood and RB blood compartments are slightly lower but nearly indistinguishable from those in Body Blood (simulation results not shown). The fT4 and fT3 concentrations in the Thyroid blood reach 17.4 and 5.9 pM respectively, which are 16 and 18.6% higher than in Body Blood, but the THBP-bound THs are only marginally higher.

Table 1 Plasma TH concentrations and percentage distribution in Body Blood compartment.

3.1.2 THBP saturation

In the Body Blood compartment, TBG is 20.29% saturated, which is contributed predominantly by T4, accounting for 19.93%, while T3 only accounts for the remaining 0.36% (Table 2). The values are comparable to the 18.4% or 20% TBG saturation by T4 in the literature (9, 10). TTR is only 0.3% saturated by THs, also consistent with the 0.16% or 0.5% values in the literature (9, 10). TH-bound TTR is also dominated by T4, which accounts for 0.2991%, while T3 only accounts for 0.0015%. ALB is barely saturated by THs, and at 0.0012%, it is comparable to the estimated 0.0016% in the literature (10). TH-bound ALB is dominated by T4, which accounts for 0.0011% and T3 only accounts for 0.00005%. The percentage saturations of THBPs in the Liver blood and RB blood are slightly lower but nearly identical to those in Body Blood (simulation results not shown). The saturations in the Thyroid blood are slightly higher than in Body Blood because of the binding of THs newly released from the Thyroid tissue (Table 2).

Table 2 Percentage saturation of THBPs by T4 and T3.

3.1.3 TH concentrations and abundance in tissue compartments

The steady-state concentrations and amounts of THs in the tissue compartments are listed in Table 3. Total T4 concentrations in the Liver tissue and RB tissue are 2.26E5 and 9.89E3 pM, respectively. These concentrations are 2.42 and 0.11-fold of plasma total T4 concentration, respectively. This Liver tissue to Plasma concentration ratio (or partition coefficient) is consistent with the near 3:1 ratio estimated in earlier studies (34). T4 amounts in the entire blood, Liver tissue, and RB tissue are 2.968E5, 3.878E5, and 6.097E5 pmol, accounting for 23%, 30%, and 47% of all extrathyroidal T4, respectively. The T4 amount in the Liver tissue is thus 1.3 times higher than that in the blood, which is consistent with the median value reported in (35). In comparison, total T3 concentrations in the Liver tissue and RB tissue are 3.51E3 and 790 pM, which are 2.07 and 0.46-fold of total T3 plasma concentration, respectively. T3 in the entire blood, Liver tissue, and RB tissue are 5.397E3, 6.018E3, and 4.867E4 pmol, accounting for 9%, 10%, and 81% of all extrathyroidal T3, which is consistent with the estimation that over 90% of extrathyroidal T3 is in extravascular compartments (36). The percentages of T4 and T3 amounts in the Liver tissue, i.e., 30% and 10% respectively, are comparable to what were reported in the literature (24, 34, 35), as summarized in Table S2.

Table 3 Concentrations, amounts and percentage distribution of THs in tissue compartments.

3.1.4 TH production and metabolism in tissue compartments

The rates of steady-state TH production and metabolism in different tissue compartments are presented in Table 4. For an average individual of 75 kg body weight in our model, T4 is produced by the Thyroid at 1.6 pmol/S (equivalent to 107.39 µg/day), which is comparable to the values reported (23, 37, 38), as summarized in detail in Table S1. At this production rate, the daily turnover of T4 in all extrathyroidal tissues is about 10.68%. Nearly 35% of T4 is eliminated in Liver and 65% in RB. T3 is produced by the Thyroid at 0.1143 pmol/S (6.43 µg/day), accounting for 22% of total T3 production in the whole body, which is close to the 23.7% estimated previously (38). The remaining T3 is produced by T4 conversion in peripheral tissues, where Liver contributes 27% and RB 51%. This nearly 1:2 ratio is close to the 15:29 ratio estimated for T3 production in liver and muscle in euthyroid condition (39). Therefore, RB contributes to nearly 2/3 of T3 produced extrathyroidally. The daily turnover of T3 in all extrathyroidal tissues is about 74%. Peripheral T4-to-T3 conversion accounts for 25% of overall T4 metabolism, among which 8.73% and 16.28% occur in Liver and RB respectively. The Liver metabolizes 40% and RB metabolizes 60% of T3 produced overall.

Table 4 TH production and metabolism rates.

3.1.5 TH blood supply and blood-tissue exchange

The rates of the steady-state TH blood supply and influx and efflux in the tissue compartments are presented in Table 5. For the Liver compartment, T4 is supplied through the blood at a rate of 1.33E3 pmol/S for Total T4 and 0.21 pmol/S for fT4. For the exchange of T4 between Liver blood and Liver tissue, the influx is 63.54 pmol/S and efflux is 62.98 pmol/S, which are much greater than the blood supply rate of fT4 and only a small fraction of the blood supply rate of Total T4. T4 efflux is 99.12% of T4 influx, indicating that only a tiny fraction of T4 coming into the Liver tissue is metabolized and the vast majority returns to blood. T3 is supplied to the Liver compartment through the blood at a rate of 24.12 pmol/S for Total T3 and 0.071 pmol/S for fT3. T3 enters the Liver tissue at a rate of 6.97 pmol/S and exits at 6.91 pmol/S. Similar to the case of T4, these rates are much greater than the blood supply rate of fT3. T3 efflux is 99.05% of T3 influx, thus only about 1% (0.0664 pmol/S) coming into the Liver tissue is metabolized. The simulation results presented here are consistent with the estimation that >99% of T4 and T3 that enter the liver return to blood unmetabolized (7). T3 is produced from T4 locally in the Liver tissue at 0.14 pmol/s, and therefore, local T3 production contributes about 68% of T3 in the Liver tissue while uptake from the blood contributes the remaining 32%.

Table 5 TH blood supply and exit, blood-tissue flux, tissue production and metabolism.

Despite a much larger volume of the RB compartment than the Liver compartment, T4 influx and efflux in RB are much slower, at 14.99 and 13.95 pmol/S respectively, which are about 22%-23% of the corresponding fluxes in the Liver. These values are about 22-23 times higher than the RB blood supply rate of fT4. T4 efflux is 93% of T4 influx, indicating only a small fraction (7%) of T4 coming into the RB tissue is metabolized. Similar to the case of T4, T3 influx and efflux in RB are much slower that in the Liver, at 0.75 and 0.70 pmol/S respectively, which are about 10% of the corresponding fluxes in the Liver. The blood-tissue exchange rates are about 3.4 times higher than the blood supply rate of fT3. T3 efflux is 94% of the T3 influx, thus 6% (0.048 pmol/S) of T3 coming into RB is metabolized. T3 is produced locally at 0.26 pmol/S in RB, which is > 5 times higher than T3 transported from the blood into RB and is about 35% of the T3 influx rate. The overall metabolism rate of T3 in RB is 41% of the T3 influx rate.

The high blood-tissue exchange rates relative to the local production and/or metabolism rates, particularly in the Liver compartment, indicate that at the steady state the tissue concentrations of THs are largely determined by the quasi-equilibrium with plasma free THs. This suggests that the steady-state tissue concentrations of T4 and T3 would be relatively insensitive to perturbations of the local metabolism as long as the plasma fT4 and fT3 concentrations stay constant, e.g., through the HPT feedback regulation. Sensitivity analysis confirmed such prediction, as all sensitivity coefficients are much lower than unity, particularly in Liver (Table S5). These results suggest an intrinsic mechanism for local T4 and T3 homeostasis.

3.2 Transient TH tissue distribution kinetics and half-lives3.2.1 Distribution phases and half-lives

With the basic steady-state characteristics of the model largely in line with the literature, we next simulated TH tracer experiments to examine the phases of plasma and tissue T4 and T3 transients and their half-lives so derived. With a T4 tracer 3 times the concentration of fT4 introduced to the Body Blood compartment (45 pM), T4 tracer is cleared in 4 distinct phases (Figures 2A, B). The 1st phase is the rapid distribution of the T4 tracer from the Body Blood to the Liver blood and RB blood compartments, which is completed in less than a minute. This step dilutes plasma T4 tracer in the Body Blood from 45 pM to 30 pM by the approximate 2:1 plasma volume ratio between Body Blood and Liver blood+RB blood. The 2nd phase is the distribution to the Liver tissue compartment, as evidenced by the peaking of T4 tracer in Liver tissue at about 2.6 h or 0.11 day (Figure 2C). The 3rd phase is the distribution to the RB tissue compartment, as evidenced by the peaking of the T4 tracer in RB tissue at 22.66 h or 0.94 day (Figure 2D). Lastly, the 4th phase is the metabolism phase. The half-life of plasma T4 tracer, calculated based on the slope of the 4th phase, is about 6.43 days.

Figure 2 Simulation of TH tracer experiments to determine distribution phases and apparent half-lives of T4 and T3 in the PBK model. For the T4 tracer simulation, 45 pM of fT4 tracer is added to the Body Blood compartment at time 0. (A) Plasma concentration of total T4 tracer in Body Blood over time. (B) Magnified view of (A) to highlight the four T4 phases as indicated. (C-D) T4 tracer concentrations in Liver tissue and RB tissue respectively. For the T3 tracer simulation, 5 pM of fT3 tracer is added to the Body Blood compartment at time 0. (E) Plasma concentration of total T3 tracer in Body Blood over time. (F) Magnified view of (E) to highlight the four T3 phases as indicated. (G-H) T3 tracer concentrations in Liver tissue and RB tissue compartments respectively.

For the simulated T3 tracer experiment, T3 tracer is also cleared in 4 distinct phases (Figures 2E, F), with the 1st phase being the distribution between blood compartments, the 2nd phase the distribution to Liver tissue, completed at 0.5 h (Figure 2G), the 3rd phase the distribution to RB tissue, completed at 7.8 h (Figure 2H), and the 4th phase the metabolism phase. The half-life of plasma T3 tracer, calculated based on the slope of the 4th phase, is about 23.56 hours.

3.2.2 Effects of blood flow rate

Several factors can interplay to determine how fast T4 and T3 are taken up by tissues, including the blood flow rate, the transporter-mediated exchange rates between the tissue proper and tissue blood, and the availability of plasma THs owing to THBP binding. Among these factors, we first explored the blood flow rate by varying the cardiac output QC. Increasing QC accelerates the 1st phase, i.e., the distribution of T4 tracer from Body Blood (Figure S2A) into Liver blood (Figure S2B) and RB blood (Figure S2C), while the opposite occurs when decreasing QC. While the blood flow rate affects the duration of the 1st phase of plasma T4 tracer clearance, the remaining 3 phases are largely unaffected (Figures S2D–F). Similar results were obtained for T3 tracer simulations with varying QC (Figure S3).

3.2.3 Effects of blood-tissue exchange rates

We next explored the exchange rates between the blood and tissues. When both k25 and k30, the two parameters governing T4 exchange in Liver, are doubled simultaneously, which leaves the basal T4 level in Liver tissue largely unchanged, the T4 tracer peak in Liver tissue becomes higher and advanced from 2.6 h (0.11 day) to 1.73 h (0.072 day) (Figure S4B). When k25 and k30 are halved simultaneously, the T4 tracer peak in Liver tissue becomes lower and delayed to 4.07 h (0.17 day). The peak time in RB tissue is unaffected except that the rising phase becomes either slightly slower or faster (Figure S4C), which results from the altered T4 tracer distribution into Liver tissue. Therefore, the T4 influx and efflux rate constants in Liver play an important role in determining how fast T4 is taken up by the Liver tissue, which in turn determines the duration of the 2nd phase of T4 tracer clearance in Body Blood (Figure S4A).

Similarly, when k21 and k28 for T4 influx and efflux in RB are doubled simultaneously, which also leaves the basal T4 level in RB tissue largely unchanged, the T4 tracer peak in RB tissue becomes higher and advanced from 22.65 h (0.94 day) to 13.68 h (0.57 day) (Figure S4F). When k21 and k28 are halved simultaneously, the T4 tracer peak in RB tissue becomes lower and delayed to 36 h (1.5 day). As a consequence, the T4 peak in Liver tissue is also affected to some extent, which becomes lower but advanced to 2.22 h (0.092 day) or higher but delayed to 3.04 h (0.13 day) as k21 and k28 are doubled or halved, respectively (Figure S4E). Therefore, the T4 influx and efflux rate constants in RB play an important role in determining how fast T4 is distributed into the RB tissue, which in turn determines mainly the duration of the 3nd phase of T4 tracer clearance in Body Blood (Figure S4D).

While affecting the T4 distribution phases, varying neither k25 and k30 nor k21 and k28 has a tangibly impact on the half-life of T4. Similar effects on the timing and duration of the 2nd and 3rd phases of plasma T3 clearance can be also observed for T3 tracer simulations when the blood-tissue exchange rate constants, k27 and k31 for Liver and k23 and k29 for RB, are varied respectively (Figure S5). Taken together, these results suggest that how fast the membrane transporters can move THs in and out of tissues play a key role in transient TH distribution kinetics.

3.2.4 Effects of THBP availability

We next explored how the availability of THBPs can affect TH clearance and half-lives. As TBG, TTR, and ALB are sequentially removed from the system, the 2nd phase and, to a lesser extent, the 3rd phases of T4 tracer clearance are progressively advanced and the half-life becomes shorter (Figure 3). When only TBG is removed (solid blue line), the Liver peak of T4 tracer (Figure 3B) is advanced from 2.6 h (0.11 day) to 1.13 h (0.047 day) and reaches a higher level, the RB peak (Figure 3C) is advanced from 22.66 h (0.94 day) to 18.96 h (0.79 day) also reaching a higher level, and the T4 half-life is shortened from 6.43 day to 5.35 day (Figure 3A). Compared with TBG, removing TTR or ALB has a similar but much minor effect (Figures 3A–C, solid orange and green lines respectively). As any two of the three THBPs are removed from the system, similar results occur, with the removal of TBG+TTR combination producing the largest effect (Figures 3D–F, dashed green line). Therefore, the effect size is positively correlated to the total abundance of T4-THBPs removed. The advancement of the Liver and RB distribution phases of blood T4 tracer clearance occurs because as THBPs are removed, fewer T4 tracer molecules are bound and more free T4 tracer becomes available in the blood; as a result, T4 tracer translocates into Liver tissue and RB tissue faster, reaching a higher peak and equilibrating with blood sooner. The shortening of the T4 half-life is because as the Total T4 in the blood is reduced with some THBPs removed, the relative abundance of T4 in the tissues becomes higher; since T4 is metabolized only in the tissues not in the blood, the overall half-life becomes shorter. In the extreme case when all THBPs are absent, the quasi-equilibrium between blood and tissues is achieved extremely fast, in only a few minutes (Figures 3D–F, dashed purple line). This suggests that the process of free T4 transportation into the tissues per se is not limiting and it becomes limiting only when plasma T4 is considerably stored in THBPs. In the absence of all three THBPs, the half-life of plasma T4 tracer decreases to 4.83 days.

Figure 3 Effects of THBP availability on T4 distribution kinetics and half-life. For the T4 tracer simulation, 45 pM of fT4 tracer is added to the Body Blood compartment at time 0. Concentrations of T4 tracer in Body Blood (A), Liver tissue (B), and RB tissue (C) over time with only one THBP removed as indicated in (A). Concentrations of T4 tracer in Body Blood (D), Liver tissue (E), and RB tissue (F) over time with two or all three THBPs removed as indicated in (D). T4 half-life is indicated for each condition.

A similar trend of phase advancement and slightly shortened half-life also occurs for T3 tracer as each THBP species is removed from the system (Figure 4), with the counterintuitive exception that in the absence of all three THBPs, the T3 half-life is lengthened to about 27 hours. This can be explained by the vastly different rate constants of T3 metabolism in the RB tissue vs. Liver tissue as detailed in Figure S6.

Figure 4 Effects of THBP availability on T3 distribution kinetics and half-life. For the T3 tracer simulation, 5 pM of fT3 tracer is added to the Body Blood compartment at time 0. Concentrations of T3 tracer in Body Blood (A), Liver tissue (B), and RB tissue (C) over time with only one THBP removed as indicated in (A). Concentrations of T3 tracer in Body Blood (D), Liver tissue (E), and RB tissue (F) over time with two or all three THBPs removed as indicated in (D). T3 half-life is indicated for each condition.

3.2.5 Effects of TH-THBP binding rates

Given the effects of THBP availability in the transient TH distribution kinetics as revealed above, we next explored the effects of how fast THs can be released by THBPs. This is done by varying the pair of association and dissociation rate constants of TH-THBP binding simultaneously by the same fold, which would leave the binding affinities unchanged. When the pairs of k1 and k2, k3 and k4, and k5 and k6 are all increased by 10-fold simultaneously, the T4 tracer Liver peak is slightly advanced from 2.6 h (0.11 day) to 2.4 h (0.1 day), while a 10-fold decrease delays the peak to 3.49 h (0.15 day) and also lowers the peak level (Figure S8B). Correspondingly the 2nd phase of the plasma T4 tracer clearance is also advanced or delayed respectively (Figure S8A). These parameter variations do not affect the T4 trace RB peak (Figure S8C). These results suggest that how fast T4 can be released by THBPs also play a role, albeit minor, in T4 Liver distribution kinetics. For the T3 tracer simulation, varying both the association and dissociation rate constants for T3 binding to THBPs, including the pairs of k7 and k8, k9 and k10, and k11 and k12, simultaneously by the same fold, has a similar but even smaller effect than the case of T4 (Figures S8D–S8F).

3.3 Effects of THBP-binding EDCs

A major concern with EDCs targeting the thyroid system is that some of them are structurally similar to the two-ringed T4 and T3 molecules and may compete with THs for target proteins including the THBPs. If an EDC is able to avidly bind to a THBP species, the THs bound to that THBP may be displaced and liberated into the plasma, causing an increase in free THs, at least transiently. To explore these situations, we assume a hypothetical reference compound X that has the same binding affinities and rate constants as T4 for TBG or TTR.

Our simulations demonstrated that for continuous, long-term exposure where plasma X reaches a steady state, its effect on the steady-state fT4 and fT3 levels is negligible regardless of the X concentrations (simulation results not shown). When X is able to bind to TBG, the steady-state plasma XTBG level increases in a typical Michaelian manner as free X increases (Figure 5A). As a reference point, at a free X concentration that is the same as basal fT4, i.e., 15 pM, T4TBG decreases by 16.57% (Figure 5B, yellow line) while Total T4 decreases by 12.4% (Figure 5B, green line). T3TBG decreases by 16.93% (Figure 5C, yellow line) while Total T3 decreases by 12.63% (Figure 5C, green line). For T4TBG and Total T4 to decrease by 10% from their basal values, the required concentration of free X is 8.5 and 11.56 pM respectively. For T4TBG and Total T4 to decrease to 50% of their basal values, the concentration of free X must increase to 75 and 150 pM respectively. Despite these wide changes in T4TBG and Total T4, fT4 barely changes, and neither do T4TTR, T4ALB, and T4 in Liver tissue and RB tissue (simulation results not shown). The free X concentrations causing 10% and 50% drops in T3TBG and Total T3 are similar to the case of T4. fT3, T3TTR, T3ALB, and T3 in Liver tissue and RB tissue only have negligible changes (simulation results not shown).

Figure 5 Steady-state plasma concentrations of TH-THBPs and X-THBPs when an EDC X competes for binding to either TBG (A-C) or TTR (D-F) with the same binding affinities and rate constants as T4. Concentration-response curves for XTBG (A), T4TBG and Total T4 (B), T3TBG and Total T3 (C). Concentration-response curves for XTTR (D), T4TTR and Total T4 (E), T3TTR and Total T3 (F). Dashed lines indicate effects when free X concentration equals fT4 concentration (15 nM).

A similar but much minor effect is observed when X binds to TTR only (Figures 5D–F). At 15 pM, X only causes steady-state T4TTR and Total T4 to decrease by a tiny 0.001%, and T3TTR and Total T3 by 0.0001%. The lessor effect is because the reserve capacity of TTR (percentage TTR unoccupied by THs) is much higher (about 99.7%) than that of TBG (about 80%) (Table 2). Therefore, at the same free X concentration, fewer TH molecules will need to be displaced from TTR than from TBG since there are so many more unoccupied TTR molecules available for X. For T4TTR to decrease to 90% and 50% of its basal value, the concentration of free X must increase to 550 and 5000 pM, respectively. T3TTR has basically the same percentage drops at these two X concentrations. Unlike the case when X can bind to TBG, here total T4 and total T3 never decrease to levels lower than 82.5% and 95% of their basal values, respectively. fT4, fT3, other bound forms, and tissue THs only have negligible changes (simulation results not shown). Taken together, continuous exposure to THBP-binding EDCs does not perturb the free TH homeostasis.

We next explored the situation of intermittent, daily exposure, where X can be quickly metabolized with a short half-life such that its plasma concentration fluctuates throughout the day. Our simulations demonstrate that when X has an overall half-life of 1 h, a daily 8-h exposure to TBG-binding X can cause large spikes of plasma free X (Figure 6A) and XTBG (Figure 6B), which are cleared nearly completely at the end of the day. When the peak level of free X reaches 15 pM, the perturbation is sufficient to cause 7% daily fluctuation of fT4 (Figure 6C) and 5.4% of fT3 (Figure 6J), 4.8% of T4 (Figure 6H) and 5.3% of T3 (Figure 6O) in Liver tissue, and 2% of T4 (Figure 6I) and 1% of T3 (Figure 6P) in RB tissue. T4-THBPs and T3-THBPs all exhibit daily fluctuations of various degrees (Figures 6D–F, G, K–N).

Figure 6 Time-course responses of X and TH species as indicated (A-P) to 8-h daily exposure to TBG-binding X that has a half-life of 1 h. The simulation of exposure to X begins from the basal steady state of the PBK model. The exposure level is set such that if the exposure is continuous for 24 h each day, it produces a steady-state plasma free X concentration of 15 pM. X binds to TBG with the same affinity and rate constants as T4.

When the daily exposure duration becomes shorter, fT4, fT3, and tissue T4 and T3 tend to spike up to higher peak levels and dip to higher trough levels, whereas when the daily exposure duration becomes longer, the opposite occurs (Figure 7). These changes are at first unexpected because with a longer exposure, more THs are expected to be displaced from their THBP-bound forms and thus free THs are expected to rise to higher levels. This counterintuitive phenomenon can be explained as follows. Because of the fast binding kinetics, the T4 and T3 displaced from TBG due to competitive X binding is initially absorbed by TTR and ALB and a quasi-steady state is quickly established in minutes, like in a tracer experiment. With a longer daily exposure duration, the displaced plasma T4 and T3 will have more time to be distributed from the blood to tissues and metabolized therein, allowing plasma fT4 and fT3 to settle back closer to their basal steady-state levels before the exposure is terminated for the day (Figures 7B, C). At the moment the exposure is terminated, X drops quickly (Figure 7A) and those TBG bound by X is thus released. These “new” free TBG molecules will mop up plasma fT4 and fT3, causing them to undershoot below the steady-state levels (Figures 7B, C). Since with a longer daily exposure duration plasma fT4 and fT3 are already settling at lower levels during the exposure, as reasoned above, the termination of exposure will cause fT4 and fT3 to drop to even lower levels. During the no-exposure period, fT4 and fT3 recover some but less so with no-exposure period being shorter. The above pattern repeats as the next-day exposure starts again. The tissue T4 and T3 changes lag behind fT4 and fT3 except for Liver T3, and they reach higher peak and trough levels with shorter exposure duration (Figures 7D–G). Exceptions occur when the exposure duration is either too short or too long. When the duration is too short (2 h or less, Figure 7 blue lines), the extra plasma T4 and T3 displaced by X do not have as much time to move into the tissues, thus the tissue peak levels tend to be lower. When the exposure duration is too long (>23 h, Figure 7 cyan lines), there is not enough time for tissue T4 and T3 to move back to the blood when the plasma fT4 and fT3 dip, resulting in a smaller undershoot in the tissues. With 24-h continuous exposure, all hormone species reach steady states without fluctuations (Figure 7 black lines).

Figure 7 Effects of different duration of daily exposure to TBG-binding X that has a half-life of 1 h on X and TH species as indicated (A-G). The model was run from the basal steady state for 52 days of exposure such that a steady-state range of fluctuation was achieved. The responses on day 51-52 (labeled as 0-48 h) are shown. The exposure level is set such that if the exposure is continuous for 24 h each day, it would produce a steady-state plasma free X concentration of 15 pM. X binds to TBG with the same affinity and rate constants as T4.

Increasing the half-life and thus the stability of X while keeping the intermittent, daily exposure duration the same decreases the amplitude of hormone fluctuations. With a half-life of 10 h, the exposure only yields steady-state fluctuations of 1.3% in fT4 and 1.0% in fT3, 1.3% and 1.1% in Liver T4 and T3 respectively, and 0.5% and 0.3% in RB T4 and T3 respectively (Figure S9). With a half-life of 100 h, the steady-state fluctuations decrease even further, producing changes of less than 0.1% for all TH species (Figure S10). Unlike the case of TBG, 8-hour intermittent, daily exposure to X that binds to TTR only produces negligible fluctuations even with a half-life of 1 h (Figure S11). The TTR-binding X induces daily fluctuations of around 0.03% in fT4 and 0.006% in fT3, 0.019% and 0.006% in Liver T4 and T3 respectively, and 0.008% and 0.0009% in RB T4 and T3 respectively. For the longer X half-lives of 10 and 100 h, the fluctuations are even smaller (Figures S12-S13).

4 Discussion

In the present study, we built a minimal PBK model that quantitatively recapitulates the production, distribution, and the metabolism of THs for an average human individual. The model is rigorously defined by utilizing and estimating parameters and measurements reported in a multitude of studies in the literature as detailed in Tables S1-S3 and their footnotes, including THBP binding affinities, association/dissociation rate constants, TH and THBP concentrations, plasma/tissue partitioning, and half-lives, etc. The model allows for the exploration of many features of TH kinetics and provides several novel biological insights into the homeostatic roles of THBPs and the impact of EDCs.

4.1 Partitioning of plasma THs among THBPs

For plasma T4, T4TBG is the predominant form (75%) followed by T4TTR and then T4ALB. For T3, T3TBG is still the most abundant (75%), but T3ALB is more abundant than T3TTR. The percentage distributions of plasma TH abundance captured by our model (Table 1) are consistent with the ranges reported in the literature as detailed in Table S2. The relative distributions of THs among these THBPs are determined by both the respective binding affinities of THs for and the abundances of THBPs. ALB is the most abundant THBP with a molar concentration > 100-fold higher than TTR which in turn is > 10-fold higher than TBG (9, 40). However, the binding affinities of T4 for these THBPs follow a reverse order with drastically different values (8, 9, 41–43), which results in an order of T4TBG>T4TTR>T4ALB for their abundance. Although the binding affinities of T3 for these THBPs also follow an order that is the converse to their abundances, the gap between the binding affinities of T3 for TTR and ALB is not large enough to compensate for the abundance difference between the two THBPs, resulting in a T3ALB concentration that is 4-fold higher than T3TTR. As a homotetramer, TTR has two binding sites for THs; however, there exists negative cooperativity. This results in a 100-fold difference in the binding affinity for the two sites (44), which justifies the use of a single site for the binding as done in our model here. Absent this negative cooperativity, the T3TTR abundance may double, but will still be lower than T3ALB.

In the present model, fT4 is 0.016% of plasma Total T4, which is slightly more than half of the 0.03% value often cited in the literature (7–9). 0.03% could be an overestimate of the free T4 fraction for a number of reasons. In our model, to increase fT4 fraction by nearly 2-fold, i.e., from 0.016% to 0.03%, it will require either (1) reducing the total THBP abundance by nearly half, (2) reducing the binding affinities of T4 for THBPs by nearly half such that plasma Total T4 will be halved, (3) doubling fT4 and reducing the binding affinities of T4 for THBPs by nearly half, or some combinations of these. Adjusting only TBG-related parameters may achieve such a result because of its dominance. All these manipulations will require dramatic departure of many parameters or hormone levels from the reference values. Moreover, measuring free T4 using the traditional immunoassays tends to overestimate its concentrations compared with using modern liquid chromatography with tandem mass spectrometry-based reference method (45). It has been reported by laboratories that free T4 fraction can range between 0.005-0.03% (46). Schussler provided an estimate of 0.02% (10), and the free T4 fraction calculated from the mean free T4 and total T4 in the thyroid profile data of the National Health and Nutrition Examination Survey (NHANES) 2007-2012 cycles is slightly over 0.01% (33). Therefore, 0.03% likely represents an upper bound of the free T4 fraction, and 0.016% in our model may be closer to the average value. In contrast, fT3 in the present model is 0.3% of plasma total T3, which is consistent with the literature 0.3%-0.4% values (7–9).

4.2 Influx, efflux, and tissue TH homeostasis

Although THs are highly lipophilic, they do not diffuse through cell membrane readily (31, 47). Earlier studies demonstrated that the transportation of THs across the cell membrane is saturable, indicating a transporter-mediated process (48, 49). Membrane transporters such as monocarboxylate transporters MCT8/10 and those in the organic anion transporter polypeptides (OATP) family are required to move THs across the cell membrane (31, 47). For both T4 and T3 in both the Liver and RB, the paired influx and efflux rates are comparable and high, and the difference, which is equal to the net influx rate of T4 or T3, is only a tiny fraction of either of the two rates. For the Liver compartment, our model is consistent with the estimation that >99% of T4 and T3 that enter the liver return to blood unmetabolized (7). With the flux rates much higher than the tissue metabolism rates, the free THs in the tissue are near equilibrium with the plasma free THs. This suggests that as long as the HPT feedback regulation is able to maintain plasma free TH concentrations relatively constant, alterations in the TH metabolism in the tissue would have minor or even negligible impact on their concentrations in the tissue, especially for T4 in the liver, as demonstrated with the sensitivity analysis (Table S5). In contrast, if the metabolism rate is a significant fraction of the unidirectional influx rate, then the tissue concentrations of THs will change markedly with the metabolism constants even though the plasma concentrations can remain unchanged. Therefore, this equilibrium mode of operation may provide some intrinsic robustness against intra-tissue perturbations. This is in addition to other known local feedback and feedforward mechanisms mediated by DIOs that regulate tissue T3 homeostasis (50–54).

4.3 TH tracer kinetics and half-lives4.3.1 TH tissue distribution phases

The T4 tracer simulation showed that the first 3 distribution phases are completed in about one day, which is consistent with experimental observations (23, 55). The 2nd phase reflects the distribution to the Liver tissue compartment, as evidenced by the peaking of T4 tracer in Liver tissue at about 2.6 h, which is similar to the observations in T4 tracer experiments (34). Likewise, the 3rd phase reflects the distribution to the RB tissue and is completed in almost a day which is on par with experimental observations (23). T3 tracer simulation exhibits a similar but faster kinetic profile which is also consistent with experimental observations (23, 38).

Several factors may interplay to determine how fast a T4 or T3 tracer moves into a tissue, including blood flow rate, tissue influx and efflux rates, plasma free TH fractions, and THBP binding kinetics. By themselves, the influx and efflux of THs between blood and tissues per se are not limiting. This is because in the absence of THBPs, the tracer in the blood can reach equilibrium with tissues in a few minutes (Figures 3, 4), far quicker than the hours it takes when THBPs are present. Examining the parameter values of the rate constants governing the influx and efflux reveals that they are high enough such that the blood-tissue equilibrium can be reached in seconds or less in the absence of other limiting factors. Without THBPs, the blood mixing or circulation rate is the limiting factor.

Interestingly enough, the influx and efflux become a limiting factor for TH tracer tissue distribution when THBPs are present. This is demonstrated by the high sensitivity of the Liver peak and RB peak time of TH tracers to changes of the influx and efflux rate constants (Figures S4, S5). When THBPs are present, because the binding kinetics is much faster than the influx/efflux kinetics, the vast majority of the tracer is first bound by THBPs. With plasma free tracer at very low level, it takes time for the tracer to be released by THBPs and move into the tissues. Increasing the influx and efflux rates will accelerate this process by clearing plasma free tracer more quickly, which will shift the association and dissociation balance of tracer-THBPs to favor unloading more tracer. Consistently, when fewer THBPs are available, less tracer will be bound, resulting in higher plasma free tracer level, which will accelerate the distribution into the tissues (Figures 3, 4). In a similar manner, because the free fraction of plasma T3 is almost 20 times higher than that of T4, the Liver peak of T3 tracer occurs much earlier than the Liver peak of T4 tracer. As far as the effect of the speed of tracer unloading by THBPs is concerned, increasing the binding rate constants from the default values does not advance the Liver and RB peaks of the tracer, indicating that the unloading is not rate-limiting in physiological conditions. However, decreasing the binding rate constants by 10-fold delays the Liver peak without affecting the RB peak (Figure S8), suggesting the speed of unloading becomes limiting only under severely altered conditions.

Taken together, the tissue influx and efflux rates are the most limiting steps in the TH tracer distribution kinetics. This is further supported by the following facts. The T4 influx rate constants k25 = 4.243 L/S for Liver and k21 = 1 L/S for RB. At Liver and RB blood flow rates of 0.0142 L/S and 0.0419 L/S respectively, the equivalent total T4 blood supply rate constants are 88 L/S and 260 L/S respectively. The influx:blood supply ratio is thus 0.048 and 0.0038 for Liver and RB respectively, which is much lower than unity. This indicates that T4 tissue uptake is a transporter-limited, not a blood flow-limited process, consistent with the conclusion above. Similarly, the influx:blood supply ratio for T3 is 0.29 and 0.011 for Liver and RB respectively, also consistent with the conclusion that T3 tissue uptake is transporter-limited.

4.3.2 TH half-lives

The half-life of T4, calculated based on the slope of the 4th decay phase in the Body Blood compartment, is 6.43 days, which is close to the 6.5 days as estimated using the T4 production rate and the total amount of T4 in the extrathyroidal whole body in the present model. In comparison, the simulated T3 tracer experiment predicts T3 half-life to be 23.6 hours, which is slightly longer than the expected 22.5 hours. While the difference is small, it is worth some further investigation. Examining the T3 tracer amounts in Liver tissue and RB tissue during the 4th, metabolism phase revealed that RB contains about 10 times more T3 tracer than Liver does (Figure S7), while at the basal steady state the RB : Liver ratio for the endogenous T3 amount is about 8:1 (Table 3). Since the T3 metabolism rate constant in RB, k33, is at least 5 times smaller than that in Liver, k35, a higher T3 tracer distribution compared to the endogenous T3 in RB indicates that the overall T3 tracer metabolism in the whole body would be slower than the endogenous T3, leading to a longer half-life of T3 tracer. The reason T3 tracer has a higher RB : Liver ratio than the endogenous T3 is due to the following. Assume the RB : Liver ratio of T3 tracer amount is 8:1 at some time point, i.e., the same as the endogenous T3. Because T3 tracer is eliminated in Liver with a faster rate constant than in RB, over time T3 tracer in Liver will drop faster than in RB percentage-wise. The initial 8:1 ratio is thus not sustainable without continuous supply of T3 tracer, unlike the endogenous T3 which is constantly produced locally to compensate for the faster loss in Liver. As a result, the RB : Liver ratio of the T3 tracer will decrease continuously until a higher ratio value is reached such that T3 tracer amounts in RB and Liver decline at the same rate percentage-wise (Figure S7). The reason that the T4 half-life predicted by the simulated T4 tracer experiment is largely in line with expectation is because the T4 metabolism rate constants in RB and Liver are comparable. In summary, the result suggests that clinically measured T3 half-life using radioisotope may be slightly overestimated.

As expected, by holding a significant fraction of THs in the blood where they are not metabolized, the presence of THBPs lengthens the half-lives of THs, which is demonstrated by the sequential removal of THBPs from the system. The only exception is when all THBPs are absent, the T3 half-life is lengthened considerably from 23.5 h to 27 h. Similar to the mechanism above, this lengthening occurs because of the drastically different metabolism rate constants in RB tissue vs. Liver tissue, which leads to a more predominant distribution of T3 tracer into RB tissue (Figure S6).

The half-lives of T3 and T4 are altered under pathological thyroid conditions. In general, the half-lives decrease in hyperthyroid patients, and increase in hypothyroid patients (37, 55, 56). As mentioned above, the local feedback and feedforward regulation of T4 and T3 metabolism in the tissues may be the main contributors to the altered half-lives of THs under hyper/hypothyroid conditions. In addition, hyperthyroid patients are expected to have higher TBG saturation, therefore the total T4 in the blood would not increase linearly as more T4 is produced by the thyroid gland. More T4 will tend to be fractioned into the tissues where it is metabolized, contributing to shortened half-lives. Insufficiency or lack of THBPs generally also leads to shorter half-lives of THs because of the shift of THs distribution to tissues as discussed above. Indeed, markedly shortened T4 half-lives have been observed clinically in individuals with TBG deficiency (34).

4.4 Effects of EDCs

Many EDCs are structurally similar to THs and can perturb various aspects of the thyroid system (20). Some of them are capable of competing with THs for binding to THBPs such as TTR and TBG with comparable or even higher binding affinities, including PCBs, PBDEs, PBPs, TBBPAs, dibenzo-p-dioxins, and dibenzofurans (13–18). The testing for THBP binding was usually done with in vitro T4 displacement assays, which often use high, saturating concentrations of T4 (15, 18). Under such condition, ne