A model of the cost of delaying treatment of Hashimotos thyroiditis: thyroid cancer initiation and growth

1.   Introduction

The functional unit in the thyroid gland is the follicle, which consists of follicular cells surrounding colloid ^[1]. The size of a follicle depends on the number of cells and the amount of colloid (see Figure 1). Circulating serum TSH binds to specific thyroid follicular cell surface receptors where it stimulates the production and secretion of thyroid hormones, triiodothyroxine (T3) and thyroxine (T4) into the serum ^[2]. Thyroid hormones help the body with the regulation of metabolic processes, normal growth and development. In fact, the hormones enter each cell nucleus throughout the body and regulate gene transcription. Low levels of T3 and T4 initiate the release of thyrotropin-releasing hormone (TRH) from the hypothalamus and thyroid stimulating hormone (TSH) from the pituitary. High levels of T3 and T4 inhibit the release of TRH and TSH (see Figure 2). Serum T3 and T4 are regulated by the hypothalamic-pituitary-thyroid (HPT) axis through a negative feedback control ^[1,2]. Each day the healthy thyroid gland produces $(15-30)$ $\mu$g of T3 and $(70-90)$ $\mu$g of T4 approximately ^[3,4]. Both hormones are active; however the nucleus of most cells has receptors that have ten fold greater affinity for T3. Since the body produces less T3, 80% of the T4 is converted to T3 for body's requirement ^[5].

Serum TSH regulates the growth of thyroid follicular cells by governing the internal cell division mechanisms. The normal reference interval of serum TSH levels is approximately (0.4–4) mU/L for adults with the normal function of the HPT axis ^[6]. Serum TSH level is often increased in the hypothyroidism due to the existence of the negative feedback loop in the HPT axis. Hashimoto's thyroiditis (HT) is the most common autoimmune disorder of the thyroid gland; it can cause hypothyroidism and goiter ^[7,8]. In HT, the thyroid follicular cells are damaged by the immune system through humoral and cellular mechanisms with a strong association of genes and environment ^[9]. The involvement of humoral mechanism can be verified through a simple blood test of anti-thyroid peroxidase antibodies (TPOAb) titers. More specifically, TPOAb titers can be considered as a biomarker for the tissue destructive processes associated with the hypothyroidism observed in Hashimoto's and atrophic thyroiditis.

As a result of the destruction of the cells, the functional follicle shrinks and the amount of colloid is reduced. Importantly, serum TSH levels increase gradually and rise above the upper reference limit, while thyroid hormone levels can be normal or low. Elevated serum TSH levels continuously stimulate the thyroid tissue or cells affected by the HT, which could be the main factor responsible for the development of thyroid nodules (lump or bump) and cancer ^[10]. Although most thyroid nodules are benign, some can undergo transformation and become cancerous. To be accurate, the direct relationship between serum TSH levels and thyroid cancer also exists within the reference interval of serum TSH itself. This risk intensifies further for TSH levels above the upper reference limit ^[11,12,13]. Furthermore, serum TSH levels can be used as an important biomarker for the diagnosis of HT, nodules and the presence of malignancy in a thyroid nodule ^[14,15].

Clinically, a nodule is defined as an abnormal accumulation of thyroid follicular cells at a site within the thyroid gland. The growth of cells can form a solid or fluid-filled (i.e., cystic) lump. According to current guidelines, a nodule size greater than 1 cm in diameter should be tested for one of the clinically stages - suspicious or positive for malignancy ^[16,17,18]. A small percentage of nodules become malignant, but this is considered rare since more than 90% of detected nodules are benign ^[19]. Benign thyroid nodules less than 1 cm require a follow-up with frequency of $38$ months as one study reported since benign nodule size increases slowly ^[20]. When thyroid-pituitary homeostasis is maintained, the steps leading to nodule formation are not expected to develop and the chances of nodule development are negligible. In HT, thyroid-pituitary homeostasis is disrupted and therefore TSH levels are chronically elevated, which leads to thyroid cell hypertrophy (increased size), hyperplasia (increased number of cells), and some potential for neoplasia (uncontrolled growth of cells) ^[21]. Thus, the HT patients are recommended for regular examination of nodule growth.

The ultrasound and fine-needle aspiration biopsy are tests used in thyroid clinics to identify and examine thyroid nodules ^[22,23]. The ultrasound is an imaging test in which sound waves are sent to take a picture of the gland and the enlarged lymph nodes. The ultrasound estimates the size and location of the nodule and whether it is solid, fluid-filled or mixed viz solid with a fluid component. However it cannot accurately differentiate between a benign and cancerous nodule ^[24]. Solid nodules that have little fluid or colloid are more likely to be cancerous than are fluid-filled nodules. Still, some types of solid nodules, such as hyperplastic nodules and adenomas, have too many cells, but the cells are not cancer cells. So, the fine-needle aspiration cytology (or biopsy) is the best investigation that can evaluate pre-operatively the malignancy risk of a nodule ^[23]. The biopsy removes a small sample of tissue for evaluation under the microscope. More precisely, Rizzo et al state that the fine needle aspiration biopsy were considered adequate when a smear contained at least six groups of well-preserved, well-visualized follicular cells, each group containing more than 10 cells ^[25]. Papillary microcarcinoma is frequently discovered in nodules less than 1 cm in diameter, often in tissue removed for treatment of benign disease ^[26]. The origin of these nodules appear to be related to the BRAF(V600E) mutation, and may not be initiated by TSH. Although, higher TSH may be associated with more agammaessive malignancy in these BRAF mutations ^[27].

Thyroid cancer occurs when cells in the gland undergo genetic mutations that allow the cells to grow, multiply and spread rapidly ^[21,28]. There are five types of thyroid cancer: papillary, follicular, medullary, anaplastic and thyroid lymphoma of which papillary and thyroid lymphoma are associated with the HT and higher serum TSH ^[11,29,30]. The types are based on how the cancer cells look under a microscope inspection. The papillary, follicular and anaplastic cancers arise from the mutations of follicular cells of the thyroid ^[31,32,33]. Anaplastic thyroid cancer is a rare agammaessive cancer that is very difficult to treat and has a median survival time not longer than $6-8$ months ^[33]. The medullary thyroid cancer begins in parafollicular cells called C cells, which produce the hormone calcitonin ^[34]. Thyroid lymphoma is a rare form of thyroid cancer that begins in the immune system cells in the thyroid and grows very rapidly ^[35]. Anaplastic and thyroid lymphoma occur frequently in the older population. Overall, thyroid cancer is most frequently diagnosed among people aged 45–54 ^[36] and the percent of thyroid cancer deaths is highest among people aged 75–84 ^[36]. In the United States from $2006-2015$, the number of new cases of thyroid cancer has been rising on average $3.1\%$ each year while the number of deaths has been rising on average $0.7\%$ per year ^[36].

2.   An outline of the HT model

Modeling the mechanisms of thyroid physiology coupled with disease process is a way to explore the underlying complex phenomenon of cancer ^[37,38]. In this section, we will provide a brief review of the construction and results of the Hashimoto's thyroiditis (HT) model. The modeling of HT is to explain the clinical progression of euthyroidism $\rightarrow$ subclinical hypothyroidism $\rightarrow$ hypothyroidism due to the destruction of thyroid follicles by the immune system ^[39]. The model considered four state variables namely $x$ to represent the concentration of TSH (in mU/L), $y$ to represent the concentration of FT4 (in pg/mL), $z$ to represent the functional size of thyroid gland (in L) and $w$ to represent the concentration of TPOAb (in U/mL) respectively. The HT model took the form of the system of differential equations with initial conditions:

where $x(t)\geq 0, y(t)\geq 0, z(t) > 0$ and $w(t)\geq 0$.

Initially, there was one steady state in the model for an appropriate parameter set which corresponded to the homeostasis or euthyroidism of an individual ^[38]. After the onset of Hashimoto disease or the malfunction of the immune system for a prolonged time period, the individual's steady state moved away from the euthyroidism - this behavior was captured through the bifurcation analysis of one parameter $k_7$ in the model while other parameters remained fixed. In the model, $k_7$ corresponded to the production rate of anti-thyroid peroxidase antibodies, a representation of the malfunction of the immune system. We analytically showed our main result of the model; when patients' $k_7$ value was less than $k_7^*$, serum TSH steady state value stayed within the normal reference limit whereas when $k_7$ value was larger than $k_7^*$, serum TSH steady state value increased above the upper normal reference limit while free T4 can be normal or low depending upon how high the $k_7$ value was set in the simulation. Notice that $k_7^*$ is a threshold of the system for the progression of Hashimoto disease. Figure (3) is a bifurcation diagram in which the steady state values of TSH are shown as $k_7$ varied in the model. These steady state values of TSH can be used as an input for the stochastic model of nodule formation and development.

3.   Concentration to molecules

For a given patient $k_7$ value, the above deterministic model can predict the blood concentration of TSH after $t$ years of simulation. For example, if a patient $k_7 = 18.2$ and other parameters are fixed as in paper ^[38], then the HT model would predict the blood TSH concentration as $[TSH] = 7.3$ mU/L after $4$ years of simulation. Next, one can convert the blood concentration of TSH into molecules with a conversion formula stated below ^[40,41,42]:

Suppose the blood concentration of $[TSH] = 7.3$ mU/L, then the number of TSH molecules is roughly $981786733$ in the blood (see Figure 4). One molecule of TSH binds with a receptor protein (TSHR) located on the surface of thyroid follicle cell ^[43] that results in TSH-TSHR pathway. Also, TSH-TSHR signaling pathway plays a critical role for thyroid cell growth and proliferation in HT. It usually works through signaling pathways that involve activation of proto-oncogenes and has a role in controlling cell growth and carcinogenesis ^[44]. In patients with HT, the number of TSH molecules distributed to receptor sites at the follicles (affected by HT) is expected to increase stimulation as patients move through various clinical states. Let $Y$ be a binomial random variable counting the number of successful first bindings (TSH-TSHR pathways) that can occur on given $m$ independent TSH molecules. So, when $m$ independent TSH molecules distributed to receptors either binds or misses that follows the binomial distribution and has the functional form:

where $p$ represents the probability of successful binding to a receptor and $q = 1-p$ represents the probability of unsuccessful binding to a receptor. The mean or expected value of a binomial random variable is $E(Y) = mp$ and the variance of Y is $V(Y) = mpq$ (for more details, see ^[45]). For instance, if $53796533$ TSH molecules are distributed to thyroid receptors and have 95% chance of binding to receptors, then there would be on average $51106706$ TSH molecules bound to TSHR; in other words $51106706$ ligand receptor pathways are successfully established for hormone production and secretion. In the event of large TSH signaling, each TSH-TSHR pathway poses a threat for mutation or cell growth. So, all ligand receptor pathways (i.e., TSH-TSHR pathways) are potential growth sites when there is an over-activation by TSH; in fact a mutation process can take place on these pathways resulting in nodule formation, TSH-TSHR $\rightarrow$ N (benign nodule).

4.   Input rate or distribution

As stated before, nodules appear in the areas of functional follicles that are undergoing immune attack. Assuming $n$ potential growth sites in the areas of follicles affected by the immune attack, consider the distribution of $m$ TSH molecules to those $n$ sites. By applying the Hill dynamics, we can determine a steady-state approximation for the distribution of $m$ TSH molecules to $n$ receptors with $m$ specific-binding sites ^[46,47]. We first consider the mechanism for multiple binding patterns of a receptor leading to the nodule formation:

Herein, $N$ represents a nodule and if $R_m$ represents a receptor with $m$ specific-binding sites and $m$ TSH molecules are able to bind to the same receptor, then the related dissociation constant $K_m$ is:

and the Hill equation relates the ratio ($\lambda$) of occupied receptor concentration to the total receptor concentration. That is,

where $0 \leq \lambda < 1$ and $m > 1$ represents a dimensionless parameter (the integer value) indicating the positive cooperativity of the TSHR receptor. Observe that $\lambda$ is the input rate for the nodule growth in terms of blood concentration of TSH and for some constant $K_m$. We further assume the total receptor concentration is directly proportional to total number of functional follicles per given unit volume of the thyroid gland. Similarly the occupied receptor concentration is directly proportional to the number of growth sites per given unit volume of the thyroid gland. So, we can extend this rate of nodule formation in terms of the follicles on the human thyroid gland.

Since there are approximately $315$ follicles per 1 square millimeter area of the human thyroid gland ^[48], we can further modify and extend the input rate into the surface area level of the gland. It has been estimated that one TSHR receptor has 3 distinct and discontinuous TSH-binding pockets ^[49,50]. Therefore, we choose $m = 4$ to represent positive cooperativity and the birth rate of nodule affected by HT can be denoted as,

Next, we denote $G(x)$ as the growth rate of existing nodule affected by HT for any given $x$ ^[51].

Here, the threshold value $x_{th} = [TSH]_{th}$ is obtained from the intersection of the cumulative birth and growth rate of nodules in the follicles. We computed a threshold value for the patient with $k_7 = 18.2$ by running the HT model for $4$ years. In the birth rate, the value of $K_m = 10^3$ was chosen arbitrarily to show the existence of the threshold (see Figure 5). So, $x_{th} \approx 5.648$ mU/L. The survival rate of nodule can be defined as $s(x) = G(x)-\lambda(x)$ at a given value of $x$. That is

The negative survival rate of nodule implies the birth rate is greater than the growth rate (see Figure 6).

5.   Construction of a stochastic model

We construct a stochastic model for the number of thyroid nodules for given TSH since there exists the natural variance in the population of $n$ growth sites at the follicle cells. The model will be built based on the following assumptions:

Assumptions

1. Elevated TSH ($k_7 > k_7^*$) initiates the birth of thyroid nodule in the follicle cells

2. Elevated TSH ($k_7 > k_7^*$) promotes the growth of existing thyroid nodule

3. There are $n$ sites located in the follicle cells where the nodule can be born and grow

4. Exactly one birth can occur in $(x, x+\Delta x)$ with probability $n\lambda \Delta x$ in a population of $n$ sites

5. Anti-thyroid antibodies do not destroy the nodules (no death)

6. There is no interaction between nodules, so they are Independent

5.1.   The initiation of thyroid nodules

Let $N(x)$ denote the discrete random variable for the number of thyroid nodules with associated probability function

where $i = n, n+1, n+2, ...$. The probability distribution of thyroid nodules is by letting

for some TSH value ($x$). Now, we relate the random variables $\{N(x)\}$ indexed by $x$ by defining the probability of a transition from state $i$ to state $j$, $i \rightarrow j$ in TSH $\Delta t$ as

The stochastic process $\{N(x)\}_{x = 1}^\infty$ satisfies the Markov property that means the process at $x+\Delta x$ only depends on the process at the previous step. That is, $prob\{N(x+\Delta x)|N(0), N(1), ..., N(x)\} = prob\{ N(x+\Delta x)|N(x)\}$. Since the length of TSH interval $\Delta x$ is sufficiently small, so that at most one birth can occur in $(x, x+\Delta x)$, that is to say, $i \rightarrow i+1$ or $i \rightarrow i$. The transition probabilities can be written in a compact form as

where $i = n, n+1, n+2, ...$ and $\lambda(x)$ represents the birth rate of thyroid nodule (by assumption 1) whereas $G(x) = 1-\lambda(x)$ the growth rate of existing nodule (by assumption 2). We constructed the state diagram using the above assumptions for nodule birth and growth (see Figure 7). The transition matrix for this Markov chain is

and the state equation can be written as $p(x+\Delta x) = A p(x)$. Upon performing the matrix multiplication gives the system of equations

Next, we let $\Delta x\rightarrow 0$ and obtain a system of differential equations. For $i = n$, we have

For $i = n+1, n+2, ...,$ we have

Equation (5.1) is the exponential differential equation with the initial probability $1$ at TSH $x = 0$. Since there are initially $n$ sites in the thyroid, this gives the solution

Equation (5.3) predicts the probability that the thyroid is still at size $n$ at TSH $x$, that is, the probability that no births occurred in the interval $[0, x]$. Using (5.3) with $i = n+1$ in (5.2), we obtain a first order linear differential equation

Equation (5.4) can be solved using the integrating factor $e^{(n+1)\lambda(x) x}$ and the initial condition $p_{n+1}(0) = 0$ for the solution

By using (5.5) in (5.2), we can find $p_{n+2}(x)$. Continuing on the recursive manner, we can obtain the general formula for the probability distribution of the number of thyroid nodules at given TSH $x$ value as,

5.2.   Average number of thyroid nodules

To calculate the expected value of number of nodules, we can use the definition of the expected value of a discrete random variable, that is,

Differentiating this expression and substitute $p'_i(x)$ gives the equation

We obtained the differential equation $E'(N(x)) = \lambda(x) E(N(x))$ with initial condition $E(N(0)) = n$. The solution to this equation is then the average number of thyroid nodules at TSH $x$,

Similarly, we can calculate the variance $E(N^2(x)) = ne^{\lambda(x) x}(e^{\lambda(x) x}-1)$ to provide a measure for the spread of the nodules.

5.3.   Numerical simulation

For all simulation of the HT model, we use the same initial condition (1, 13, 0.015, 1) and parameter values from Table 1. Notice that Table 1 does not include the value of $k_7$ as it varies between individuals plus it has been used mainly to provide the clinical progression of Hashimoto disease. For all simulation of the stochastic model, we use the initial nodule site, $n = 1$. Using Eqs (5.6) and (5.7), our goal in this section is to produce the statistical distribution of number of thyroid nodules and the mean number of thyroid nodules for given TSH value at each of the clinical states, euthyrodisim, subclinical and overt hypothyroidism. So, we let $k_7 = 2 < k_7^*$ that corresponds to euthyroidism state of an individual, then run the HT model to obtain TSH output for 4 years with time step-size of $1$ week. The output of TSH did not change from the initial value 1 mU/L, so inputting this value into the stochastic model with initial number of thyroid nodule, $n = 1$ predicts its probability distribution (see Figure 8).

Next, we let $k_7 = 10 > k_7^*$ in the HT model, which corresponds to subclinical hypothyroidism state of an individual, then simulate the model for 4 years with time step-size of $1$ week and obtain the output of TSH values (see Figure 9). Using these TSH values as input, we obtained the birth, growth and survival rate of thyroid nodules (see Figure 10). With all this information, the stochastic model can predict the probability distribution for number of thyroid nodules for various TSH values = 1 mU/L, 1.0142 mU/L, 1.5224 mU/L and 2.4116 mU/L respectively (see Figure 11) and the average number of thyroid nodules as TSH value varies from 1 to 4 mU/L (see Figure 12). Lastly, we let $k_7 = 18.2 > k_7^*$ in the HT model, which corresponds to overt hypothyroidism state of an individual, then simulate the model for 4 years with time step-size of $1$ week and obtain the output of TSH (see Figure 4). For this case, we already obtained the birth and growth rates with the threshold value $x_{th} \approx 5.648$ using TSH output as input (see Figure 5). So, the stochastic model can predict the probability distribution for number of thyroid nodules for various given TSH value = 1 mU/L, 1.5113 mU/L, 3.633 mU/L and 5.1896 mU/L respectively (see Figure 13) and the average number of nodules for given TSH value varies from 1 to 7.2 mU/L (see Figure 14).

6.   Growth of the largest thyroid nodule

After the birth of a nodule, we can focus on its size and growth at areas of the thyroid gland affected by HT. The size of largest thyroid nodule in terms of volume has the reference range approximately (0.27–1.71) mL ^[19]. The size of largest thyroid nodule depends on the functional size of the thyroid gland and its rate of birth. We can compute the size of largest thyroid nodule (in mL) denoted by $L(x, z)$ from the following distribution.

where $z$ is the functional size of the thyroid gland (in mL) can be obtained from the HT model for any autoimmune patient. The marginal probability density function can be computed for simulation by taking the derivative on both sides of Eq (6.1) with respect to $x$, that is,

Similarly, the marginal probability density function can be computed by taking the derivative on both sides of Eq (6.1) with respect to $z$, that is,

To measure change in the size of largest thyroid nodule, we analyze the nodule size distribution for given $x = x_0$ and $z$ varies or for given $z = z_0$ and $x$ varies that can be written in the following manner:

where $L{(x_0, z_0)}$ is the baseline size of the largest thyroid nodule. Significant growth has been defined as an increase of $20\%$ of the baseline size of the nodule ^[19].

6.1.   Numerical simulation of nodule growth

Using the output from the HT model in $(6.1)$, we can obtain the distribution of nodule size (in mL) for the specific TSH value or the functional size of thyroid gland, the results can be visualized through a histogram. Generally, the existing nodule is likely to grow, mature and possibly become cancerous in the event of large signals of TSH. So, our goal in this section is to generate the frequency of the nodule size with the output from the HT model for the same patients used in the previous section. As with the first case scenario, we simulate the HT model for 4 years with the value of $k_7 = 10$ and obtain the output of TSH and the functional size of the thyroid gland. By using the output of functional size of thyroid gland in (6.1), we can produce the distribution of thyroid nodule size for a given value of $x$ (see Figure 15). As we expected, the nodule size increases significantly when TSH value goes up. As with the second case scenario, we simulate the HT model for 4 years with the value of $k_7 = 18.2$ and obtain the output of TSH and the functional size of the thyroid gland. We noted again as TSH increases, the nodules become large in size which displays the significant growth that can be seen through the histogram (see Figure 17). Similarly, we can produce the distribution of thyroid nodule size for a given value of the functional thyroid size (say $z = 15$ mL) and using the output of TSH in (6.1). We observed the frequency of largest thyroid nodule size increases as the effect of increasing TSH (see Figure 16).

6.2.   The incidence function of cancer nodules

Let $x_N$ be the TSH value that a cancer nodule develops for the first time given that there are no cancer nodules at $x_{0}$. We now denote $h_{x_N}(x)$ be the probability density function (pdf) of ${x_N}$ and $\lambda_{x_N}(x)$ the cumulative distribution function (cdf) of ${x_N}$ so that $h_{x_N}(x) = \dfrac{d\lambda_{x_N}}{dx}$. Then, we define the incidence function $I(x)$ of cancer nodule onset at a given $x$ as

Using $\lambda_{x_N}(x)$ from Section (4), we have $h_{x_N}(x) = \dfrac{4K_mx^3}{(K_m+x^4)^2}$ and $s_{x_N}(x) = \dfrac{(K_m+x^4)}{K_m-x^4}$. Therefore,

For the hypothyroid patient with $k_7 = 18.2$, we plotted the incidence function for TSH value below and above the threshold value (see Figure 18). When $x < x_{th}$, the cancer incidence on the nodule increases however when $x > x_{th}$ the cancer is initiated on the nodule. Next, the probability distribution function at $x_N$ can be expressed in terms of the incidence function of the cancer nodule,

Using the formulas $(6.2)$ and $(6.3)$ for given $I(x)$, one can derive the pdf and the survival function. Moreover, we can determine the average number of cancer nodules at $x_N$ over the interval $x_0 \leq x_N \leq x$. That is,

7.   Discussion and conclusion

In this article, we studied via modeling an impact of delaying treatment in patients with Hashimoto's thyroiditis until the clinical signs and symptoms appear and become evident. Serum TSH is increasing in Hashimoto disease, which is the main clinical indicator and a substitution factor for the destruction process initiated by the immune system. In general, the patients with this disease go through one of the clinical progressions, from euthyroidism to euthyroidism, euthyroidism to subclinical hypothyroidism, or euthyroidism to subclinical hypothyroidism and then overt hypothyroidism. Basically, these are three clinical mechanisms typically seen in Hashimoto patients. A consequence of increasing serum TSH encourages the birth and growth of thyroid nodules on these patients, which can then become suspicious and malignant. So, we used serum TSH as a continuous input into the stochastic model to mimic the nodule birth and growth process when no treatment was administered to control the rise of serum TSH.

We first simulated our deterministic model of Hashimoto's thyroiditis and the output of the clinical progression depending upon the patient-specific parameter $k_7$ value. More precisely, we simulated the model for 4 years with a time step-size of 1 week and obtained the output of TSH and the functional size of the thyroid gland for specific $k_7$ value progressing to euthyroidism, subclinical hypothyroidism and overt hypothyroidism. Second, we used the weekly TSH output values to construct the input distribution of the birth and growth rate for the stochastic model, which can then produce the probability distribution for the thyroid nodules' number, size (in terms of mL) and growth for a given TSH value of a patient. The input distribution of thyroid nodule was determined through the Hill dynamics with positive cooperativity. Third, we calculated the survival rate of a thyroid nodule for the given TSH value of a patient. Finally, we obtained a patient-specific threshold or cutoff value ($x_{th}$) of TSH at the intersection of the birth and growth rate of thyroid nodules. The threshold on TSH values can play an important role in the initiation of thyroid microcarcinoma.

In summary, two mathematical models have been used to understand the role of TSH in thyroid oncogenesis and the increased risk of cancer. The models are patient-specific; it can help us study the association of Hashimoto's thyroiditis and papillary carcinoma through serum TSH and the functional size of the thyroid gland. The results involving the rate of production of TPOAb ($k_7$) and threshold on TSH can help with the clinical decision such as ordering the fine-needle aspiration cytology to evaluate the malignancy risk of a nodule. For example, when the production rate of TPOAb is above the critical value ($k_7 > k_7^{*}$) and the serum TSH value is below the threshold ($x < x_{th}$), the stochastic model can predict the distribution of thyroid nodule size and its growth. On these nodules, the incidence of thyroid carcinoma increases as serum TSH increases but there is no need for the fine-needle aspiration examination. When $k_7 > k_7^{*}$ and $x > x_{th}$, the nodules may already became cancerous and the results recommend the fine-needle aspiration examination. Before the start of hormone replacement therapy, the two models can be used for testing assumptions that exist between the thyroid carcinoma and Hashimoto's thyroiditis. The major contribution of this paper is setting up the foundation for the patient-specific TSH cutoff value for monitoring the malignancy risk of papillary carcinoma.

Acknowledgments

We would like to thank the reviewers for their comments and suggestions to improve this manuscript. We also thank the University of Wisconsin-Whitewater for their summer support of this project.

Conflict of interest

We have no conflict of interest in this paper.