Mathematical modeling and simulation of thyroid homeostasis: Implications for the Allan-Herndon-Dudley syndrome

Introduction A mathematical model of the pituitary-thyroid feedback loop is extended to deepen the understanding of the Allan-Herndon-Dudley syndrome (AHDS). The AHDS is characterized by unusual thyroid hormone concentrations and a mutation in the SLC16A2 gene encoding for the monocarboxylate transporter 8 (MCT8). This mutation leads to a loss of thyroid hormone transport activity. One hypothesis to explain the unusual hormone concentrations of AHDS patients is that due to the loss of thyroid hormone transport activity, thyroxine (T 4) is partially retained in thyroid cells. Methods This hypothesis is investigated by extending a mathematical model of the pituitary-thyroid feedback loop to include a model of the net effects of membrane transporters such that the thyroid hormone transport activity can be considered. A nonlinear modeling approach based on the Michaelis-Menten kinetics and its linear approximation are employed to consider the membrane transporters. The unknown parameters are estimated through a constrained parameter optimization. Results In dynamic simulations, damaged membrane transporters result in a retention of T 4 in thyroid cells and ultimately in the unusual hormone concentrations of AHDS patients. The Michaelis-Menten modeling approach and its linear approximation lead to similar results. Discussion The results support the hypothesis that a partial retention of T 4 in thyroid cells represents one mechanism responsible for the unusual hormone concentrations of AHDS patients. Moreover, our results suggest that the retention of T 4 in thyroid cells could be the main reason for the unusual hormone concentrations of AHDS patients.

The entrance of T 4 in the blood stream is modeled as first order lag element, taking into account dilution (α T ) and clearance (β T ) of T 4 . Furthermore, a dead time (τ 0T ) is consider to account for diffusion processes. Additionally, the time constant of the first order lag element (τ 1T ) is explicitly shown.
Only a small fraction of T 4 is available as F T 4 , since most of it is bound to plasma binding proteins. This phenomenon is taken into account by considering the concentrations of thyroxine-binding globulin (T BG), transthyretin (T BP A), and their dissociation constants (K 41 /K 42 ) as displayed in the bottom right corner.
The peripheral production of T 3 can be seen in the middle of Figure 1. The periphery summarizes the net effects of peripheral organs like the liver or the kidney. This means that we do not explicitly consider the liver and the kidney in the model. We rather summarize their net effects in the periphery. An explicit representation would certainly be helpful in several ways, e.g., to investigate whether the low T 4 concentrations of AHDS patients can be explained by an accumulation of T 4 in the kidney. However, an explicit representation would go along with new parameters for which numerical values are needed and new states that must be considered. This renders the model even more complex and induces more uncertainty, because the parameters would not necessarily be uniquely identifiable.
In these peripheral organs, D1 and D2 convert F T 4 into F T 3 , a process which is again modeled by means of a Michaelis-Menten kinetics with G D1 , G D2 , K M 1 and K M 2 . This peripheral production is finalized by considering a first-order lag element, taking into account dilution (α 31 ) and clearance (β 31 ). The binding of T 3 to T BG and the respective dissociation constant (K 41 ) are also visible.
The production of T 3 in the pituitary is considered explicitly in the top right-hand corner of Figure 1. It depends on the maximal activity of D2 (G D2 ), the respective dissociation constant (K M 2 ), the dilution (α 32 ), the clearance factor (β 32 ) and the dead time (τ 3Z ). The central T 3 (T 3z ) mainly serves as a feedback signal to the pituitary.
Inside the pituitary, which is illustrated in the top of the scheme, many different processes take place. One part of the T 3z binds to the intracellular-binding-substrate (IBS), which does not serve as a feedback signal. The damping constant D R models the fact that T 3 must bind to specific thyroid hormone receptors in order to influence the concentration of T SH. The remaining constants (G R and L S ) stand for the maximum gain of the pituitary receptors in relation to thyroid hormones and a breaking constant, respectively.
Furthermore, the influence of T RH is illustrated in the top of Figure 1. The receptors for T RH at the pituitary are taken into account with a damping constant (D H ). The ultra-short feedback loop of T SH on its own secretion (compare (2)) is considered with the respective dilution (α S2 ) and clearance (β S2 ) factors. The maximum secretory capacity of the pituitary (G H ) is shown. The entrance of T SH into the bloodstream is modeled with a first order lag element, with dilution (α S ) and clearance (β S ).
After having derived the relationships between the different hormone concentrations, one must find the numerical values of all introduced parameters. The parameter values used within this work are listed in Section S9 below. Ideally, we would like to use only human parameters. However, this is not possible for all parameters, consider, e.g., the damping constant of T RH at the pituitary D H , which is not measurable in humans. Therefore, we partially exploit experimentally determined murine numerical parameter values. This results in a model which contains human and murine numerical parameter values. Nevertheless, this approach is meaningful for two reasons: first, there will always be a variation in the exact numerical parameters even for humans only. Second, when applying the mathematical model, we pursue the objective to get insight about the mechanisms of the AHDS, which is also possible when using murine and human numerical parameter values jointly. This is the case since (slightly) different parameter values lead to the same qualitative behavior of hormone concentrations in our model (compare also the sensitivity analysis in (3)).

S2 STATE TRANSFORMATION
Before stating the formal problem of the constrained parameter optimization for healthy individuals and AHDS patients in Sections S3 and S4, respectively, it is useful to introduce the system's differential equations for the Michaelis-Menten modeling of the membrane transporters T SH(t) T SH(t)+k Dio with the relationships T RH(t) = T RH 0 (1 + 0.6cos(2π(t/86400))) (S10) For the linear modeling of the membrane transporters, one must replace the term in equations (S1) and (S2) by k l T 4,th . (S12) As mentioned in the main part of this paper, we perform a constrained parameter optimization in order to identify the parameters G T , G D1 , G T 3 and G M T (or k l ). In other words, we are looking for the configuration of parameters which fits best the given real measured hormone data under the condition that the (steady-state or dynamic) differential equations (S1) -(S6) are satisfied.
One challenge of this approach is that the order of magnitude of the differential equations is very different. The value of F T 3 is of order 10 −12 , whereas the order of T SH is 1. In order to obtain a numerically well conditioned problem, we hence perform the following state transformation with (S14)

Frontiers
This transformation matrix leads to (transformed) hormone concentrations in the identification (F T 3 , F T 4 , T 4,th and T SH) that are in the order of magnitude of 10 0 . The dynamics of the transformed system arė where f (x) denotes the right hand side of the differential equations (S1) -(S6) and x, z as defined in (S13). When the mentioned values are plugged in, the expressioṅ Frontiers is obtained. The transformed differential equations are T SH(t)+k Dio T SH(t)+k Dio T SH(t)+k Dio T SH(t)+k Dio (S22)

Frontiers
Furthermore, the relationships of the hormone concentrations need to be considered. Using this state transformation, all variables are in the same order of magnitude.

S3 PARAMETER ESTIMATION FOR HEALTHY INDIVIDUALS USING DYNAMIC HORMONE MEASUREMENTS
As mentioned in the main part, we here explain the parameter estimation for healthy individuals using measured mean dynamic hormone concentrations F T 3,meas (t), F T 4,meas (t), and T SH meas (t). To this end, we use the dynamic hormone concentrations of healthy individuals documented in (4). The concept is to minimize the normalized quadratic error between measured (dynamic) hormone concentrations and simulated hormone concentrations. In other words, we want to find the configuration of the G T , G D1 , G T 3 , and G M T parameters that explains best the given measured (dynamic) hormone concentrations. In mathematical terms, we minimize the objective function , and T SH model (t) are the (back-transformed) numerical solutions to the transformed system of differential equations (S17) -(S22) 1 with the constraints that G T ≥ 0, G D1 ≥ 0, G T 3 ≥ 0, and G M T ≥ 0. The constants F T 3,meas , F T 4,meas and T SH meas are the mean hormone concentrations of all patients and all time points. The cost function sums up the normalized quadratic difference lasting 15 days. To generate dynamic hormone profiles lasting 15 days, we duplicate the 24 hours profile, see (4), 15 times. This longer hormone profile is necessary to guarantee that the transient and the steady-state behavior is considered in the estimation. This approach implies the assumption that the measured dynamic hormone concentrations remain unchanged throughout 15 days.
After solving this optimization problem, we observed the following: different initial guesses of the optimization problem lead to different minima with the same cost. Consequently, there are different parameter configurations of G T , G D1 , G T 3 and G M T that all explain equally well the given dynamic hormone measurements. This problem was also observed, when only steady-state hormone measurements are used for the identification of parameters of the model, compare (3). To solve this problem, we incorporate the information into the cost function that 20 % of T 3 is produced inside the thyroid gland and the remaining 80 % in peripheral organs, according to (5, 6) 2 . This is realized by penalizing the difference to this relation in the cost function. On the right hand side of (S29), we add 200|0.8 − PR peri (t)| 2 , where PR peri (t) is defined as the part of T 3 which is produced in the periphery, i.e., Hence, the updated cost function is Since the remaining production of T 3 will necessarily take place in thyroid cells (compare eq. (S3)), we do not need to incorporate a term related to the production of T 3 in thyroid cells. Using this additional information, we obtain unique optimal parameter values.
After having determined the unique optimal parameter configuration, we simulate the hormone concentrations 100 times and artificially corrupt the simulated concentrations by some noise following a normal distribution with µ = 0 and σ = 0.1. For each dataset, we estimate the parameters applying the aforementioned procedure. The results are shown in Table 1 of the main part.
Further, to quantify the uncertainty of the parameter estimation by means of individual hormone measurements, we determine the parameters G T , G D1 , G T 3 and G M T individually for 27 different healthy individuals, documented in (4). This individual estimation leads to the results presented in Table S1 and Section S5.
The estimation related to the linear modeling of the membrane transporters works completely analogous. One must simply replace the G M T T 4,th K M T +T 4,th term by k l T 4,th and estimate k l instead of G M T . The results are presented in Table 2 of the main part and Table S3 (regarding the individual parameter estimation).

S4 PARAMETER ESTIMATION FOR AHDS PATIENTS USING STATIC HORMONE CONCENTRATIONS
As mentioned in the main part of this paper, we perform a constrained parameter optimization using mean steady-state hormone concentrations in order to estimate the parameter G M T (or k l ) for AHDS patients. In other words, we are looking for the configuration of parameters which fits best the given real measured hormone data under the condition that the steady state equations are satisfied. In formal words, the task is to minimize the objective function subject to the nonlinear equality constraints The bounds are meaningful given that the parameter G M T as well as the hormone concentrations physiologically only make sense when they are positive. We only deal with steady-state expressions, therefore, the time dependence of the variables is neglected. Note that we plugged in the steady state expressions of the differential equations of T 3c and of T SH z into the steady state expression of T SH. This leads to condition (S38). The detailed steps how to reach (S38), which are cumbersome but straightforward, are shown in Section S8 for completeness. The results are presented in Table 1 of the main part.
We solve this optimization problem for the mean steady-state hormone concentrations of 13 different AHDS patients given in (7,8,9,10,11). This enables us to perform parametric bootstrapping. After having determined the optimal value of G M T , we simulate the hormone concentrations and corrupt the F T 3 , F T 4 and T SH concentrations by some noise following a normal distribution with µ = 0 and σ = 0.1. In this way, we generate dynamic hormone concentrations of AHDS patients. Therefore, the optimal parameters of the generated hormone concentrations are determined by the procedure outlined in Section S3 with the only difference that we need to estimate only G M T and not additionally G D1 , G T , and G T 3 .

S5 UNCERTAINTY QUANTIFICATION BASED ON INDIVIDUALLY ESTIMATED PARAMETERS
In this section, we show the results regarding the quantification of the uncertainty of the estimated parameters when these are estimated individually using the approach explained in Sections S3 and S4. First, we show the results for the Michaelis-Menten modeling of the membrane transporters and, second, we show the analogous results for the linear modeling of the membrane transporters.

S5.1 Michaelis-Menten Modeling
The results of the uncertainty quantification of the G D1 , G T 3 , G T and G M T are shown in Table S1. In general, the uncertainty quantification based on individually estimated parameters yields similar results as the uncertainty quantification based on bootstrapping in the main part. Except for the G T 3 parameter, the numerical values of the parameters are not subject to large uncertainties. Here, we can additionally show that the difference in the mean of the G M T parameter is significant (based on a two-sample t-test with a significance level of 5 %).
Next, the computed steady-state hormone concentrations are illustrated in Table S2. One can observe that the hormone concentrations of healthy individuals and of AHDS patients are not subject to large variations. All hormone concentrations have a coefficient of variation which is below 40 %. On one side, we observe the characteristic hormone concentrations of AHDS patients. The mean F T 3 concentration is 49 % times higher, the mean F T 4 concentration is 41 % lower, and the mean T SH concentration is 38 % higher for AHDS patients compared to healthy individuals. Interestingly, the mean T 4,th concentration is approximately 800 % higher for AHDS patients compared to healthy individuals. Moreover, the mean hormone concentrations differ significantly from healthy individuals to AHDS patients (again based on a two-sample t-test with a significance level of 5 %).
Finally, as in the main part, we simulate the hormone concentrations of healthy individuals and AHDS patients, as shown in Figure 2, using the mean parameters of Table S1. Since the mean parameters are similar to the one determined using bootstrapping, the hormone concentrations in Figure 2 are similar to the ones shown in Figure 3 of the main part.

S5.2 Linear Modeling
In this section, we consider the linear modeling of the membrane transporters. As visible in Table S3, the uncertainty quantification of the parameters yields once again similar results compared to the application of bootstrapping, compare Table 2 Table S1. concentrations do not differ considerably from the Michaelis-Menten modeling approach. In Table S4, we observe once again the characteristic hormone concentrations of AHDS patients, a (46 %) higher F T 3 and a (36 %) higher T SH concentration together with a (40 %) lower F T 4 concentration compared to healthy individuals. The T 4,th concentration is approximately 800 % for AHDS patients compared to healthy individuals. Finally, in Figure 3, the simulation results of the dynamic hormone concentrations are shown.
In conclusion, we want to point out that we applied two different approaches to quantify the uncertainty (bootstrapping and individual parameter estimation) of the parameters and that both yielded similar results in terms of how uncertain which parameters are.   (1, 3). The remaining parameters are estimated through a constrained parameter optimization approach and shown in Table S3.

S6 MCT8-MEDIATED T 3 TRANSPORT
As mentioned in the main part, we here consider additionally membrane transporters for T 3 . To this end, we proceed similarily as for the membrane transporters for T 4 . We introduce a new state called T 3,th which represents the T 3 concentration in thyroid cells. Its associated differential equation is T SH(t)+k Dio T SH(t)+k Dio The production of T 3 in thyroid cells depends on (i) the conversion of T 4 into T 3 in thyroid cells and (ii) and a direct T 3 synthesis based on the T SH concentration. The dilution factor α th is the same as in the differential equation of T 4,th (since we assume that the intrathyroidal T 3 has the same volume of distribution as the intrathyroidal T 4 ). The choice of β th,T 3 is more complicated. From (12), we can deduce that the intrathyroidal half-life of T 4 is approximately 44 hours, which is much shorter than the plasma half-life of T 4 corresponding approximately to 7 days. Assuming that the intrathyroidal half-life of T 3 is reduced in the same way (compared to the plasma half life of T 3 ), its half-life would correspond to approximately 6 hours and the corresponding clearance rate to β th,T 3 = 3.06 · 10 −5 s −1 . Once again, slightly different numerical parameter values will lead to the same qualitative behavior 3 .
Next, the differential equation of the peripheral T 3 concentration (called T 3p ) must be adapted. We now consider positively the part of T 3,th which is transported out of the thyroid cells and the peripheral conversion of T 4 into T 3 , yielding containing the new parameters G M T,T 3 and K M T,T 3 . Once again, the value of K M T,T 3 has already been determined experimentally in (13) and its value is applied here (although this is a simplification since different mutations will lead to different sensitivities). The block diagram in Figure 4 illustrates the incorporation of membrane transporters for T 3 (compare the block with a green frame) and T 4 (compare the block with a blue frame).
The parameter G M T,T 3 describing the maximum activity of the T 3 transport must be estimated for healthy individuals. We estimate this parameter by following the same concept, as described in Section S3 (using cost function (S33)) with the only difference that we additionally identify the parameter G M T,T 3 . Subsequent simulations for one healthy individual subject are shown in Figure 5. In this figure, one can see that the modeled hormone concentrations still fit to the measured hormone concentrations. Additionally, one can still observe the circadian rhythm of F T 3 . Hence, our model remains a valid approximation of the pituitary-thyroid feedback when considering M CT 8-mediated T 3 transport.
Even though the hormone concentrations for healthy individuals can be accurately represented considering membrane transporters for T 3 , the parameters become structurally non-identifiable. In other words, substantially different (of several orders of magnitude) parameter configurations lead to the same hormone concentrations. This difference to the estimation results when only membrane transporters for T 4 are considered is due to the additional parameter that we need to estimate, i.e., G M T,T 3 , while using the same real hormone concentrations measurements. In other words, we want to estimate an additional parameter without further data or information. Therefore, we do not consider membrane transporters for T 3 in the main part, although this approach is a simplification since it has been shown that the MCT8 transports T 3 in humans, compare (13). Nevertheless, this is a meaningful approach, since the T 3 export is most likely not harmed in M CT 8-deficiency (14). This intriguing observation could be explained by the existence of further membrane transporters as the MCT10 (which transports T 3 and occurs in the thyroid gland (15)) or the recently discovered SLC17A4 (which also transports T 3 (16)).

S7 STABILITY ANALYSIS
As mentioned in the main part of this paper, we performed a local stability analysis. To this end, we introduce Lyapunov's indirect method. Then, we document and discuss the results from a physiological point of view.

Lyapunov's Indirect Method
Lyapunov's indirect method can be formulated as follows (17): let x = 0 be an equilibrium point/hormone concentration for the nonlinear systemẋ where f : D → R n is continuously differentiable and D is a neighborhood of the origin. Let then 1. the origin is locally exponentially stable if Re λ i < 0 for all eigenvalues of A.
2. The origin is unstable if Re λ i > 0 for one or more of the eigenvalues of A.
As Lyapunov's indirect method is only a criterion for local exponential stability, no statement is possible regarding the region of attraction of the equilibrium point. The region of attraction of an equilibrium point is the set of all initial states (here: hormone concentrations), for which the solution of the nonlinear system (S47) asymptotically converges to the equilibrium point.

Results Lyapunov's Indirect Method
The stability analysis concerns the equilibrium hormone concentrations documented in Tables 1 and 2 of the main part, which were obtained using a constant concentration of T RH. The application of Lyapunov's indirect method leads to the results illustrated in Table S5. One can see that the eigenvalues are negative for all cases. Consequently, the computed equilibrium hormone concentrations are all locally exponentially stable.

Discussion Stability Analysis
After having documented these (rather technical) results, we can focus on their medical interpretation.
Local asymptotic stability means that if the hormone concentrations have values different from their equilibrium hormone concentrations, which are still in the region of attraction, the hormone concentrations will converge to their equilibrium hormone concentrations. Unfortunately, with the proposed method it is impossible to state how large the region of attraction is. Hence, it is impossible to quantify the maximal deviation of the hormone concentrations for which they converge to their equilibrium hormone concentrations.
This result holds true for healthy individuals and for AHDS patients. This is particularly interesting since AHDS patients have a highly perturbed pituitary-thyroid feedback loop. From a systems theoretic point of view, their equilibrium points are still locally asymptotically stable.
In the future, it would be highly valuable to determine the mentioned region of attraction. A possible method that goes along with a quantification of the region of attraction is Lyapunov's direct method (17). However, its application is highly challenging if the differential equations are complex as it is the case here.
The aforementioned results hold for a constant T RH concentration. Since the real T RH concentration in humans follows a pulsatile course, one might ask how this time-varying T RH concentration impacts the stability results. To this end, we note the following. As discussed above, the application of Lyapunov's indirect method shows that the equilibrium hormone concentrations (with constant T RH input) are locally exponentially stable. This implies that the (nonlinear) system (S1) -(S9) comprising the pituitary-thyroid feedback loop is locally input-to-state stable (ISS) (which follows, e.g., by applying Lemma 4.6 of (17)). This property (ISS), which has become one of the standard methods to quantify a certain notion of robustness of nonlinear systems, implies that small/bounded (but potentially time-varying) inputs (or deviations from a constant equilibrium input) result in small/bounded states (or deviations from the equilibrium state). In other words, small (potentially time-varying) deviations of T RH around a constant equilibrium value (e.g., a pulsatile course) result in small deviations of the hormone concentrations.