May 20, 2018, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Correcting the analysis for Interaction Between the Thyroid and Pituitary Glands

From MathBio

Jump to: navigation, search


This project was based on correcting an integration error found in a recent article, "A mathematical model of pituitary–thyroid interaction to provide an insight into the nature of the thyrotropin–thyroid hormone relationship" by Melvin Leow [1]. The article was based on the regulating system for the production of thyroid hormone. This project used the original deterministic model and redid the analysis to come up with another solution. Then numerical analysis was performed to compare the original result to the modified one. In addition, a parameter fit of the data given in the paper was performed in an attempt to get approximate values for the various parameters.

The net result was minimal, as while the original solution was incorrect, a simplified solution found in the article could also be reached from the correct solution, albeit under assumptions differing from the original article.

Biological System Summary

The system under consideration

Please see the original wiki page covering this paper for a full explanation of the biological system under consideration. The basic concept is the self-regulating relation between the thyroid stimulating hormone (TSH) produced by the pituitary gland and the two hormones produced by the thyroid, L-3,5,3′-triiodothyronine (T3) and L-3,5,3′,5′-tetraiodothyronine (T4). Only one of the hormones needs to be considered in the model, as T3 is the active hormone, but the generation of T3 is mostly due to the transformation of T4 to T3. (T4 is in fact more commonly measured, and so is used to represent the level of thyroid hormones).

The model is therefore a relation between T, representing the level of thyroid hormones, and TSH. The model created by the author was a simple direct relation between T and TSH, without regard to time. It was designed just to show the interplay between the levels of the two hormones, in order to help physicians and clinicians gain a better understanding of the relationship between the two hormones that they observe in patients.

Differential equation

The variables of interest:

H: levels of TSH

T: levels of T4

k: a proportionality constant

\alpha _{0},\beta  : stimulatory factor of the thyroid. (The amount of T that the thyroid produces depends on the current level of T, in addition to the level of H, T_{{new}}=\alpha _{0}e^{{-\beta T}}H. The rate at which new thyroid hormone is produced decreases as T increases, -{\frac  {dT_{{new}}}{dT}}=\beta \ T_{{new}}, giving us the exponent).

The differential equation being solved is as follows:

-{\frac  {dH}{dT}}=kH(T+\alpha _{0}e^{{-\beta T}}H)

= k * Amount of H * [ (Ambient levels of T) + (New T being produced due to the level of H) ]

The analysis of this equation up until the integral is given in the wiki page summarizing this article.

The step at which an error occurs is when, with v={\frac  1H}, we have

{\frac  {e^{{-k\left.T^{2}\right/2}}}{H}}-{\frac  {1}{H_{0}}}=\int _{0}^{T}k\alpha _{0}e^{{-\left(\beta t+\left.\left(kt^{2}\right)\right/2\right)}}dt

This integral is unfortunately done incorrectly in the article. The steps to correcting this integral follow. First, some algebraic manipulation to get it into a more usable form.

\int _{0}^{T}k\alpha _{0}e^{{-\left(\beta t+\left.\left(kt^{2}\right)\right/2\right)}}dt=\int _{0}^{T}k\alpha _{0}e^{{-\left(\beta \left/{\sqrt  {2k}}\right.+\left({\sqrt  {k}}t\right)/{\sqrt  {2}}\right)^{2}}}e^{{\left.\beta ^{2}\right/2k}}dt

Using the definition of the error function, ErrorFunction(x) = erf(x) = {\frac  {2}{{\sqrt  {\pi }}}}\int _{0}^{x}e^{{-t^{2}}}\,dt, and using a u-substitution, we have that the integral is

\left.=k\alpha _{0}e^{{\left.\beta ^{2}\right/2k}}{\frac  {{\sqrt  {\pi }}}{2}}\left({\frac  {{\sqrt  {2}}}{{\sqrt  {k}}}}\right){\text{erf}}\left(\beta \left/{\sqrt  {2k}}\right.+\left({\sqrt  {k}}t\right)/{\sqrt  {2}}\right)\right|_{0}^{T} ={\sqrt  {{\frac  {\pi }{2}}}}{\sqrt  {k}}\alpha _{0}e^{{\left.\beta ^{2}\right/2k}}\left({\text{erf}}\left((\beta +kT)\left/{\sqrt  {2k}}\right.\right)-{\text{erf}}\left(\beta \left/{\sqrt  {2k}}\right.\right)\right)

Note: The error function is difficult to analyze by itself. It is best to use a computer numerical approximation program, in this case Wolfram Mathematica. The actual approximation used comes from the integral for the series for e^{{-x^{2}}},

erf(x)={\frac  {2}{{\sqrt  {\pi }}}}\int _{0}^{x}e^{{-t^{2}}}\,dt={\frac  {2}{{\sqrt  {\pi }}}}\int _{0}^{x}\sum _{{n=0}}^{{\infty }}{\frac  {\left(-t^{2}\right)^{n}}{n!}}dx={\frac  {2}{{\sqrt  {\pi }}}}\sum _{{n=0}}^{{\infty }}{\frac  {(-1)^{n}t^{{2n+1}}}{n!(2n+1)}}

The actual solution is then just a couple of algebra steps away.

H(T)={\frac  {\left(2e^{{-{\frac  {k*T^{2}}{2}}}}H_{0}\right)}{\left(2-\alpha _{0}e^{{{\frac  {\beta ^{2}}{2k}}}}H_{0}{\sqrt  {k}}{\sqrt  {2\pi }}{\text{ erf}}\left[{\frac  {\beta }{{\sqrt  {2}}{\sqrt  {k}}}}\right]+\alpha _{0}e^{{{\frac  {\beta ^{2}}{2k}}}}H_{0}{\sqrt  {k}}{\sqrt  {2\pi }}{\text{ erf}}\left[{\frac  {\beta +kT}{{\sqrt  {2}}{\sqrt  {k}}}}\right]\right)}}

Derivative Comparisons

The solution from the article, H={\frac  {H_{0}\beta (\beta +kT)}{e^{{kT^{2}/2}}(\beta +kT)(\beta +kH_{0}\alpha _{0})-k\beta H_{0}\alpha _{0}e^{{-\beta T}}}} , and the corrected solution above look quite different. I wanted to see if they did indeed lead to different relations between H and T. Both equations do have the property that H(T=0)=H_{0} However, if you take the derivative of both equations (the full result is quite messy, and so won’t be shown here), you find that for the original equation

{\text{Original }}{\frac  {dH}{dT}}(T=0)=-{\frac  {\alpha _{0}H_{0}^{2}k\left(\beta ^{2}+k\right)}{\beta ^{2}}}

whereas for the correct equation,

{\text{Corrected }}{\frac  {dH}{dT}}(T=0)=-\alpha _{0}H_{0}^{2}k

This agrees with the original differential equation,

-{\frac  {dH}{dT}}=kH(T+\alpha _{0}e^{{-\beta T}}H)

evaluated at T = 0. In fact, if you compare {\frac  {dH}{dT}} and -kH(T+\alpha _{0}e^{{-\beta T}}H) for the corrected equation, you'll find that they are the same.

Graphical Comparisons

The article used an approximation for their H equation, assuming that T was somewhat large. The approximate equation is

TSH\propto {\frac  {1}{e^{{k(freeT4)^{2}/2}}}}

Here is the plot created from this equation (assuming k = 1, and that the proportional relation is direct equality) used on the original wiki page.

Approx HTeqn Graph.jpg

Finding whether the new equation follows a similar pattern will show whether or not the original relation demonstrated the correct concepts. The only problem is a lack of knowledge about the various parameters. Here is a graph for what a complete guess on the parameters looks like:

ParameterGuess HTeqns.jpg

Obviously neither of the equations plotted above match up with the result from the paper. However, as was stated earlier, the parameters were only guesses. ( \alpha =1,k=1,H_{0}=10,\beta =0.1). These guesses can be improved, both by varying the parameters (and then looking at the results), and by trying to match the graph to the known possible values for T and H. The highest value given for H in the paper was above 70 mIU/L, the lowest was below measuring levels (effectively zero). The normal ranges are 0.3 - 4 mIU/L. For T4, the standard range is 10-20 pmol/L, with the extremes in the examples ranging from less than 2 to more than 64 pmol/L. So ideally, the relations will be able to intersect those values (assuming the parameters are in the correct units). Below is a graph of the two equations together for a given set of corrected parameters.

(\alpha =0.1,k=0.05,H_{0}=70,\beta =0.1). Note that the original shape is only followed to a limited degree (in favor of the graph having the proper range of values).

BetterParamGuess HTeqns.jpg

What is striking is the relative importance of \beta for the similarity of the two graphs. Below is a graph with the same parameters as the one directly above, except \beta =0.02.

BadB HTeqns.jpg

The maximum relative error between the two equations (Max {\frac  {H_{{\text{corrected}}}-H_{{\text{original}}}}{H_{{\text{corrected}}}}}) for \beta =0.1 is 15.4% (at T = 8.47), whereas the maximum relative error for \beta =0.02 is 57.8% (at T = 3.42)

Parameter Fit

So far we have some evidence as to the paper possibly matching up with the corrected form. It all depends on the parameters. The paper does give two examples for patients. One patient was being treated with varying levels of antithyroid drugs. For the second patient, the author uses the fact that hysteresis is present in the relationship between T and TSH to analyze the data, without telling fully what the model with hysteresis would be. These patients gave data that were excellent examples for their particular situation, but the information contained in their data made for a very poor fit for the model given above. Wolfram Mathematica was again used. The results were not pretty, but they showed a little information. In particular, the value for b was estimated as at least 0.16, which we already know leads to relatively small errors between the two forms of the model. Below is a graph for one of the fits.

DataFit HTeqns.jpg

For this fit, the parameter values found were

\alpha =0.000371,k=0.0715,H_{0}=77.5,\beta =1.86

The Approximation

The reason for the similarities between the two graphs can be partially explained by the behavior of the error function. As x becomes large, erf(x) quickly approaches 1. (For instance, erf(2) = 0.995, erf(3) = 0.99998. Thus for large values of {\frac  \beta {{\sqrt  {2k}}}},

{\text{erf}}\left({\frac  {\beta }{{\sqrt  {2k}}}}\right)\approx {\text{erf}}\left({\frac  {\beta +kT}{{\sqrt  {2k}}}}\right), and so

H(T)={\frac  {\left(2e^{{-{\frac  {k*T^{2}}{2}}}}H_{0}\right)}{\left(2-\alpha _{0}e^{{{\frac  {\beta ^{2}}{2k}}}}H_{0}{\sqrt  {k}}{\sqrt  {2\pi }}{\text{ erf}}\left[{\frac  {\beta }{{\sqrt  {2k}}}}\right]+ae^{{{\frac  {\beta ^{2}}{2k}}}}H_{0}{\sqrt  {k}}{\sqrt  {2\pi }}{\text{ erf}}\left[{\frac  {\beta +kT}{{\sqrt  {2k}}}}\right]\right)}} \approx {\frac  {\left(2e^{{-{\frac  {k*T^{2}}{2}}}}H_{0}\right)}{2}}=c*e^{{-{\frac  {k*T^{2}}{2}}}} which was the original approximation given for H(T).

The original justification for H(T)\approx c*e^{{-{\frac  {k*T^{2}}{2}}}} was for what happens when T is moderately large. The math works out, however this justification has some possible flaws, depending on for what values of T this is a good assumption. Placing too strong of an assumption on the values of T would inherently limits the model to cases without low levels of T (hypothyroidism). Of course, the same could be said for placing restrictions on the values of \beta without knowing what the actual values of \beta might be.

Still, looking at the original equation found in the article, H={\frac  {H_{0}\beta (\beta +kT)}{e^{{kT^{2}/2}}(\beta +kT)(\beta +kH_{0}\alpha _{0})-k\beta H_{0}\alpha _{0}e^{{-\beta T}}}}, we can see that for \beta >>kH_{0}\alpha _{0}, the model again reduces to the simplified form above. H\approx {\frac  {H_{0}\beta (\beta +kT)}{e^{{kT^{2}/2}}(\beta +kT)(\beta )-k\beta H_{0}\alpha _{0}e^{{-\beta T}}}}={\frac  {H_{0}}{e^{{kT^{2}/2}}-{\frac  {k\beta H_{0}\alpha _{0}}{(\beta +kT)(\beta )}}e^{{-\beta T}}}}\approx c*e^{{-{\frac  {k*T^{2}}{2}}}}

Thus, given large \beta , the results from the original article hold.

One should also note that H(T)=H_{0}e^{{-{\frac  {k*T^{2}}{2}}}} is simply the solution to the differential equation {\frac  {dH}{dT}}=-kHT,H(0)=H_{0}, which is what the original differential equation approaches as \beta becomes large.

Graves' Disease

Grave's disease is a case of hyperthyroidism wherein an additional hormone, TSHRAb (TSH receptor autoantibody), stimulates the production of thyroid hormone. The level of thyroid hormone of concern in our equation is normally T_{{total}}=T_{a}+\alpha _{0}e^{{-\beta T_{a}}}H, the level of ambient free thyroid hormone plus the amount of thyroid hormone currently being produced.

The amount of secreted hormone in Grave's disease changes in two ways. First, \alpha _{0}e^{{-\beta T_{a}}} becomes \alpha _{0}e^{{-\beta T_{a}}}+g({\text{TSHRAb}}) , where g is a constant that depends on the effects TSHRAb actually has on the thyroid's TSH receptors (which varies from patient to patient). Then H is replaced by TSHRAb.

Thus, the principle change in Grave's disease is T_{{total}}=T_{a}+({\text{TSHRAb}})*\left(\alpha _{0}e^{{-\beta T_{a}}}+g({\text{TSHRAb}})\right)

This simplifies to T_{{total}}\approx T_{a}+g({\text{TSHRAb}})^{2}, when we consider that in Grave's disease, T_{a} is going to be very large (hyperthyroidism). (Also note that g is going to be nonzero for most patients).

At this point, the paper substitutes the above equation for T_{{total}} for where T is in their simplified model, getting

H(T)={\frac  {C}{e^{{k\left(T_{a}{}^{2}+2T_{a}g({\text{TSHRAb}})^{2}+g^{2}({\text{TSHRAb}})^{4}\right)/2}}}}

But this is really an easy differential equation to solve:

-{\frac  {{\text{dH}}}{{\text{dT}}}}=kH\left(T+g({\text{TSHRAb}})^{2}\right),H(0)=H_{0}

H(T)={\frac  {H_{0}}{e^{{k\left(T+g({\text{TSHRAb}})^{2}\right){}^{2}/2}}}}

This matches up with what the article originally found.


The original article certainly had a mistake in it. However, given certain assumptions and approximations, it appears that the results remain valid. Whether or not those assumptions are justified is difficult to determine without data from healthy patients. It is also quite clear that the result that is found via those approximations is quite trivial, as using the approximations on the differential equation makes it incredibly basic.

See Also

For the Mathematica file associated with this project, which contains the code for the above plots and additional plots (static and dynamic), click here.

For the wiki page summarizing the article, click here.


  1. Leow, M., 2007. A mathematical model of pituitary–thyroid interaction to provide an insight into the nature of the thyrotropin–thyroid hormone relationship. Journal of Theoretical Biology 248, 275-287.