Mathematical Modeling of Thyroid Size and Hypothyroidism in Hashimoto’s Thyroiditis

Abstract

This chapter is devoted to studying the physiology of the pituitary-thyroid axis and thyroid size in autoimmune thyroiditis via modeling. The pituitary-thyroid axis consists of a feed forward and backward loop in humans, which is responsible for maintaining the body’s metabolism. Under a disease situation, the dynamics of the axis becomes more complex and unique among patients. Hashimoto’s autoimmune thyroiditis disrupts the normal operation of the axis by slowly destroying the thyroid follicle cells through complex immune mechanisms. So, the size of thyroid and the axis operation are fully, partly, or totally not functional in this disease. Basically, the patient situation in the disease process is unique in describing the diffused goiter and/or a clinical symptom of hashitoxicosis, euthyroidism, subclinical hypothyroidism, or overt hypothyroidism. Using patient-specific modeling, we can predict the hidden dynamics of the natural history of autoimmune thyroiditis and test hypothesis on the operation of axis. In addition, we unfold case studies of three patients from the thyroid literature through the modeling viewpoint and describe their hidden dynamics.

1. Introduction

The normal operation of the pituitary-thyroid axis depends on the levels of thyroid stimulating hormone (TSH) and thyroid hormones, triiodothyronine (T3) and thyroxine (T4) [1]. Serum TSH is produced and released by the pituitary gland in response to the low levels of free thyroid hormones in the serum. Circulating TSH in turn stimulates the thyroid to produce and secrete T3 and T4 into the serum. When thyroid hormones reach their highest levels, it inhibits the production of TSH, which describes the normal day to day operation of the pituitary-thyroid axis. The axis is commonly referred as an important negative feedback loop in the endocrine system. The normal function of this loop is essential for the body’s metabolic rate—it affects how quickly the cells in our body use the energy stored within itself [2, 3]. For the purpose of this modeling work, we use only serum free T4 for the levels of thyroid hormones as frequently asked to measure in the thyroid clinics for all patients with Hashimoto’s autoimmune thyroiditis (see Figure 1 ).

Figure 1.

The pituitary-thyroid axis is shown in this picture. It is commonly referred as a negative-feedback control.

The normal operation of the axis can be tested and verified clinically via one measurement of TSH and free T4 from the blood serum. Laboratory tests are important in diagnosing conditions of the thyroid gland [2]. The result of the blood test is determined approximately based on the normal reference range of TSH and free T4 is (0.4–4) mU/L and (7–18) pg/mL, respectively as recommended by the American Thyroid Association [2]. As the physiology of the axis is governed by TSH and free T4, the clinical state of the axis has been defined based on these values [3]. Suppose TSH and free T4 levels falls within the normal reference range; the state of the axis is said to be clinically normal. Suppose TSH levels falls above 4 mU/L but free T4 levels falls within the reference range; the state of the axis is said to be clinically subclinical hypothyroidism. Suppose TSH levels are above 4 mU/L and free T4 levels below 7 pg/mL; the state of the axis is said to be clinically hypothyroidism (an underactive thyroid gland). Keeping free T4 levels within the normal reference range is very important [3]. Lower levels of free T4 can cause hypothyroidism that results in several health problems including obesity, joint pain, infertility, slow metabolism, puffy face, constipation, stiffness, dry skin, depression, fatigue and higher heart rate [4].

Hyperthyroidism is another clinical state of the axis, in which free T4 levels stay above 18 pg/mL while the levels of TSH stay below the reference range 0.4 mU/L [5]. This state of the axis occurs when the thyroid gland over produces and secretes the thyroid hormones either in response to chronic TSH stimulation or due to Graves’ disease. Transient hyperthyroidism is a phenomenon that refers to the leakage of stored thyroid hormones from the gland, typically called Hashitoxicosis or thyroid burst [6]. It happens due to the Hashimoto’s autoimmune thyroiditis. The bursting of thyroid gland can happen at any clinical stage of the patients.

Hashimoto’s autoimmune thyroiditis is one of the immune disorders hosted by the thyroid gland [6, 7]. Under this disease process, the thyroid follicular tissue undergoes the slow destruction by the immune system and thereby the thyroid struggles to produce enough hormones for the body’s requirement. Currently, the incidence rate of Hashimoto’s autoimmune thyroiditis is estimated to be 300–500 cases per 100,000 individuals per year. In this disease, the interaction of the immune system with the thyroid is highly complex involving the cellular and humoral mechanisms. The involvement of humoral mechanism can be verified through a simple blood test for anti-thyroid antibody biomarkers. More specifically, the elevated titers of anti-thyroid antibodies, thyroid peroxidase antibodies (TPOAb) and/or thyroglobulin antibodies (TGAb) may indicate the malfunction of the immune system [8]. For simplicity in patient-specific modeling, we choose TPOAb as a biomarker representing the presence of the Hashimoto’s autoimmune thyroiditis.

In 1912, Hakaru Hashimoto published the description of four cases of diffused goiter which were showing the signs of infiltration of immune cells and antibodies in the thyroid gland [9]. All these goiters appeared differently from the colloid goiters that occur due to insufficient iodine in-take from the diet. In his time, the idea of autoimmune disease had not been established in the medical literature. He isolated these cases and asked the future clinicians to explore the mechanism of this new form of goiter [10, 11]. In 1936, the description of goiter due to lymphocytic thyroiditis was rediscovered in the United States and was labeled as Hashimoto’s thyroiditis. It has been characterized as an organ-specific autoimmune disease and swelling might be the result of the chronic stimulation of the thyroid gland by the serum TSH [12, 13]. Now, scientists have described autoimmune thyroiditis with atrophy, a variance from the original description of Hakaru Hashimoto.

The pituitary gland has no clue about the undergoing immune process in the thyroid gland and does not know how to adapt to the changing environment and the abnormal function of the thyroid [14, 15]. It simply tells the thyroid to keep up with the release of thyroid hormones for its signal. By responding through bulging its size, the poor thyroid gland needs to accomplish the task of producing enough levels of hormones if it could. What initiates the immune process against the thyroid gland? The answer is still largely unknown to the researchers, but they all speculate the immune process might be initiated through a complex combination of the genetic and environmental pollution [15]. Herein, we care about the functional size of thyroid gland as opposed to the exact size of the thyroid gland, which is treated as a hidden compartment [16]. The exact size of the thyroid gland includes the size of the parathyroid glands, the isthmus, capillaries, blood flow and so on whereas the functional size includes the follicular cells that contains thyroglobulin molecules in which hormones are stored. Basically, the functional size of the gland counts the part of the gland that is able to make and secrete hormones; obviously this information is not available in the clinical setting.

A measurement of TPOAb titers is required from the blood test to diagnose autoimmune thyroiditis and afterward, TPOAb titers is not usually measured as its levels not helpful in the treatment. As mentioned above, the axis function is determined from one-time measurement of TSH and free T4, respectively. However, the functional size of the thyroid gland cannot be measured through the laboratory experiments or with any medical tools, which does not exist in practice at least so far. Even if it exists, then it would be very expensive in terms of cost and time. However, the mathematical model can help by telling us the functional size from the clinical measurements of specific patient and from the physiological point of view. Moreover, the functional size of the gland is different for everyone, such as boy’s thyroid size versus men’s thyroid size or girl’s thyroid size versus women’s thyroid size and so on. Using the model and the clinical measurements of TSH and free T4, we will determine the initial functional size of the gland for a given patient. Goiter due to this disease occurs more frequently in women than men [17]. Using the model, we describe the diffused goiter as the functional size increases to keep up with the normal production of hormones.

In reality, it is hard to perform experiments on patients in the laboratory setting and collect data to understand the underlying dynamics of the pituitary-thyroid axis in autoimmune problem. Hypothetically speaking, if the experiment is possible, then one needs to consider patients’ time, cost, safety, and variability. Two patients with similar characteristics from societies might react to an experiment in an unusual way [18, 19]. The inter and intra variability among patients causes a major problem for treatment and challenges the obtained information from labs. Using the model, we can carry out an experiment for a specific patient and the hidden dynamics of the axis can be investigated at the microscopic level for shorter or longer time-period. In general, the model is implemented through a computer program and the parameters are tuned to a specific patient and then the experiments can be performed with the test hypothesis [20]. Also, we can identify sensitive and insensitive parameters that are responsible for the physiology of the axis under autoimmune thyroiditis, that causes the hypothyroidism and goiter. This is an effective and modern way to explore the complex interaction of the immune system to the thyroid and its consequence on the negative feedback loop [21, 22].

The hormone TSH has a half-life of 1 hour in serum on the average scale. Similarly, the hormone free T4 has a half-life of 7 days in serum on average, the thyroid per-oxidase antibodies (TPOAb) has a half-life of 24 hours in serum on average and the half-life of functional size of the thyroid gland is not known and probably varies among individuals. Several time scales have been involved in the disease process. Based on the fundamental principle of rate laws, the model will be developed here to predict the history of patient-specific dynamics of the axis and thyroid size in autoimmune thyroiditis [23, 24, 25]. The model can unfold the consequences of the presence of antibodies titers in the blood. Basically, it can replace humans but mimics the system that provides convenience and flexibility for scientists to run experiments on the computer. Information obtained through the dynamics might be useful in treating patients and improving the accuracy of the levothyroxine treatment [26]. For instance, the levothyroxine drug can be targeted and administered to the specific patient so that the levels of TSH and free T4 maintained within the normal reference range, which in turn can avoid other health problems.

The remaining of this chapter is divided into three main sections. Section 2 introduces the development and construction of a coupled model describing the interactions of two subsystems (the axis and humoral immune system) and reduces the coupled model based on the physiological assumptions. Section 3 analyzes the solutions of the reduced model qualitatively and numerically for various parameter values. Section 4 describes the case studies of three patients, validity of the model and provides the conclusion of this work.

2. Model outline

The model is constructed as a system of four ordinary differential equations based on the following known physiological assumptions [20, 21].

2.1 Assumptions

1. The pituitary gland is diseased free, so the feed forward loop is intact.

2. Total TSH receptors concentration does not change during Hashimoto’s autoimmune thyroiditis.

3. Serum TSH stimulates the growth of functional thyroid and the production and secretion of thyroid hormones.

4. The humoral immune system uses serum TPOAb to attack the thyroid and those titers can be used as a biomarker for the level of the anti-thyroid immune activity.

5. The patient does not demonstrate central or peripheral resistance to thyroid hormone.

6. Serum free T4 have much slower dynamics in the disease process compared to serum TSH, serum TPOAb and the functional size of the gland.

We denote $x\left(t\right)$ , $y\left(t\right)$ and $w\left(t\right)$ to represent the concentrations of serum TSH (mU/L), serum free T4 (pg/mL) and serum TPOAb (U/mL), respectively. We let $z\left(t\right)$ represent the functional size of the gland. According to the fundamental principle of rate laws, the rate of change of concentration of TSH over time $t$ is the secretion rate from the pituitary gland minus the elimination rate of TSH via kidney through the unspecified mechanism. The rate of change of concentration of free T4 over time $t$ is the secretion rate of free T4 from the thyroid minus the elimination rate of free T4 via kidney through unspecified mechanism. Similarly, the rate of change of concentration of TPOAb over time $t$ is the production rate due to the abnormal interaction thyroid gland and immune system minus the elimination rate of TPOAb via kidney through unspecified mechanism. Finally, the rate of change of the functional size of thyroid gland over time $t$ is the growth rate minus the destruction rate of the gland.

Next, we write a coupled model that represents the interaction of two subsystems: the negative feedback loop and the humoral immune system.

$$\frac{\mathit{dx}}{\mathit{dt}}=k_1-\frac{k_{1}y\left(t\right)}{k_a+y\left(t\right)}-k_{2}x\left(t\right),x\left(0\right)=x_0$$ E1

$$\frac{\mathit{dy}}{\mathit{dt}}=\frac{k_{3}z\left(t\right)x\left(t\right)}{k_d+x\left(t\right)}-k_{4}y\left(t\right),y\left(0\right)=y_0$$ E2

$$\frac{\mathit{dz}}{\mathit{dt}}=k_5\left(\frac{x\left(t\right)}{z\left(t\right)}-N\right)-k_{6}z\left(t\right)w\left(t\right)z\left(0\right)=z_0$$ E3

$$\frac{\mathit{dw}}{\mathit{dt}}=k_{7}z\left(t\right)w\left(t\right)-k_{8}w\left(t\right)w\left(0\right)=w_0$$ E4

where $x\left(t\right)\ge 0,y\left(t\right)\ge 0,z\left(t\right)>0$ and $w\left(t\right)\ge 0$ . The normal average value for each state variable is given in Table 1 . The model has 11 positive parameters unchanged over time, in which each parameter has a unique physiological meaning and their corresponding values listed in Table 2 . The model (1)–(4) captured the dynamics of two subsystems coupled through the functional size of the thyroid gland. Using assumption (8) that serum free T4 has much slower dynamics compared to serum TSH, serum TPOAb and the functional size of the thyroid gland, we can take $\mathit{dy}/\mathit{dt}=0$ and obtain a reduced model consisting of three differential equations and one algebraic equation:

Name	Normal value	Normal range	Source	Unit
$x\left(t\right)$	1	$0.4-\left(2.5-4\right)$	Literature [2]	mU/L
$y\left(t\right)$	13	$7-18$	Literature [2]	pg/mL
$z\left(t\right)$	0.015	$0.005-0.125$	Literature [16]	L
$w\left(t\right)$	0	$0-<200$	Dataset	U/mL

Table 1.

Variable names, normal values, ranges, sources and units.

Name	Normal value	Normal range	Source	Unit
$k_1$	$5000$	$>4000$	Literature [27]	mU/L day
$k_2$	$16.6$	N/A	Literature [3]	1/day
$k_3$	$90.11$	$4-257$	Simulation	Pg/mL L day
$k_4$	$0.099021$	N/A	Literature [3]	1/day
$k_5$	$1$	N/A	Simulation	$L^3$ /mU day
$k_6$	$1$	N/A	Simulation	mL/U day
$k_8$	0.035	N/A	Literature [28]	1/day
$k_a$	$0.043$	$0.02-0.06$	Calculation	pg/mL
$k_d$	$0.05$	N/A	Simulation	mU/ $L^2$

Table 2.

Parameter names, normal values, ranges, sources and units.

$$\frac{\mathit{dx}}{\mathit{dt}}=\frac{k_1{k}_4{k}_a\left(k_d+x\left(t\right)\right)}{k_4{k}_a{k}_d+k_4{k}_{a}x\left(t\right)+k_{3}z\left(t\right)x\left(t\right)}-k_{2}x\left(t\right),x\left(0\right)=x_0$$ E5

$$z=\frac{k_{4}y\left(k_d+x\right)}{k_{3}x}=f\left(x,y\right)$$ E6

Using (Eq. (6)), the functional size can be determined for patients if a result of TSH ( $x$ ) and free T4 ( $y$ ) value is known from the blood test, which in turn can be used as an initial condition on (Eq. (3)) for the model simulation. The graph of (Eq. (6)) is shown on the Figure 2 , which illustrates the functional size of thyroid gland in terms of varying TSH and free T4 values. Using the reduced model, we can describe the clinical progression of overt hypothyroidism and thyroid size such as goiter and atrophy.

Figure 2.

The functional size of the thyroid gland can be calculated from the graph of this function of TSH and free T4.

3. Model analysis

As a first step in the process of analyzing the model, we solve for steady states of the reduced model (Eqs. (3)–(5)) from the following equations:

$$\frac{k_1{k}_4{k}_a\left(k_d+x\right)}{k_4{k}_a{k}_d+k_4{k}_{a}x+k_3\mathit{zx}}-k_{2}x=0$$

$$k_5\left(\frac{x}{z}-N\right)-k_6\mathit{zw}=0$$

$$k_7\mathit{zw}-k_{8}w=0$$

which leads to diseased-free and diseased steady state solutions besides the initial condition. The first steady state (diseased-free) denoted as $\left(z_1,w_1,x_1\right)$ where $z_1=x_1/N,w_1=0$ and the value of $x_1$ is given by the cubic equation:

$$a_1{x_1}^3+a_2{x_1}^2+a_3{x}_1+a_4=0$$ E7

where

$$a_1=\frac{k_2{k}_3}{N}>0$$

$$a_2=k_2{k}_4{k}_a>0$$

$$a_3=k_4{k}_a\left(k_d{k}_2-k_1\right)<0\text{since}{k}_1>k_d{k}_2$$

$$a_4=-k_4{k}_a{k}_1{k}_d<0$$

By Descarte’s rule of signs, Eq. (7) has one positive real solution, so the reduced model has diseased-free state in the positive octant for the system parameters. In fact, this steady state solution lives on the surface of the function $z=f\left(x,y\right)$ given by Eq. (6). The operation of the pituitary-thyroid axis is healthy and fully functional in the presence of diseased-free steady state solution.

Next, the second steady state is the diseased state solution denoted as $\left(z_2,w_2,x_2\right)$ where

$$z_2=\frac{k_8}{k_7}$$

$$w_2=\frac{k_7{k}_5}{k_6{k}_8}\left(\frac{k_7{x}_2}{k_8}-N\right)$$

and the value $x_2$ is given by the quadratic equation:

$$b_1{x_2}^2+b_2{x}_2+b_3=0$$ E8

where

$$b_1=\left(k_2{k}_4{k}_a+\frac{k_2{k}_3{k}_8}{k_7}\right)>0$$

$$b_2=k_4{k}_a\left(k_d{k}_2-k_1\right)<0\text{since}{k}_1>k_d{k}_2$$

$$b_3=-k_4{k}_a{k}_1{k}_d<0$$

By Descarte’s rule of signs, Eq. (8) has one positive real solution when $x_2>\frac{N{k}_8}{k_7}$ so the reduced model has diseased state in the positive octant for the system parameters.

Remark: When $N=k_7{x}_2/k_8$ , the second steady state (diseased state) of the system (Eqs. (3)–(5)) has $w_2=0$ . So, it must coincide with the first steady state (diseased-free state) in the positive quadrant if

$$N=\frac{k_7{x}_2}{k_8}=\frac{k_7{x}_1}{k_8}$$

3.1 Definition: bifurcation parameter

We let $N^{*}$ be the unique value of N, where

$$N^{*}=\frac{k_7{x}_0}{k_8}=\frac{k_7{x}_2}{k_8}=\frac{k_7{x}_1}{k_8}$$

the system undergoes a bifurcation. We call $N^{*}$ is a bifurcation value of the system (Eqs. (3)–(5)) and $N$ is the bifurcation parameter. We can take the initial state of the system at where the diseased-free and diseased steady states merge together. Using this definition, we can also calculate the parameter $k_7$ (say ${k_7}^{*}$ ) provided the values of $k_8$ and $x_0$ is known.

3.2 Equation of tangent plane

We first define the level surface (S) through a function of three variables from Eq. (6) as

$$g\left(x\left(t\right),y\left(t\right),z\left(t\right)\right)=\frac{k_{4}y\left(t\right)\left(k_d+x\left(t\right)\right)}{k_{3}x\left(t\right)}-z\left(t\right)=0$$ E9

The normal surface vector $\overrightarrow{n}$ can be calculated from the gradient of this function

$$\overrightarrow{n}=⟨\frac{\partial g}{\partial x},\frac{\partial g}{\partial y},\frac{\partial g}{\partial z}⟩$$

In particular, the normal vector of the surface $S$ at initial state $P\left(x_0,y_0,z_0\right)$ is a vector perpendicular to the tangent plane of $S$ at $P$ is given by

$$\overrightarrow{n}=⟨\frac{-k_4{k}_d{y}_0}{k_3{x_0}^2},\left(\frac{k_4{k}_d}{k_3{x}_0}+\frac{k_4}{k_3}\right),-1⟩$$

The equation of the tangent plane at point $P\left(x_0,y_0,z_0\right)$ is given by

$$\mathit{ax}+\mathit{by}+\mathit{cz}=d$$

where

$$a=\frac{-k_4{k}_d{y}_0}{k_3{x_0}^2},b=\left(\frac{k_4{k}_d}{k_3{x}_0}+\frac{k_4}{k_3}\right),c=-1\wedge$$

$$d=a{x}_0+b{y}_0+c{z}_0$$

Next, we take the implicit differentiation of Eq. (9) with respect to time $t$ and substituting the initial point $P\left(x_0,y_0,z_0,w_0\right)$ , which results in the following equation for $N$ :

$${\left(\frac{\partial g}{\partial x}\bullet \frac{\mathit{dx}}{\mathit{dt}}+\frac{\partial g}{\partial y}\bullet \frac{\mathit{dy}}{\mathit{dt}}+\frac{\partial g}{\partial z}\bullet \frac{\mathit{dz}}{\mathit{dt}}\right)}_P=0$$

$$N=\frac{-k_4{y}_0\left(-k_1{k}_a{k}_d-x_0\left(-k_2{k}_d+k_4\left(k_d+x_0\right)\right)\left(k_a+y_0\right)\right)}{k_5{k}_3{x}_0^2\left(k_a+y_0\right)}+\frac{x_0}{z_0}-\frac{\left(k_4+k_6{w}_0\right)z_0}{k_5}$$ E10

Notice that $N$ is independent of the parameters $k_7$ and $k_8$ . So, we can use this value of $N$ as $N^{*}$ . This threshold value can be uniquely determined for every patient based on their initial blood test results and if we know other parameters $k_1,k_2,k_3,k_4$ , $k_5,k_6,k_a$ and $k_d$ . We have estimated the values of all other parameters from the literature or through calculation (see Table 2 ).

3.3 Linear stability

As the bifurcation parameter $N$ varies from the threshold value $N^{*}=k_7{x}_0/k_8$ , the number of the steady states changes, and its behavior near the steady states (linear stability) changes as well. Moreover, as the bifurcation parameter varies, the appearance and disappearance of steady states can be seen on the tangent plane at the initial state of the function $z=f\left(x,y\right)$ . For this system (Eqs. (3)–(5)), there are three cases, namely when $N<N^{*}$ , $N=N^{*}$ and $N>N^{*}$ . Figure 3 illustrates a nice overview of how the solutions become stable and unstable as parameter $N$ changes from 0 to 300.

Figure 3.

The bifurcation diagram shows the steady state values of TSH as the parameter N varies from 0 to 300. BP means the branch or bifurcation point where another equilibrium curve passes through and the system switches its stability. The parameter value of N at BP is 66.7. H is a neutral saddle, but not the bifurcation point for the equilibrium curve. At H, there is a homoclinic orbit for the system.

Case 1: $N<N^{*}$

When $N<N^{*}$ , the system (Eqs. (3)–(5)) has two steady states, which undergoes a saddle-node bifurcation—disappearance of the saddle diseased state. One can imagine if we start the system from the initial steady state, which splits into the diseased-free and diseased state. The disease state is asymptotically stable which means all its eigenvalues are purely negative real. However, the diseased-free state is saddle because two of its eigenvalues are positive real and one of its eigenvalues is a negative real. See [22] for proof of the stability of the diseased-free and diseased steady state.

As the value of $N$ decreases from the threshold value and when it is close to zero, there exists a homoclinic orbit associated at the initial state which forms the boundary for the asymptomatically stable interior diseased state. The orbit captures the hidden dynamics of the diffused goiter and hashitoxicosis (see Figures 4 and 5 ).

Figure 4.

When the value of N = 10 and ${\mathit{k}}_{\mathit{7}}=\mathit{2.3345}$ , the hidden dynamics of the axis and the thyroid size reveal hashitoxicosis and goiter. The model simulation is performed for 2 years from the initial state: TSH = 1 mU/L, free T4 = 13 pg/mL, and thyroid size = 0.015 mL and TPOAb = 10 U/mL.

Figure 5.

When N = 0.01 and ${\mathit{k}}_{\mathit{7}}=\mathit{2.3345}$ , the hidden dynamics of the thyroid size capture the development of goiter and return to normal size.

Case 2: $N=N^{*}$

When $N=N^{*}$ , the system (Eqs. (3)–(5)) has one steady state (diseased-free), which is asymptotically stable. Since there is only one steady state in the system, the operation of the pituitary-thyroid axis in the presence of anti-thyroid peroxidase antibodies can be tested. The numerical simulation of the model shows TPOAb slowly approaches zero as time increases over months (see Figure 6 ).

Figure 6.

The parameter value is set to $\mathit{N}=\mathit{66.7},{\mathit{k}}_{\mathit{7}}=\mathit{2.3345}$ and the initial state is at euthyroidism (1, 0.015, 10). The levels of anti-thyroid peroxidase decrease slowly whereas TSH and the functional thyroid size remains unchanged over time, which indicates the normal operation of the axis with fully functional thyroid size.

Case 3: $N>N^{*}$

When $N>N^{*}$ , the system (Eqs. (3)–(5)) has one steady state (diseased-free), which is asymptotically stable. As the value of $N$ increases beyond the threshold value, the diseased free state is moving on the function $z=f\left(x,y\right)$ describing the progression of thyroid size and the function of the axis. One can imagine if we start the system from the initial state of normal thyroid size, then the system may approach goiter or atrophy depending upon the parameter value of $N$ while keeping other parameters fixed. On the other hand, varying the parameter $k_7$ , the clinical mechanisms euthyroidism $\to$ subclinical hypothyroidism $\to$ overt hypothyroidism can be explained while keeping $N$ fixed [20]. See Figure 7 that shows the clinical progression of subclinical hypothyroidism and thyroid size.

Figure 7.

Clinical progression from euthyroidism to mild subclinical hypothyroidism when $\mathit{N}=\mathit{80},{\mathit{k}}_{\mathit{7}}=\mathit{2.3345}$ with initial state (1, 0.015, 10).

3.4 Exploration of parameter curve

Suppose that an individual is diagnosed with Hashimoto’s autoimmune thyroiditis; we can use the reduced model to explain the physical and clinical symptoms occurring for this individual in the course of this disease. Having this disease means everyone has a unique behavior and knowing that may be very helpful in managing the course. In fact, certain parameter values in the reduced model are responsible for the uniqueness in patient behavior. More specifically, we have identified through the stability analysis that the parameter $N$ can explain the size of the thyroid (goiter, normal or atrophy) whereas $k_7$ can explain the clinical progression of the disease [21, 22]. Putting these two parameters together, one can test the hypothesis on the patients’ symptoms. See Figure 8 that shows the parameter curve in terms of $N$ and $k_7$ that is responsible for different thyroid size and clinical conditions.

Figure 8.

The sample parameter curve of a Hashimoto’s patient is shown here. The hypothesis testing can be done for different ${\mathit{k}}_{\mathit{7}}$ and $\mathit{N}$ values from the curve.

Suppose a patient’s information is given; then a threshold $\left(k_7^{*},,,N^{*}\right)$ can be found for the patient using Eq. (10) and the definition of the bifurcation parameter. If another information is available in the future, then another threshold can be obtained. To be exact, the threshold is an ordered pair found on the curve given by the bifurcation definition of parameter (see Figure 8 ). To simulate physical and clinical symptoms, we take two test ordered pairs from the parameter curve below and above the threshold $\left(k_7^{*},,,N^{*}\right)=2.3345,66.7$ ) while other parameters came from Table 2 . With the first test point for $\left(k_7,,,N\right)=\left(\mathrm{0.9904,12.65}\right)$ , we simulated the model for a day and found the symptoms of the disease are hashitoxicosis and goiter (see Figure 9 ). Similarly, by keeping the second test point as (5.8, 397.9), the one-day simulation showed mild atrophy and overt hypothyroidism (see Figure 10 ). Finally, by keeping the third test point away from the curve as (4.5, 150), the one-day simulation showed subclinical hypothyroidism and normal size of the thyroid ( Figure 11 ).

Figure 9.

The result of physical and clinical symptoms is goiter and hashitoxicosis when test values $\mathit{N}=\mathit{12.65}$ and ${\mathit{k}}_{\mathit{7}}=\mathit{0.9904}$ are taken from the parameter curve.

Figure 10.

The result of physical and clinical symptoms is mild atrophy and overt hypothyroidism, when test values $\mathit{N}=\mathit{397.9}$ and ${\mathit{k}}_{\mathit{7}}=\mathit{5.8}$ are taken from the curve.

Figure 11.

The result of physical and clinical symptoms are normal thyroid size and subclinical hypothyroidism, when test values of $\mathit{N}=\mathit{4.5}$ and ${\mathit{k}}_{\mathit{7}}=\mathit{150}$ are not taken from the curve.

4. Case studies of three patients

Using patients’ information from the peer-reviewed published article, we will predict patients’ natural history of the disease. More precisely, the reduced model can describe thyroid size and clinical progression from euthyroidism to subclinical or overt hypothyroidism for each patient given below. These patients’ information was provided by Salvatore Benvenga originally and were already published in the article [21]. The data consists of TSH, free T4 and TPOAb information whose normal reference ranges are (0.4–2.5) mU/L, (7–18) pg./mL and (0–200) U/mL, respectively. Using Eq. (6) and data, we have computed the functional size of the thyroid gland. Using Eq. (10), data and the bifurcation definition, we have computed $N$ and $k_7$ for all patients given below (see Tables 3 – 5 ).

Time (months)	TSH (mU/L)	FT4 (pg/mL)	TPOAb (U/mL)	Thyroid size (mL)	$N^{\ast }$	${k_7}^{\ast }$
0	0.8282	13.5	50	15.73	51.865	2.1918
8	0.93	13	55	15.1	60.952	2.2939
30	1.178	11.8	159	13.5	84.998	2.5254

Table 3.

Patient 1 information.

Time (months)	TSH (mU/L)	FT4 (pg/mL)	TPOAb (U/mL)	Thyroid size (mL)	$N^{\ast }$	${k_7}^{\ast }$
0	1.54	11.78	3810	13.37	64.301	1.4614
31	1.57	12.36	1310	14.02	93.662	2.088
39	1.54	11.97	1480	13.58	93.294	2.1203

Table 4.

Patient 2 information.

Time (months)	TSH (mU/L)	FT4 (pg/mL)	TPOAb (U/mL)	Thyroid size (mL)	$N^{\ast }$	${k_7}^{\ast }$
0	1.46	15.16	164	17.23	81.908	1.9635
2	1.56	13.74	153	15.58	97.724	2.1925
3	1.85	13.09	191	14.77	122.4	2.3157
5	2.06	11.67	482	13.14	150.5	2.557
22	4.61	10.47	773	11.63	387.39	2.9411
39	5.04	7.63	537	8.467	590.65	4.1018

Table 5.

Patient 3 information.

Assuming all patients had diseased-free state (1, 13, 0.015, 0) at some point, the data provides us the diseased state from which one can back trace the physical and clinical conditions caused by the disease and experienced by the patients. We have taken three untreated patients with a known biomarker for the presence of the autoimmune thyroiditis. Patients 1 and 2 visited the clinic three times whereas patient 3 visited the clinic six times due to mild goiter or other symptoms. We do not have any information about their thyroid sizes in the data.

In Table 3 , the clinical measurements, the functional thyroid size and parameters are listed for patient 1 which form the diseased states. The model simulation started from each diseased state for 10 years to back trace the natural symptoms of patient 1 due to autoimmune thyroiditis. The course of disease revealed thyroid size changes from mild goiter to normal while the function of axis remained normal (see Figure 12 ). In Table 4 , the clinical measurements, the functional thyroid size and parameters are listed for patient 2. The model simulation is done for 1 year to back trace the symptoms from the course of the disease. The first disease course revealed a development of subclinical hypothyroidism with normal functional thyroid size (see Figure 13 ). Similarly, patient 3 diseased states are found in Table 5 . The model simulation is done only for 2 days to capture the early course of the disease. Surprising results are seen in the early course of the disease (see Figure 14 ). This patient had developed both clinical and physical symptoms at different points in time. In particular, the first course shows mild goiter and mild hashitoxicosis whereas the sixth course shows the overt hypothyroidism and mild atrophy.

Figure 12.

Three courses of the natural history of the disease is shown here for patient 1. Patient 1 visited the clinic three times and was not required any treatment based on their TSH and free T4 values. The model is simulated for 10 years in which the first course of disease shows thyroid size started with mild goiter and then returned to normal.

Figure 13.

Three courses of the disease is shown here for the second patient. This patient visited the clinic three times and was not required treatment based on their TSH and free T4 measurements. The model is simulated for 1 year to observe the hidden dynamics. The course of the disease showed patient had developed subclinical hypothyroidism with normal thyroid size.

Figure 14.

Six courses of the disease has been shown here for the third patient. This patient visited the clinic six times and was not required any treatment based on their TSH and free T4 measurement. The model is intentionally simulated for 2 days to capture the early dynamics of the patient. The model simulation shows this patient had developed both clinical and physical symptoms at different points in time.

5. Conclusion

The human body is made up of so many subsystems such as the pituitary-thyroid axis and the immune system. These subsystems do not disrupt the function of each other in healthy people. A function of the pituitary-thyroid axis is to secrete the appropriate levels of thyroid hormones and take care of body’s metabolism whereas the function of the immune system is to protect the body’s organs such as the thyroid gland and remove foreign substances that try to invade the organs. Hashimoto discovered the abnormal interaction of the immune system to the thyroid gland, which resulted in the disruption of the physiology of the axis. Later this diseased condition has been named as Hashimoto’s autoimmune thyroiditis. Hashimoto noticed the destruction of thyroid follicular cells through physical and clinical symptoms in four middle aged female patients.

A keystone of the functional thyroid gland is the follicular cells, which die due to the aggressive and destructive attack by the immune system. Hashimoto patients discovered the disease sometimes by themselves due to discomfort in the neck (small goiter) or accidentally during annual checkup by the family physicians. Consequences of Hashimoto disease can be classified into physical and clinical symptoms at various stages of the disease. Both types of symptoms occur in a sequential manner from one extreme to another. More precisely, the physical symptom runs from goiter $\mathbf{\to }$ atrophy whereas the clinical symptom runs from hashitoxicosis $\mathbf{\to }$ overt hypothyroidism. In order to describe these symptoms mathematically, we chose four state variables namely the size of functional thyroid gland to describe the physical symptom, TSH and free T4 to describe the clinical symptoms and thyroid peroxidase antibodies as a representative for the presence of Hashimoto disease. Modeling the disruption of the axis through a coupled model is the key to unlock the hidden dynamics experienced by the patients. The hidden dynamics can be seen in patients through physical and clinical symptoms.

The literature description of Hashimoto disease begins with a gradual swelling of the thyroid gland and development of mild clinical condition, euthyroidism or subclinical hypothyroidism and subsequent gradual progression of overt hypothyroidism. Small goiter and hashitoxicosis are the very early stages of the disease and typically go untreated and hidden in the view of patients and physicians. Overt hypothyroidism is the irreversible end clinical state (where levothyroxine treatment is needed) whereas atrophy is irreversible end physical state of the disease. Basically, the mechanism involved in the progression of the disease is unique and sequential. For instance, some patients may have the disease course of goiter and euthyroidism, some may have goiter and the clinical progression euthyroidism $\mathbf{\to }$ subclinical hypothyroidism and some patients may have gradual progression from normal thyroid size to atrophy and euthyroidism to overt hypothyroidism.

Herein, we have developed and used patient-specific model to describe all possible mechanisms involved in the autoimmune thyroiditis. This can be achieved using two parameters $\mathit{N}$ and ${\mathit{k}}_{\mathbf{7}}$ . To be precise, the parameter $\mathit{N}$ and ${\mathit{k}}_{\mathbf{7}}$ can trace the thyroid size and clinical progression. We validated the model with three patients’ data by describing their natural history of the disease.

Acknowledgments

I would like to dedicate this chapter to Dr. Stephen J. Merrill for his guidance on mathematical modeling of Hashimoto’s autoimmune thyroiditis.

Conflict of interest

We declare there is no conflict of interest on this chapter.

Notes/thanks/other declarations

I would like to thank Dr. Salvatore Benvenga for providing Hashimoto’s patients data for the modeling work.