Mathematical Modeling of the Pituitary–Thyroid Feedback Loop: Role of a TSH-T3-Shunt and Sensitivity Analysis

Despite significant progress in assay technology, diagnosis of functional thyroid disorders may still be a challenge, as illustrated by the vague upper limit of the reference range for serum thyrotropin (TSH). Diagnostical problems also apply to subjects affected by syndrome T, i.e., those 10% of hypothyroid patients who continue to suffer from poor quality of life despite normal TSH concentrations under substitution therapy with levothyroxine (L-T4). In this paper, we extend a mathematical model of the pituitary–thyroid feedback loop in order to improve the understanding of thyroid hormone homeostasis. In particular, we incorporate a TSH-T3-shunt inside the thyroid, whose existence has recently been demonstrated in several clinical studies. The resulting extended model shows good accordance with various clinical observations, such as a circadian rhythm in free peripheral triiodothyronine (FT3). Furthermore, we perform a sensitivity analysis of the derived model, revealing the dependence of TSH and hormone concentrations on different system parameters. The results have implications for clinical interpretation of thyroid tests, e.g., in the differential diagnosis of subclinical hypothyroidism.


INTRODUCTION
In recent years, the mathematical modeling of human thyroid hormone homeostasis via the hypothalamic-pituitary-thyroid feedback loop has received an increasing amount of attention. Starting from early phenomenological models, more precise models have been developed based on molecular and pharmacokinetic data, see, e.g., Ref. (1)(2)(3)(4)(5)(6) for recent surveys on existing modeling approaches. These mathematical models can give important insight into the functionality of the hypothalamic-pituitary-thyroid axis and can be used to simulate the dynamic behavior of thyroidal hormone concentrations under different (euthyroid and non-euthyroid) conditions, and sometimes also for clinical decision-making (7). Furthermore, in Ref. (8), a method is proposed to compute personalized euthyroid setpoints that can be used for individualized diagnosis and treatment of thyroid diseases. While this static model is appealing due to its simplicity (only two parameter values have to be estimated), it does not consider any dynamic phenomena in the HPT axis, which are, however, of great importance for a deepened understanding of the HPT axis and ultimately the development of personalized optimal medication strategies. Another drawback is the absence of any consideration of T 3 , which has been shown to be significant not only as a key actor in the hypothalamic-pituitary-thyroid feedback loop (4,9) but also in maintaining a good quality of life (5).
The main objective of this paper is an improved mathematical modeling of the HPT axis in order to obtain a more detailed understanding of the dynamic phenomena occurring in thyroid hormone homeostasis. In particular, as a first contribution, we extend the model originally developed in Ref. (1,2) in order to incorporate new insights obtained through several recent clinical studies. In particular, we incorporate a direct TSH-T 3 path inside the thyroid, accounting for the central T 3 production by the thyroid. Existence of such a TSH-T 3 -shunt was hypothesized and demonstrated in several experiments and clinical observations (10)(11)(12)(13)(14)(15). In Ref. (10), it was shown that L-T 4 -treated athyreotic patients exhibit decreased FT 3 concentrations despite normal free thyroxine (FT 4 ) levels, which would not be the case if peripheral FT 3 was mainly produced by deiodination of peripheral FT 4 . Furthermore, the sum activity of step-up deiodinases (G D ) is positively correlated with the TSH concentration (11) and with the thyroidal volume (12) and significantly decreases after thyroidectomy. These observations suggest that besides the peripheral T 4 /T 3 conversion, also TSH-stimulated deiodinases inside the thyroid contribute to the total T 3 production. In our work, we show that the extended model including such a TSH-T 3shunt is in good accordance with various clinical observations. For example, we show that the FT 3 concentration shows a clear circadian pattern, as was observed in vivo in Ref. (16). Notably, this is not the case in the previous model, which did not include the TSH-T 3 -shunt.
As a second main contribution of this paper, we perform a sensitivity analysis of the derived model. Loosely speaking, the (first-order) sensitivities are a measure for how "sensitive" certain system states (i.e., TSH or hormone concentrations) are with respect to changes in certain parameters (such as, e.g., the thyroid's secretory capacity G T ). These sensitivities reveal structural insight into the functionality of the hypothalamic-pituitary-thyroid axis and can provide explanations for certain clinical observations. For example, we show that the sensitivity of TSH with respect to G T is much higher for low values of G T (i.e., in hypothyroidism) than for high values of G T (i.e., in hyperthyroidism). This fact can be used to explain why in clinical practice, TSH concentrations may significantly vary beyond the upper limit of the reference range despite normal thyroid function.
The remainder of this paper is structured as follows. Section 2 presents the extended mathematical model and discusses the identification of the required (additional) parameters. In Section 3, we show simulation results of the derived model and discuss the observed properties (such as the existence of a circadian rhythm in FT 3 concentrations). A sensitivity analysis of TSH, FT 4 , and FT 3 concentrations with respect to different parameters is performed in Section 4. Finally, we conclude the paper in Section 5.

PRESENTATION OF THE EXTENDED MODEL AND PARAMETER IDENTIFICATION
As outlined above, several clinical observations have led to the hypothesis that a direct, TSH-stimulated path exists for T 3 production inside the thyroid, which we now incorporate into the mathematical model from Ref. (1,2). The extended model including this TSH-T 3 -shunt is illustrated in Figure 1, see Section S1 in the Supplementary Material for a mathematical description of the underlying differential equations. To this end, both intrathyroidal conversion of T 4 into T 3 via type 1 and 2 5 ′ -deiodinases as well as a direct synthesis of T 3 are modeled (see upper three blocks in the "Thyroid" block in Figure 1). Both mechanisms are stimulated by TSH and modeled via nonlinear Michaelis-Menten-Hill kinetics, see Section S1 in the Supplementary Material for further details.
Most of the parameters of the extended model can be taken from the model in Ref. (1,2), where the parameters have been estimated according to known physical quantities (such as the half-life period of certain substances, etc.) or have been identified using data measured in vivo. A detailed listing of these parameters can be found in the Tables S1-S3 in Supplementary Material. Some of the parameters were calibrated according to average population data, and hence the resulting model can be interpreted to be a functional model of some generic euthyroid subject. Clearly, personalized model identification would be highly valuable for individualized clinical decision-making and the development of personalized optimal medication strategies. For this, however, sufficient data such as individual dynamic trajectories of hormone concentrations would be needed to avoid overfitting. We note that, while the present report deals mainly with average population data, the observed phenomena are in good accordance with individual samples (9). For the extended model, the new parameters G T3 and k have to be determined. Also, the sum activity of the type 1 5 ′ -deiodinase, G D1 , has to be re-estimated. This is the case since the extended model considers the additional T 3 secretion inside the thyroid, while in the original model, G D1 was calibrated by only considering peripheral T 3 production, and hence G D1 had been estimated too high. In order to obtain the parameters G T3 and G D1 , a least squares estimation was performed, fitting the FT 3output of the presented model to FT 3 measurements of a clinical study. Although there is no unique solution in case that only single FT 3 measurements are available, it provides a set of optimal parameters, which could be further reduced to a unique solution if additional measurements were available (compare the detailed discussion below). In order to perform the least squares estimation, the equilibrium FT 3 level predicted by the extended model, in the following denoted by FT 3,eq , can be computed in dependence of the parameters by solving a cubic polynomial (see Section S1 in the Supplementary Material for a more detailed description). For this computation, we set TRH to a constant value (later, for the dynamic analysis TRH is varying in a sinusoidal fashion). This equilibrium value is then fitted in a least-squares sense to real measurement data resulting from 1,121 untreated patients of the clinical study in Ref. (11). In particular, this is FIGURE 1 | Block diagram of the thyrotropic feedback control loop with an additional TSH-T 3 -shunt, adapted from Ref. (1,2). Except for G T 3 , k, and G D1 , all parameters were adopted from the model in Ref. (1,2). The parameters G D1 and G T 3 were estimated to obtain an optimal (in a least squares sense) fit to measured in vivo FT 3 -concentrations. To this end, the value of k was normalized to 1 mU/l. achieved by minimizing the following cost function with respect to the parameters G T3 , k, and G D1 : Here, FT 3,i denotes the measured FT 3 -concentration of the i-th patient, and FT 3,eq (G T3 , k, G D1 ) is the equilibrium FT 3 -level predicted by the model depending on the parameters G T3 , k, and G D1 . The other parameters that are needed to compute FT 3,eq are adopted from Ref.
(1) (see Tables S1-S3 in the Supplementary Material). The optimal solution to the above optimization problem can be determined analytically and is given by FT 3,eq (G T3 , k, G D1 ) = FT 3 , where FT 3 is the mean value of the 1,121 FT 3 measurements. Using the derived formula for FT 3,eq (G T3 , k, G D1 ) (see Section S1 in the Supplementary Material), this results in different (infinitely many) optimal parameter combinations for G T3 , k, and G D1 . For example, normalizing k to 1 mU l (which will be used in the following), the optimal parameter combinations for G T3 and G D1 can be seen in Figure 2.
Different (optimal) parameter combinations for G T3 and G D1 result in different fractions of thyroidal and peripheral T 3 production. For example, G D1 = 22 nmol s and G T3 = 394 fmol s approximately lead to 80% T 3 production from peripheral conversion of FT 4 and approximately 20% T 3 production from intrathyroidal secretion, corresponding to the values suggested by Ref. (17,18). On the other hand, also, a higher or lower fraction of intrathyroidal T 3 production is possible, depending on the values of G T3 and G D1 . In particular, higher values of G T3 and lower values for G D1 result in a higher fraction of intrathyroidal T 3 production FIGURE 2 | Set of optimal (in a least-squares sense) parameters G T 3 and G D1 when normalizing the parameter k to 1 mU/l. Due to the affine dependence of FT 3, eq (G T3 , k, G D1 ) on G T 3 and G D1 , the set of optimal parameters is contained in a one-dimensional affine subspace of R 2 . and vice versa. For the dynamic simulation of the model and the sensitivity analysis in the following sections, we (mostly) use the values G D1 = 22 nmol s and G T3 = 394 fmol s , and we comment when certain results qualitatively change if other parameter values for G D1 and G T3 are used.
The above discussed non-uniqueness in the optimal parameter fit is due to the fact that the model is not fully identifiable given the measured data. Namely, FT 3 is the only hormone that is affected by the parameters G T3 , k, and G D1 , and we only have stationary measurements available. Furthermore, in the above estimation, we made the simplifying assumption that peripheral and thyroidal deiodinase activities (G D1 and G D2 ) are the same, which might in general not be the case. Identifying the corresponding parameters separately would result in a possibly better parameterized model, which is, however, again not possible given only the stationary FT 3 measurements. On the other hand, if we had additional data such as dynamic hormone concentration trajectories or additional measurements (e.g., intrathyroidal hormone concentrations), the above described non-uniqueness in the parameter estimation could be removed and also different parameter values for thyroidal and peripheral deiodinase activity could be identified, allowing for a more exact parameterization of the model. This would be an interesting topic for future research, however, such in vivo data are difficult to obtain and are typically not available. Moreover, the presented model does not consider membrane transport processes between thyroidal and peripheral tissue. Incorporating such processes by means of a compartment model would further increase the quality of our model, yet, this would yield additional parameters, which had to be identified. Nevertheless, as we will show in the following sections, the extended model with the parameters as identified in this section is a clear improvement compared to the previous model, allowing for a better reproduction and interpretation of various clinically observed phenomena.

DYNAMIC PROPERTIES OF THE EXTENDED MODEL
In the following, we simulate the extended model and analyze and interpret the obtained results. First, some simulation runs are carried out to illustrate the role of the TSH-T 3 -shunt in obtaining a circadian rhythm in the FT 3 -concentration. Afterwards, we investigate the delay between TSH and FT 3 , which has been observed in several clinical studies [e.g., Ref. (16,19)].

Dynamic Simulation
As detailed in the previous section, the intrathyroidal T 3 secretion is composed of two mechanisms, namely intrathyroidal conversion of T 4 into T 3 via type 1 and 2 5 ′ -deiodinases (upper middle and right block inside the thyroid in Figure 1) as well as a direct synthesis of T 3 (upper left block inside the thyroid in Figure 1). In the dynamic simulation using the parameters as identified in Section 2, the intrathyroidal contribution to the total T 3 secretion rate was composed as follows: Hence, with the parameters identified in Section 2, the main thyroidal source to T 3 -production is direct T 3 -synthesis via Michaelis-Menten-Hill kinetics, represented by the block "T3 Synthesis" in Figure 1. On the other hand, deiodination by type 2 5 ′ -deiodinases has a negligible effect only, since the sum activity of type 2 5 ′ -deiodinases is much smaller compared to that of type 1 5 ′ -deiodinases. In case that a different optimal combination of parameters G T3 and G D1 is used (compare Section 2), the above results change accordingly, i.e., a higher value of G D1 yields a higher contribution of the deiodination by type 1 5 ′ -deiodinases to the thyroidal T 3 -production. However, this also causes a change in the ratio between thyroidal and peripheral T 3 production, as discussed in Section 2. Figure 3 shows simulated FT 3 -plots, where we further investigated the effect of the TSH-T 3 -shunt on the dynamic behavior of FT 3 . 1 In particular, Figure 3A shows simulation results using the previous model from Ref.
(1) without the TSH-T 3 -shunt whereas in Figure 3B, the full TSH-T 3 -shunt as described in the previous section is included. For each of the two scenarios (i.e., for the corresponding models), we separately identified the (in a least-squares-sense) optimal parameter(s): G D1 for the model corresponding to Figure 3A and G D1 as well as G T3 for the model corresponding to Figure 3B. The exact values of these parameters for the different model configurations can be seen in Table S4 in Supplementary Material. We observe that the TSH-T 3 -shunt causes a clear circadian oscillation in the FT 3 concentration, which is not (or only very weakly) present without considering intrathyroidal T 3 secretion. Such a circadian rhythm in FT 3 concentration has been observed in vivo in several clinical studies [see, e.g., Ref. (16,19)], and hence our simulation results again support existence of the TSH-T 3shunt.
Quantitatively, the oscillation amplitude of the measured in vivo FT 3 concentration in Ref. (16) is approximately six times as big as the amplitude observed in the simulated model (see Figure 3B). This difference might be due to the assumptions we made for the identification in Section 2 (same values for G D1 , G D2 , K M1 , and K M2 inside the thyroid and the peripheral tissue). Namely, if thyroidal deiodination activity and/or G T3 were higher than computed in Section 2, without increasing the peripheral deiodination activity as well, we would obtain a larger oscillation amplitude in FT 3 concentration. Nevertheless, the fact that a clear circadian pattern arises in the simulations when including the TSH-T 3 -shunt into the model is a clear indicator supporting both its existence as well as the fact that thyroidal T 3 secretion is stimulated by TSH.

Delay of FT 3 w.r.t. TSH
The authors in Ref. (16) make the observation that in vivo FT 3measurements follow a clear circadian pattern, which is approximately 90 min delayed w.r.t. TSH; this number can also vary between different individuals (19). As already mentioned in the previous section, the FT 3 -level obtained by the model in Figure 1 including intrathyroidal T 3 -secretion shows a clear circadian pattern. In this section, we investigate how the delay between TSH and FT 3 in the presented model is influenced by this newly incorporated mechanism.
A dynamic simulation with the same setup as in Section 3.1 yields the following: when incorporating the TSH-T 3 -shunt into the model, FT 3 is delayed w.r.t. TSH by approximately 6 h, whereas the delay amounts to 13 h in the previous model, which did not incorporate this mechanism. These observed values can be explained as follows. The phase shift between FT 3 and TSH in our model mainly results from the first order lag elements α iω+β modeling peripheral T 3 and T 4 secretion (i.e., the ones with parameters α 31 , β 31 , and α T , β T , respectively in Figure 1). The phase shift of the output signal of such a first order lag element for a given sinusoidal input signal with frequency ω depends on the parameter β and is given as follows: In our case, ω = 2π T where T = 86,400 s is the circadian period of 1 day. The delay between the output and input signal is now computed by simply relating the phase shift to the period length T: delay = −phase· T 2π . For the given parameter values α 31 , β 31 , and α T , β T of T 3 -and T 4 -generation, respectively, we obtain a delay of approximately 5.5 and 6 h, respectively.
The above observed delay of FT 3 w.r.t. TSH can now be explained as follows. In the previous model not including the TSH-T 3 -shunt, the circadian oscillation has to pass through both first order lag elements for peripheral T 4 and T 3 production, resulting in a high delay w.r.t. TSH. On the other hand, the fraction of T 3 secreted inside the thyroid does not exhibit the delay caused by peripheral T 4 production and hence exhibits a much shorter delay w.r.t. TSH. Interestingly, the observed delay of total T 3 (approximately 6 h) mainly seems to be determined by the shorter one resulting from intrathyroidal T 3 production, although approximately 80% of the total T 3 -production results from peripheral FT 4 -deiodination and only 20% from intrathyroidal secretion. The reason for this is that as explained above, the circadian rhythm of FT 3 is mainly induced by intrathyroidal T 3 secretion. Namely, the ratio of the amplitude and the mean value equals 0.3% for the peripheral T 3 production rate PR(T 3 , peripheral) and 23% for the thyroidal T 3 production rate PR(T 3 , thyroid). Thus, the phase of FT 3 is almost solely characterized by the phase of thyroidal T 3 -production, and hence the delay of FT 3 w.r.t. TSH is determined by the phase shift of only one first order lag element when the shunt is included, compared to two without the shunt. To conclude, the inclusion of the TSH-T 3 -shunt into the model significantly reduces the delay of FT 3 w.r.t. TSH. While the absolute numbers are still too high compared to the observed in vivo delays (16,19), this is again a clear indicator for the existence of the TSH-T 3 -shunt.

SENSITIVITY ANALYSIS AND STATIONARY DEPENDENCIES
In this section, we perform a sensitivity analysis of the previously presented mathematical model of the hypothalamic-pituitary-thyroid feedback loop (see Figure 1). Sensitivity analysis is a tool for determining how a certain parameter influences the trajectories resulting from simulation of the model, i.e., from the solution of the underlying system of differential equations, and in particular, how "sensitive" these trajectories are with respect to certain parameter changes. In the following, we give a brief non-formal introduction to sensitivity analysis and refer to the Section S2 in Supplementary Material for a more complete and formal description.
To define sensitivities, consider the following vector-valued ordinary differential equation with parameter vector p: The first-order sensitivity function 2 is now defined as S(t) = ∂x(t,p) ∂p p=p 0 , where p 0 is some nominal (constant) parameter value. The sensitivity function S(t) is a time-dependent matrix with as many rows as the dimension of x and as many columns as the dimension of p. Under some assumptions (smoothness, existence of solutions, . . .), it can be shown that S satisfies the following differential equation, which is solved simultaneously with the state equation (3), see Ref. (20).
The initial sensitivity, i.e., S(t 0 ), is set to zero since the states' initial values are independent of the parameters. The above presented mathematical model of the hypothalamic-pituitary-thyroid feedback loop (see Figure 1) includes 36 parameters. With 5 states (pituitary TSH and T 3 as well as peripheral TSH, T 4 , and T 3 concentrations), this makes a total of 180 different sensitivity curves -for one specific nominal parameter configuration p 0 . In the following, we only analyze a few interesting curves to obtain some new insights. Of course, if desired, one could analogously analyze further sensitivities of other state and parameter pairs. In order to be able to employ the standard sensitivity analysis tools described above, the time delays in the hypothalamic-pituitary-thyroid (HPT) axis model are neglected.

Sensitivity of T 4 w.r.t. G T
We start by examining the sensitivity of peripheral T 4 with respect to the thyroid's secretory capacity G T . In Figure 4, several plots of the sensitivity ∂T 4 ∂G T (t) are shown with different G T -values, corresponding to different parameter values p 0 in equation (4). 3 Note that the sensitivity curves do not exhibit large variations over the day, i.e., only show a small circadian oscillation. It can be seen that different values of G T result in different sensitivities of T 4 with respect to G T . For example, a low value of G T (which can be seen as a simple modeling of hypothyroidism) causes an increase of the sensitivity, whereas a high G T value (which can be seen as a simple modeling of hyperthyroidism) causes a decrease of the sensitivity. This means that larger fluctuations in T 4 can be expected at the lower end of its euthyroid reference range (compare Section 4.3). This observation can compactly be expressed for a wide range of G T -values by investigating the stationary sensitivity (i.e., lim t→∞ ∂T 4 ∂G T (t)). Figure 5 shows it as a function of the parameter G T . A comparison between Figures 4 and 5 shows that the stationary sensitivity, indeed, is the limit of the sensitivity curve for t → ∞.

Sensitivity of TSH w.r.t. TRH
It is interesting to observe that the Ultra-Short-Feedback loop (i.e., the lower left part inside the pituitary in Figure 1, compare Ref. (2)) has a significant influence on the sensitivity of TSH w.r.t. TRH. Figure 6 shows the curves of this sensitivity for different values of S s . It can be seen that an increase in S s causes a decrease in the sensitivity.
In the considered HPT axis model, TRH is treated as a timedependent input that comes from the hypothalamus. A disturbance in the system could lead to a change in the TRH concentration arriving at the pituitary. Apparently, the Ultra-Short-Feedback increases the robustness of the TSH production w.r.t. changes in portal TRH. If the additional feedback is absent (i.e., S s = 0), the sensitivity is significantly higher than in the nominal case (S S = 100 l mU ).

Stationary Dependencies of TSH and T 4 on G T
In the following, the influence of the thyroid's secretory capacity G T on the equilibrium concentrations of TSH and T 4 is analyzed.
To this end, we solve the system's stationary equations (i.e., for t → ∞) for the different hormones and plot the resulting equilibrium hormone levels as functions of G T . The slopes of these functions are exactly the entries of the stationary sensitivity matrix lim t → ∞ S(t). For example, the stationary sensitivity lim t→∞ ∂T 4 ∂G T (t) for a given value of G T is equal to the derivative of the curve T 4 (G T ) w.r.t. G T , which we treat in the following.
The curves of T 4 and TSH depending on G T are shown in Figure 7. One can see that the parameter G T can be used as a measure of hypo-or hyperthyroidism (1). The equilibrium T 4concentration increases almost linearly with G T . Furthermore, we have high TSH-levels for low values of G T , i.e., in hypothyroidism, and vice versa. This is a well-known fact, which is usually used in clinical decision-making for the determination of subclinical thyroid diseases. Another interesting fact is that the magnitude of the sensitivity of TSH w.r.t. G T (which is the slope of Figure 7B) is high for low values of G T and vice versa. This means that TSH is much more sensitive to fluctuations in the thyroid's secretory capacity if G T is low. This fact can be used to interpret the following clinical observation. In practice, TSH concentrations may be misleading, especially, if located slightly above the vague upper limit of the reference range. A reason for this could be that as discussed above, TSH is much more sensitive to fluctuations in the thyroid's secretory capacity (e.g., due to different iodine supply and other influences) at the upper limit of its (euthyroid) reference range than at its lower limit.

Sensitivity of FT 3 w.r.t. G T
Finally, we perform a sensitivity analysis for FT 3 w.r.t. the parameter G T . Figure 8 shows the sensitivity curve, where in Figure 8A, the extended model including the TSH-T 3 -shunt was used, whereas Figure 8B uses the previous model without the shunt. It can be seen that the sensitivity of FT 3 w.r.t. G T decreases when the shunt is included. This is to be expected since in the extended model, a direct synthesis of T 3 (upper left block inside the thyroid in Figure 1) is included, which is independent of  G T , i.e., from the thyroid's secretory capacity for T 4 . Another interesting fact is that we can observe a small circadian rhythm in Figure 8A whereas the plot ( Figure 8B) seems not to be affected by this. This confirms the observations we made in Section 3.1, namely that incorporating intrathyroidal T 3 -secretion causes a circadian rhythm in FT 3 and hence also in the sensitivity w.r.t. G T . As for T 4 and G T , we can also analyze the stationary sensitivity lim t→∞  outcomes of a clinical study lead to the observation that the dependency of T 3 -generation on G T is lower than that predicted by a model, which does not include a TSH-T 3 -shunt. Indeed, comparing Figures 9A,B, one can see that the sensitivity of FT 3 w.r.t. G T significantly decreases when the shunt is incorporated into the model, i.e., T 3 production is less sensitive to fluctuations in the thyroid's secretory capacity if the shunt is included. As already mentioned above, this seems plausible since we now have a completely G T -independent path from TSH to FT 3 (upper left block inside the thyroid in Figure 1).

CONCLUSION AND OUTLOOK
In this work, a mathematical model of the hypothalamic-pituitary-thyroid feedback loop was extended to include TSH-stimulated intrathyroidal T 3 -secretion. The hypothesis of the existence of such a TSH-T 3 -shunt has been brought forward in various recent clinical studies. Our results show that the hypothesized mechanism can indeed explain various clinical findings. In particular, we have shown that intrathyroidal T 3 -secretion results in a clear circadian pattern of peripheral FT 3 , which has been observed in vivo in, e.g., Ref. (16,19), and which is not the case without the incorporation of such a TSH-T 3 -shunt. Also, a sensitivity analysis revealed that the sensitivity of peripheral FT 3 with respect to the thyroid's secretory capacity for T 4 is indeed lower when including intrathyroidal T 3 -secretion into the model, in accordance with the clinical study of Ref. (10). While the present report deals primarily with technical aspects of the thyroid pituitary feedback regulation, a better understanding of the underlying control system is of high clinical interest and relevance. Currently, clinical diagnosis and treatment of thyroid disease heavily relies on an indirect approach assessing the pituitary TSH response rather than circulating free thyroid hormones, FT 3 and FT 4 (21). The application is based on the underlying assumption that pituitary TSH in equilibrium at all times provides an accurate mirror image of the peripheral hormones. However, recent evidence has challenged this simplistic tenet suggesting that the HPT axis is a much more dynamic system than has been previously thought (5,22). In particular, the interrelationships between FT 3 , FT 4 , and TSH are less constantly fixed, rather conditional and contextualy adaptive (5,22). Mathematical modeling presented in this study confirms and advances the theoretical framework that is emerging from recent clinical studies. Given the high prevalence of subclinical thyroid disorders in the population, being as high as 10% in middle aged women, the epidemiological and therapeutic implications are substantial. From the performed sensitivity analysis in the present study, important insights into the functionality of the HPT axis have been obtained. These include the robustification of TSH production through the ultrashort feedback loop in the pituitary, as well as a possible explanation why in clinical practice, diagnosis of wrong subclinical hypothyroidism is much more common than diagnosis of wrong subclinical hyperthyroidism.
In particular, the upper reference limit for TSH has been a matter of fierce debate for a decade (23). According to our models, the issue appears to be more fundamentally rooted. This relates to a substantial error rate, depending on the statistical analytical technique used, in the conventional disease classification based solely on statistical TSH abnormality (24). The relative variability in TSH rises even further with higher TSH concentrations in subclinical hypothyroidism (25). Recent guidelines have raised the clinical threshold for therapeutic intervention in subclinical hypothyroidism to a TSH level of 10 mU/l, whereas the laboratorybased disease definition continues to rely on the upper reference limit of approx. 4 mU/l (21). A better understanding and refined mathematical expression of hypothalamic-pituitary regulation in allostatic reactions (22), in thyrotropic insufficiency (26), and in situations of imminent thyroid failure (9,27) as well as in their interactions may help reconcile this discrepancy that poses a considerable challenge to clinical decision-making.
Future work should focus on the further extension of the model to include currently unmodeled phenomena and mechanisms, such as, e.g., non-classical thyroid hormone signaling (28) and compartment models for the incorporation of membrane transport processes, which are increasingly understood as a regulatory element in their own right (29,30). We would also aim to define the steady-state more narrowly and precisely for individual subjects under different conditions in an attempt to reduce the high uncertainty surrounding TSH measurements at the upper limit of its reference range. In general, obtaining further insight into the overall functionality of the hypothalamic-pituitary-thyroid feedback loop and developing suitable and detailed enough mathematical models might eventually pave the way for designing optimal medication strategies for various non-euthyroid states of human hormone homeostasis.

AUTHOR CONTRIBUTIONS
JB and MM drafted the manuscript. JB performed the simulations using Matlab/Simulink. Figure 1 was designed by JD and modified by JB, all other figures were created by JB. The deidentified data used for the parameter identification was provided by RH. All authors read and approved the manuscript.