Sensitivity Analysis in Mathematical Models of the

Hypothalamus-Pituitary-Thyroid Axis

CLARA HORVATH, ANDREAS KÖRNER

Intitute of Analysis and Scientific Computing

TU Wien

Wiedner Hauptstraße 8

Vienna, AUSTRIA

Abstract: Mathematical models are promising and important for advancing the current medical practice in the field

of endocrinology. To assess the reliability of the range of mathematical models describing the hypothalamus-

pituitary-thyroid axis and to establish their applicability in clinical decision support, we conducted a local and

global sensitivity analysis of the model. Thyroid regulation in euthyroid and diseased individuals may be quan-

tified and dynamic behavior predicted through mathematical models, thereby revolutionizing the current clinical

practice. We investigated the influence of model parameters of a selected mathematical model utilizing ordinary

differential equations describing the HPT-axis. Motivated by a graphical depiction of the varying influence of the

model parameters, feasible methods such as a local sensitivity analysis are conducted. Furthermore, to account

for the influence of parameters on the output variance of the considered model, the theory of Sobol’ indices is uti-

lized. Although the system of differential equations describing the hormone concentrations of thyroid-simulating

hormones and unbound Thyroxine has similar equation structures, the results of the sensitivity analyses varied

according to the equation.

Key-Words: Mathematical Modeling, Hypothalamus-Pituitary-Thyroid Axis, Local Sensitivity Analysis, Global

Sensitivity Analysis, Ordinary Differential Equations

Received: May 25, 2024. Revised: June 29, 2024. Accepted: September 12, 2024. Published: October 8, 2024.

1 Introduction

Thyroid diseases affect an estimated 200 million

people globally, with an additional 60% of those

affected believed to be undiagnosed, making this an

increasingly important area of research [1], [2]. A

promising approach to gain insight into the functional

interactions of the thyroid gland lies within the field

of mathematical modeling. Throughout the years of

research varying mathematical modeling approaches

were utilized to describe aspects of the endocrine

system such as the functionality of the thyroid gland.

In general, the feedback loop of hormonal signals

starts in the hypothalamus, continues to the pituitary

gland and terminates in a specific target gland. In

the case of the thyroid gland, models consisting of

systems of differential equations were established

to describe this process. A prominent model was

introduced by DiStefano et al. [3], which describes

the course of thyroid hormone concentrations over

time using a large system of differential equations.

Additionally, a sub-model was incorporated to ex-

plore the effects of oral medication in the case of

hypothyroidism, a condition in which the thyroid

gland does not produce enough thyroid hormones.

Berberich et al. [4] investigated how their model

accounts for the influence of the circadian rhythm on

the functionality of the thyroid gland. Furthermore,

the interaction between various organs, such as liver

and heart are also included to better understand

hormone homeostasis. Making [4] a rare exception,

a sensitivity analysis was conducted as well. A very

recent contribution includes Pompa et al. [5] with

their work on compartment models. The function-

ality of the thyroid gland is described by including

concentrations of thyrotropin-releasing hormone

and thyroglobulin with additionally considering

the effect of iodine in the modeling process. Also

utilizing differential equations, Leow [6] proposed

an exponential connection between hormones of the

pituitary and thyroid gland. Utilizing approaches

from control theory, Sharma et al. [7] introduce an

approach to introduce an automated dosage recom-

mendation of levothyroxine for patients suffering

from Hashimoto’s Thyroiditis.

However, thorough analysis of the mathematical

properties of these models, although crucial for their

clinical reliability, is missing. To better understand

the underlying properties of the considered mathe-

matical models, we conducted a sensitivity analysis of

a model introduced in [8]. Based on this analysis, the

reliability of a model and its usability in clinical de-

cision support can be established. To offer a compre-

hensive treatment of the topic, the biological mech-

anisms underlying the endocrine system and the the-

oretical background of the mathematical approaches

will be presented in the following sections.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

313

Volume 21, 2024

2 Physiological Background

The central aspect of the endocrine system and start-

ing point of signaling pathways is located in the brain

and consists of the hypothalamus and pituitary gland

leading to the stimulation of target glands such as the

thyroid gland. This sequence of stimuli will be ref-

ered to as the Hypothalamus-Pituitary-Thyroid-axis

(HPT-axis) and is responsible for maintaining overall

homeostasis of the human body. Additionally, vital

functions including brain development, metabolism,

and regulating the cardiac cycle are associated with

the HPT-axis, which makes the presence of diseases

all the more consequential. The production of thy-

roid hormones is initialized by the hypothalamus and

its release of Thyrotropin-Releasing Hormone (TRH)

in the hypophysiotropic neurons. Afterwards, TRH

is transported through blood vessels in the infundibu-

lum, the connection between the hypothalamus and

pituitary gland. Thereby, the thyrotrophs located in

the anterior pituitary are stimulated, which leads to the

production of Thyroid-Stimulating Hormones (TSH).

The thyroid gland is located above the collar and be-

low the larynx, consisting primarily of follicular cells.

TSH connects to the outside membrane of these cells

and initializes the transcription of a protein called

thyroglobulin, which after additional processes leads

to the release of the thyroid hormones Triiodothyro-

nine (T3) and Thyroxine (T4). Serum concentrations

of the thyroid hormones are regulated by a negative

feedback loop, which means that the production of

TRH and subsequently TSH will cease if a certain

threshold is reached in the blood stream. According to

[9], thyroid hormones T3 and T4 are produced in the

ratio 1:10 and occur in the blood as bound or unbound

moities. Through this connection to plasma proteins

called Thyroxin Binding Globulin (TBG) the half-life

of thyroid hormones is prolonged to approximately 1

day for T3 and 7 days for T4. About 1% is freely

distributed in the blood, therefore unbound, making

them biologically more active.

3 Mathematical Models of the

Thyroid Gland

Since Danzinger and Elmergreen pioneered the use

of non-linear differential equations to model the

negative-feedback control of the HPT-axis in 1954

[10], numerous mathematical models have been de-

veloped to better understand the complex dynamics

of this endocrine system. More recent publications in-

clude [11], [2], that account for diseases related to the

thyroid gland e.g. hypothyroidism, Grave’s disease or

Hashimoto’s thyroidism in the introduced mathemati-

cal models. Among these models, introduced models

extend the base models by accounting a more general

description of the HPT-axis.

3.1 Physiological Considerations

In the context of mathematically modeling the HPT-

axis, free T4 (FT4) is of particular importance as

this hormone is transformed through deodinase to T3

in target cells such as liver or kidneys, which is re-

sponsible for the transcription of hormone-dependent

genes. Unlike TSH, TRH is not considered in math-

ematical models as acquiring data for validation is

hardly possible without invasive procedures. Thus,

the Hypothalamus-Pituitary complex (HP-complex)

is considered as single compartment and represented

by the concentration of TSH in the mathematical

models, see Fig. 1.

Fig. 1: Illustration of the aspects considered in mod-

eling the HPT-axis.

Furthermore, the functionality of the thyroid gland

is represented by the hormone concentration of FT4.

There are different opinions in the scientific commu-

nity about the normal range of hormone concentra-

tions of TSH and FT4 in blood samples [12]. The

considered normal ranges by the scientific commu-

nity vary only slightly in the case of both hormones,

therefore we choose the following as normal ranges

given in Table 1.

Hormone Concentration Normal range Unit

TSH 0.4−4mU/L

FT4 7−18 pg/mL

Table 1: Normal ranges and units of TSH and FT4.

Based on the model structure in Fig. 1, the follow-

ing two models are considered and further analyzed

through a sensitivity analysis.

3.2 Underlying Model

The first model we will consider, consists of four dif-

ferential equations and describes the negative feed-

back loop of the HPT-axis in the case of the au-

toimmune disease Hashimoto’s thyroiditis [13]. This

disease is characterized by the presence of thyroid

antibodies in blood samples, which attack and de-

stroy the follicular cells of the thyroid gland and

herefore cause hypothyroidism. According to [14], [15],

Hashimoto’s thyroiditis is the most common cause of

hypothyroidism, affecting approximately 30% of pa-

tients with thyroid disorders. Patients with hypothy-

roidism can experience fluctuating body weight, fa-

tigue, weakness and mood swings, which can signifi-

cantly reduce their quality of life.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

314

Volume 21, 2024

The model presented in [13] aims to incorporate

significant characteristics of the disease Hashimoto’s

thyroiditis into the model the HPT-axis, which focus

on the following assumptions stated in [13].

1. Anti-thyroid peroxidase antibodies recognize

and attack the functional thyroid gland (precisely

thyroid peroxidase (TPO) enzymes in the follic-

ular cells).

2. The damaged part of the gland is no longer func-

tional (active) in secreting thyroid hormones.

3. TSH stimulates the functional (active) part of the

thyroid gland for growth and hormonal secretion.

4. TSH disappears from the blood through a non-

specific excretion mechanism.

5. TSH distributes uniformly throughout the func-

tional part of the gland.

6. The hypothalamus–pituitary function is intact.

7. The blood concentration of iodine is sufficient for

synthesis of hormonal production.

8. The total TSH receptor concentration in the

thyroid gland is unaffected in the presence of

Hashimoto’s thyroiditis. [16]

Based on this the system of differential equations is

given by

dt=k1−k1y

ka+y−k2x, x(t0) = x0,

dt=(k3z)x

kd+x−k4y, y(t0) = y0,

dt=k5(x

z−N)−k6zw, z(t0) = z0,

dt=k7zw −k8w, w(t0) = w0,

(1)

with initial conditions x0, y0, z0, w0≥0at time t0is

introduced. The variables x, y, w represent the hor-

mone concentrations of TSH, FT4 and Anti-Thyroid

Peroxidase Antibodies (TPOAb) at time tin days re-

spectively. Furthermore, zaccounts for the functional

(active) size of the thyroid gland since Hashimoto’s

thyroiditis reduces the size of the thyroid gland. The

model also includes the parameters ki,N∈

i∈ {1, . . . , 8, a, d}. These parameters were derived

from literature and simulation studies. A stability

analysis of the system of differential equations (1)

has already been carried out in [13]. Furthermore,

through bifurcation a parameter is identified which in-

dicates the transition from euthyroid states to disease

states.

3.3 Reduced Model

The second mathematical model, derived from (1) and

introduced in [8], aims to describe the HPT-axis in a

healthy (euthyroid) state, without considering the ef-

fects of specific thyroid disorders. The original sys-

tem of four differential equations (1) is reduced by

omitting the equations in zand wresulting in the sys-

tem of two differential equations given by

dt=k1−k1y

ka+y−k2x, x(t0) = x0,(2)

dt=k3x

kd+x−k4y, y(t0) = y0,(3)

with initial conditions (x0, y0)Twith x0, y0≥0.

The variables xand ykeep their original meaning,

and the model parameters get reduced to k1,k2,k3,

k4,ka,kd∈

+. Since the variable zis no longer

considered, only the term k3z

kd+xin the second equation

changes, whereas the rest remains unchanged. This

is compensated by k3.

In [17], the reduced model has been the subject

of much investigation as to whether it can be used

to present clinical data through calibration of selected

model parameters. Furthermore, the above model is

similar to the differential equations introduced in [2],

however, it does not include time-dependent model

parameters. Similar to the analysis in [2], we perform

a sensitivity analysis using the following mathemati-

cal theory.

4 Sensitivity Analysis

As mathematical models become increasingly com-

plex, it becomes more challenging to understand how

small errors or perturbations in the input data affect

the simulation results produced by these models. The

theory of sensitivity analysis deals with the changes

in simulation results due to the variation of model pa-

rameters, which grants essential information regard-

ing model development and design optimization [18].

Particularly in the context of calibrating model pa-

rameters to fit a mathematical model to real-world

data, identifying model parameters that have a signif-

icant impact on model output is a critical step in the

optimization process.Therefore, the theoretical back-

ground to conduct a local and global sensitivity anal-

ysis is introduced.

4.1 Local Sensitivity Analysis

The underlying approach of this subsection is based

on [19] and introduces the theory of local sensitiv-

ity analysis, which refers to calculating derivatives of

a solution to a differential equation with respect to

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

315

Volume 21, 2024

model parameters. The following representation of

a system of differential equations is given by

∂z(t, θ)

∂t =f(t, z(t, θ), θ),(4)

with the function z:I×G→

n, z =z(t, θ)for

I⊆

, G ⊆

mopen and f:I×

n×G→

ncon-

tinuous and locally Lipschitz continuous with respect

to z. The solution zand additionally the function f

are dependent on the vector θ∈G, which represents

the model parameters. Furthermore, we assume that

the solution zis partially differentiable in θ, which

allows the approximation

zk(t, θ)−zk(t, θ∗)≈

∑

j=1

∂zk(t, θ∗)

∂θj

(θj−θ∗

j)(5)

for a fixed parameter vector θ∗and k= 1, . . . , n

using Taylor’s theorem. Thereby, the local sensitivity

of the state vector zwith respect to model parameters

θjfor j∈ {1, . . . , m}at time tis defined as

∂zk(t, θ∗)

∂θj

,(6)

where θ∗is a given parameter. According to [19], a

simple approach to calculate model sensitivities is the

use of numerical differentiation given by

∂z(t, θ)

∂θj

=z(t, θ + ∆θj)−z(t, θ)

∆θj

+O(∆θj)(7)

for j= 1, . . . , m, where θ+ ∆θjstands for the

mathematical expression of only adding the value

∆θjto the j-th component of θ. However, this

method has some drawbacks as the accuracy of the

approximation (7) strongly depends on the choice

of ∆θj. As mentioned in [20], if ∆θjis chosen

unnecessarily large, the corresponding error term of

the approximation becomes indecently large as well.

If ∆θjis set too small, the error can increase, as well

due to floating point cancellations in the calculation.

To avoid these errors, Matlab toolbox for Sensi-

tivity Analysis for ODEs and DAEs is used, as it em-

ploys a suitable strategy described in [20] for the cal-

culation of ∆θj. The Matlab function sens_sys is

used to compute the local sensitivities defined in (7).

In the case the a local sensitivity analysis, only in-

formation about the influence of a single parameter

is provided while the remaining parameters are fixed.

To gain additional insight into the collective influence

of parameters the global sensitivity analysis based on

Sobol’ indices is introduced in the following section.

4.2 Global Sensitivity Analysis

In contrast to local sensitivity, global sensitivity pro-

vides information about the collective influence of

multiple model parameters. In that regard, the method

of Sobol’ is used, which conducts a global sensitivity

analysis on the basis of variance decomposition, see

[21], [22]. We consider a rudimentary modeling ap-

proach and assume, that each model can be described

Y=f(X) = f(X1, . . . , Xn),(8)

where f:

n→

, and Yrepresents the model out-

put, which serves the simple propose of introducing

the theory assumed to be a scalar, but higher dimen-

sional output is permitted as well. Regarding the input

X= (X1, . . . , Xn)∈[0,1]n, elements within the n-

dimensional unit hypercube are considered, where the

values (X1, . . . , Xn)are uniformly distributed within

[0,1]n. This condition may be assumed without loss

of generality since any input space can be transformed

to [0,1]n. Furthermore, we assume that the Hoeffding

decomposition for the model output Yholds, and is

given by

Y=f0+

∑

i=1

fi(Xi) + . . .

∑

i=1

∑

j=i+1

fij (Xi, Xj) + . . .

+f1,...,n(X1, . . . , Xn).

(9)

A necessary condition for (9) to be fulfilled, is

given by

∫1

fi1,i2,...,is(Xi1, Xi2, . . . , Xis)dXik= 0,(10)

for 1≤i1< i2<· · · < is≤kand

ik∈ {i1, i2, . . . , is}. This results in f0=C, where

Cis a constant. Additionally, (9) is unique and the

summands are mutually orthogonal. The functions

fi1,i2,...,isare obtained by

f0=E(Y), fi=EX∼i(Y|Xi)−E(Y),

fij =EX∼ij (Y|Xi, Xj)−fi−fj−E(Y),(11)

where higher-order terms follow the same pattern.

The expression X∼irefers to the set of all variables

except Xi. Using the expressions defined in (11), the

partial variances are denoted to

Vi=V(fi(Xi)) = VXi(EX∼i(Y|Xi)),

Vij =V(fij (Xi, Xj)),(12)

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

316

Volume 21, 2024

in which case again higher order expression are ob-

tained by following the same pattern. Based on (12)

and assumption (9) an expression for the variance de-

composition is achieved as

V(Y) =

∑

i=1

Vi+

n−1

∑

i=1

∑

j=i+1

Vij +· · · +V1,...,n.

Through the use of Sobol’ sensitivity indices Si, Sij

the variance contributions to the total variance of indi-

vidual parameters and the variance accounting for the

interaction between parameters are determined. The

Sobol’ indices are given by

Si=Vi

V(Y), i = 1, . . . , n

Sij =Vij

V(Y),1≤i < j ≤n

(13)

which denotes the ratio to of the corresponding partial

variance to the total variance and are determined us-

ing a Python script. The calculation of Sobol’ indices

make it necessary to calculate a large variety of dif-

ferent numerical solutions for varying sets of param-

eters. Thereto, Python libraries such as SALiB and

SciPy including functions like saltelli.sample

and solve_ivp, are used.

5 Results

This section is based on [23] and presents the results

of the local and global sensitivity analysis, starting

with a graphical motivation. Understanding the ef-

fects of changes in the output of (2) and (3) due to

variation in model parameters provides useful infor-

mation about the extent to which the model is able to

represent clinical data. First, a graphical motivation

is given to investigate the parameter influence on a

fundamental level.

5.1 Graphical Analysis

In this approach, the model parameters are changed

individually and the corresponding solutions are

plotted. The reduced model (2) and (3) is considered

with initial conditions (x0, y0)T= (1,13)Tfor the

starting concentrations of TSH and FT4 respectively.

They are chosen in accordance with [13] and are

situated within the normal ranges of the respective

hormone. In the course of this graphical motivation

the initial conditions of the differential equations (2)

and (3) are not varied, only the underlying model

parameters. The initial values of the model parame-

ters ki∈

+, i ∈ {1,2,3,4, a, d}are assigned to

k1= 5000, k2= 16.63, k3= 1.1262, k4= 0.099,

ka= 0.0434, kd= 0.0021 and taken from literature

sources in [13] or in the case of the parameter k3

results from [17] were taken into account. Since (2)

and (3) originate from a larger system of differential

equations, the omitted quantity of the thyroid size

denoted as zin (3) is compensated by the value of

k3. This consequence has been thoroughly discussed

in [17] through the calibration of model parameters

including k3to clinical data. To conduct a graphical

presentation of the solutions of (2) and (3), the

parameters are individually increased by a factor

s∈[0.1,1000]. For each set of model parameters,

(2) and (3) are solved using the Matlab function

ode15s. Each model parameter is increased twice,

therefore resulting in three different solutions for

each hormone concentration and each parameter.

The resulting curves are displayed in Fig. 2 and

Fig. 3.

The solution associated with the starting parame-

ter set is represented by a continuous line, and those

resulting from modified parameters by a dashed line.

An increase of the model parameters k1, k4, ka, kd

leads to an increase of the resulting solutions for TSH.

The exact opposite can be observed for the model pa-

rameters k2, k3. Even though the model parameters

k3, k4, kdare only contained in (3), their increase re-

sults in a noticeable change in the solutions of (2).

0 10 20 30 40 50 60

k1= 5000

k1= 7000

k1= 9000

0 10 20 30 40 50 60

k2= 16:63

k2= 36:63

k2= 56:63

0 10 20 30 40 50 60

k3= 1:1262

k3= 1:3762

k3= 1:6262

0 10 20 30 40 50 60

k4= 0:099

k4= 0:109

k4= 0:119

0 10 20 30 40 50 60

ka= 0:0434

ka= 0:1434

ka= 0:2434

0 10 20 30 40 50 60

kd= 0:0021

kd= 0:1021

kd= 0:2021

Time

Hormone Concentration

Fig. 2: Solutions of (2) with singularly varied model

parameters ki, i ∈ {1,2,3,4, a, d}. The default

values are k1= 5000, k2= 16.63, k3= 1.1262,

k4= 0.099, ka= 0.0434, kd= 0.0021.

By comparison to the above observation, the

model parameters k1, k2, kaare only included in (2)

and their alteration lead to no changes in the solutions

of (3). This lies in contrast to the model parameters

k3, k4, kd, which influence the solution of both hor-

mone concentrations, whereas the exact opposite can

be observed for k4and ka.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

317

Volume 21, 2024

0 10 20 30 40 50 60

k1= 5000

k1= 7000

k1= 9000

0 10 20 30 40 50 60

k2= 16:63

k2= 36:63

k2= 56:63

0 10 20 30 40 50 60

k3= 1:1262

k3= 1:3762

k3= 1:6262

0 10 20 30 40 50 60

k4= 0:099

k4= 0:109

k4= 0:119

0 10 20 30 40 50 60

ka= 0:0434

ka= 0:1434

ka= 0:2434

0 10 20 30 40 50 60

kd= 0:0021

kd= 0:1021

kd= 0:2021

Time

Hormone Concentration

Fig. 3: Solutions of (3) with singularly varied model

parameters ki, i ∈ {1,2,3,4, a, d}. The default

values are k1= 5000, k2= 16.63, k3= 1.1262,

k4= 0.099, ka= 0.0434, kd= 0.0021.

A first indication of the influence of model param-

eters on model output is provided by the graphical

representation of different solutions with various un-

derlying parameter sets.

5.2 Results of Local Sensitivity Analysis

To support the conclusions drawn from the previous

section, a local sensitivity analysis is conducted. The

Matlab function sens_sys included in the toolbox

Sensitivity Analysis is utilized to calculate the local

sensitivity introduced in (6) of the reduced model.

The resulting local sensitivities ∂x

∂kjand ∂y

∂kjfor

j∈ {1,2,3,4, a, d}are shown in Fig. 4 and the fol-

lowing influences on the solution are deduced.

Starting with the model parameters k1and ka, the

local sensitivities ∂x

∂k1,∂y

∂k1,∂x

∂kaand ∂y

∂kashow the

same influence on the model output as indicated in

the graphical motivation. We observe that the indi-

vidual variation of the model parameters k1and ka

only lead to changes in the trajectories representing

the hormone TSH whereas no change is observable

in the trajectory representing FT4. Continuing with

the model parameter k2, similar behavior to k1and

kais noticeable. No change regarding the hormone

concentration of FT4 is deduced from ∂y

∂k2and an in-

crease in the parameter value leads to a decrease of

the trajectory representing TSH. In the case of k4and

kd, there is almost no influence on the resulting so-

lution of TSH. However, the individual variation of

these parameters leads to a decrease in the trajectory

representing FT4. When k4is varied, the local sen-

sitivity ∂y

∂k4even indicates a significant change in the

model output.

0 20 40 60

#10-4 k1

@k1

0 20 40 60

-0.06

-0.04

-0.02

@k2

0 20 40 60

@k3

0 20 40 60

-100

-50

@k4

0 20 40 60

@ka

0 20 40 60

-10

-5

@kd

Time

Fig. 4: Local sensitivities for all parameters

k1, k2, k3, k4, ka, kdincluded in (2) and (3). The sen-

sitivities for the differential equations (2) and (3) rep-

resenting the concentration of TSH and FT4 are rep-

resented in blue and red, respectively.

Lastly, similar to k4and kathe variation of k3does

not alter the trajectory representing TSH, however an

increase of the parameter value leads to an increase in

the solution of FT4.

5.3 Results of Global Sensitivity Analysis

The goal of conducting a global sensitivity analysis

is to quantify the relative importance of each input

variable on the value of an assigned output response

[24]. First and second order Sobol’ indices are

calculated to quantify the contribution of each model

parameter, or a combination of them, to the response

variance decomposition [25].

The first step is to decide which parameters to in-

clude in the Sobol’ analysis. The model parameters

k1, k2and k4are not included in the analysis as they

are attributed to physiological properties such as se-

cretion and excretion rate, therefore only the parame-

ters k3, ka, kdare included. Additionally, the contri-

bution of the initial conditions (x0, y0)Tis accounted

for in the following global sensitivity analysis. In

[17], the first measurement of a times series of clinical

data were used as initial values for subsequent simula-

tions of (2), (3), which varied widely. For this reason

it is important to understand how the different initial

values will affect the further course of the simulation.

For the five selected parameters k3, ka, kd, x0, y0the

following intervals for variation given in Table 2 are

chosen due to the results of the graphical motivation

and the respective normal ranges for the concentra-

tions of TSH and FT4.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

318

Volume 21, 2024

Parameter Variation Intervall

k3[1,2]

ka[0.01,2]

kd[0.001,2]

x0[7,18]

y0[2.5,4]

Table 2: Intervals considered for variation of the pa-

rameters k3, ka, kdand initial conditions (x0, y0)Tin

the global sensitivity analysis.

Next, the Python function saltelli.sample is

utilized to generate a certain amount of different sets

of parameters k3, ka, kd, x0, y0. This amount is de-

termined by N·(2D+2), where Nis a chosen value,

in this case N= 512 and Dstands for the num-

ber of parameters considered in the analysis, therefore

D= 5. This results in a total of 6144 different sets

of parameters. To generate a numerical solution of

(2) and (3), the Python function solve_ivp is used to

solve the mathematical model for each set of param-

eters generated by the function saltelli.sample.

For the solution of the mathematical model a time

interval of t∈[0,1000], unit in days, but only the

last evaluation point is considered. The resulting first

and second order Sobol’ indices, determined using

the function sobol.analyze, are displayed in a heat

map, see Fig. 5.

Fig. 5: Global sensitivities for model parameters

k3, ka, kdand initial values (x0, y0)Tregarding (2)

and (3). The Sobol’ indices are displayed in a

heatmap, the left graphic displays the results for TSH

and the right for FT4.

The Sobol’ indices for (2) and (3) are displayed

in the left and right graphic of Fig. 5, respectively.

We start with discussing the Sobol’ indices for (2).

Interestingly, the variation of the initial conditions

(x0, y0)Tin combination with any other considered

parameter does not influence the output variance

of the mathematical model in the case of TSH. In

addition, the influence of the parameter k3is rather

insignificant in comparison to kaand kd. The first

order Sobol’ indices, displayed generally on the

diagonal, of kaand kdstand out in particular. To a

small extent, the combined variation of these two

parameters lead to a variation of the output variance

of (2).

The results of the global sensitivity analysis in the

case of (3) representing the hormone concentration of

FT4 displayed on the right hand side of Fig. 5 differ

greatly from the results for TSH. In combination

with the model parameter kd, both initial conditions

(x0, y0)Tinfluence the response variance decomposi-

tion. Furthermore, the single variation of kdresults in

a significant alteration of the output model variance,

which is concluded from the first order Sobol’ index.

Even in combination with k3influential effects are

observable.

Additionally, some further numerical examples re-

garding the global sensitivity analyses are presented

in Fig. 6.

0 20 40 60

time

90.0% region

median value

model simulation

0 20 40 60

time

90.0% region

median value

model simulation

Fig. 6: Mean model response is displayed by the or-

ange line, the shaded region on blue represents 90%

of the simulation results and the dotted line additional

simulation results.

Utilizing the Matlab function sbiomodel, the re-

duced mathematical model is implemented as Sim-

Biology Model and afterwards the first- and to-

talorder Sobol indices are calculated using the func-

tion sbiosobol. These numerical results show that

even a small variation of the parameters k3, ka, kdin

accordance with Table 2 lead to a significant variation

of the model output.

6 Discussion

The primary objective of this scientific research

comprises a local and global sensitivity analysis of a

system of differential equations aiming to model the

dynamics of the HPT-axis. The conducted analysis

revealed key model parameters which strongly influ-

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

319

Volume 21, 2024

ence the model output. The hormone concentration of

TSH is described through secretion minus excretion

rate, where the secretion rate is described through

two terms. The first term k1represents the maximum

secretion rate of TSH and the second term k1y

ka+y

accounts for the inhibition rate of FT4 present in the

bloodstream using Michaelis-Menten kinetics [26].

The result that k1is a highly influential parameter

is very plausible from a physiological point of view.

Furthermore, variation of the model parameter ka,

which is also included in the term attributed to

Michaelis-Menten kinetics, led to a change in the

model output. Interestingly, the excretion rate k2has

no effect on the model output in the case of TSH.

Future research is necessary to relate the finding

to physiological characteristics of the HPT-axis.

According to the global sensitivity analysis and the

resulting first order Sobol’ indices, kaand kdequally

influence the output variance of the mathematical

model in the case of TSH, even though kdis only

included in the equation of FT4. This finding is to

a certain extent in contradiction with the results of

the local sensitivity analysis. The discrepancy may

arise due to an insufficient sample number for the

calculation of the Sobol’ indices.

In the case of FT4, the hormone concentration

is described using two terms, the secretion minus

the excretion rate. In contrast to the excretion rate

of TSH, the results of the local sensitivity analysis

suggest that the excretion rate of FT4 k4is highly

influential on the model output regarding FT4. In

addition, the model parameters that are included in

the Michaelis-Menten term also show an influence to

a certain extent. This result aligns with the findings of

the global sensitivity analysis, as the Sobol’ indices

suggest that the model parameter kdhas a strong

influence on the output variation in the case of FT4.

In addition, the second order Sobol’ indices of kdand

all other considered parameters except kashowed a

significant influence. Further research is required to

account for the influence of varying initial conditions

on model output in the case of FT4 vs. TSH and the

resulting physiological and clinical implications.

The findings concerning the model parameter

k1align with previous studies [2], which partially

also conducted a sensitivity analysis of a very

similar system of differential equations describing

the HPT-axis and incorporating thyroid diseases

e.g. Hashimoto’s thyroiditis, and Graves disease.

However, in [2] only a partial sensitivity analysis is

conducted for selected parameters, which are in addi-

tion time-dependent and dependent on the considered

thyroid disease. The identification of significant

model parameters may lead to implications for

understanding thyroid disorders. We concluded that

small changes in certain parameters result in drastic

changes in the output of the mathematical model.

By understanding the characteristics of the model

parameters, we aimed to put mathematical aspects

into context with physiological and pathological

conditions. However, in depth studies are necessary

to gain further insight. Additionally, through clinical

trails, this insight could suggest treatment strategies

by managing thyroid-related diseases by targeting

sensitive parameters. Despite the high potential

and fundamental findings regarding the sensitivity

analysis, further research is necessary to bridge the

gap between model parameters and physiological

and pathological characteristics. Additionally, the

simplifications made in the model by assuming

constant parameter values may overlook dynamic

changes in the endocrine system.

In conclusion, the analysis of mathematical mod-

els of the HPT-axis is essential to understand its regu-

latory mechanisms and connect physiological charac-

teristics with properties of the derived model. These

findings enhance our understanding of thyroid regu-

lation and pave the way for targeted therapeutic inter-

ventions. The outlook is to analyze models targeting

different hormonal axes, such as the hypothalamus-

pituitary-ovary axis, to combine different techniques

in order to accurately depict clinical data.

Acknowledgment:

The authors would like to thank their research

partner Prof. Dr. Michael Krebs from the Medical

University of Vienna for his continued guidance in

the field of endocrinology. Additionally, Lana Međo,

Alexander Edthofer and Tünde Papp for proof

reading and contributing their excellent advice.

References:

[1] J.M. Kuyl, ”The evolution of thyroid function

tests”, Journal of Endocrinology, Metabolism

and Diabetes of South Africa, Vol.20, No.2,

2015, pp.11-16.

[2] B. Yang, X. Tang, M. J. Haller, D. A. Schatz,

and L. Rong, ”A unified mathematical model of

thyroid hormone regulation and implication for

personalized treatment of thyroid disorders”,

Journal of Theoretical Biology, Vol.528, 2021.

[3] M. Eisenberg, M. Samuels, and J. J. DiStefano

III, ”Extensions, Validation, and Clinical

Applications of a Feedback Control System

Simulator of the

Hypothalamo-Pituitary-ThyroidAxis”, Thyroid :

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

320

Volume 21, 2024

Official journal of the American Thyroid

Association, Vol.18, 2008, pp.1071–1085.

[4] J. Berberich, J. W. Dietrich, R. Hoermann, and

M. A. Müller, ”Mathematical modeling of the

pituitary–thyroid feedback loop: role of a

TSH-T3-shunt and sensitivity analysis”,

Frontiers in endocrinology, Vol.9, 2018, pp.1-11.

[5] M. Pompa, et al., ”A physiological mathematical

model of the human thyroid”, Journal of

Computational Science, Vol.76, 2024.

[6] M. K. Leow, ”A mathematical model of

pituitary–thyroid interaction to provide an

insight into the nature of the thyrotropin–thyroid

hormone relationship”, Journal of Theoretical

Biology, Vol.248, No.2, 2007, pp.275-287.

[7] R. Sharma, V. Theiler-Schwetz, C. Trummer, S.

Pliz., M. Reichhartinger, ”Automatic

Levothyroxine Dosing Algorithm for Patients

Suffering from Hashimoto’s Thyroiditis”,

Bioengineering, Vol.10, No.6, 2023.

[8] S. Goede, ”General review on mathematical

HPT modeling General Review on Mathematical

Modeling in the Hypothalamus Pituitary Thyroid

System”, submitted for publication.

[9] M. Gekle, et al., Taschenlehrbuch Physiologie,

Georg Thieme Verlag, 2015.

[10] L. Danzinger, and G. L. George, ”Mathematical

theory of periodic relapsing catatonia”, Bulletin

of Mathematical Biophysics, Vol.16, 1954,

pp.15-21.

[11] M. C. Eisenberg, F. Santini, A. Marsili, A.

Pinchera, and J. J. DiStefano, ”TSH Regulation

Dynamics in Central and Extreme Primary

Hypothyroidism”, Thyroid, Vol.20, No.11, 2010,

pp.1215-1228.

[12] J. S. Neves, et. al., ”Thyroid hormones within

the normal range and cardiac function in the

general population: the epiporto study”,

European Thyroid Journal, Vol.10, No.2, 2021,

pp.150-160.

[13] B. Pandiyan, S. J. Merrill, and S. Benvenga, ”A

patient-specific model of the negative-feedback

control of the hypothalamus-pituitary-thyroid

(HPT) axis in autoimmune (Hashimoto’s)

thyroiditis”, Mathematical medicine and

biology: a journal of the IMA, Vol.31, No.3,

2014, pp.226-258.

[14] F. Ragusa, et. al., ”Hashimotos’ thyroiditis:

Epidemiology, pathogenesis, clinic and therapy”,

Best Practice & Research Clinical

Endocrinology & Metabolism Vol.33, No.6,

2019.

[15] L. Chiovato, F. Margi, A. Carlé,

”Hypothyroidism in Context: Where We’ve

Been and Where We’re Going”, Advances in

Therapy, Vol.36, No.2, 2019, pp.47-58.

[16] H. Tamaki, et. al., ”Low prevalence of

thyrotropin receptor antibody in primary

hypothyroidism in Japan”, The Journal of

Clinical Endocrinology and Metabolism Vol.71,

1990, pp.1382-1386.

[17] C. Horvath, A. Körner, and C. Modiz,

“Data-based model identification of the

hypothalamus-pituitary-thyroid complex”,

EUROSIM Congress 2023, to be published.

[18] M. Koda, A. H Dogru, and J. H. Seinfeld,

”Sensitivity analysis of partial differential

equations with application to reaction and

diffusion processes”, Journal of Computational

Physics, Vol.30, No.2, 1979, pp.259-282.

[19] C. Rackauckas, et al., ”A comparison of

automatic differentiation and continuous

sensitivity analysis for derivatives of differential

equation solutions”, in Proc. IEEE High

Performance Extreme Computing Conference

(HPEC), Waltham, 2021, pp.1-9.

[20] R. L. Burden, J. D. Faires, and A. M. Burden,

Numerical Analysis, Cengage Learning, 2015.

[21] A. Saltelli, et el., ”Variance based sensitivity

analysis of model output design and estimator

for the total sensitivity index”, Computer

Physics Communications, Vol.181, No.2, 2010,

pp.259-270.

[22] J. Nossent, P. Elsen, and W. Bauwems, ”Sobol’

sensitivity analysis of a complex environmental

model”, Environmental Modelling and Software,

Vol. 26, No.12, 2011, pp.1515-1525.

[23] C. Horvath, ”Modelling and Analysis of the

HPT-Complex”, M.Sc. thesis, Intitute of

Analysis and Scientific Computing, TU Wien,

Vienna, 2023.

[24] K. Cheng, Z. Lu, Y. Zhou, Y. Shi, and Y. Wei,

”Global sensitivity analysis using support vector

regression”, Applied Mathematical Modelling,

Vol.49, 2017, pp.587-598.

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

321

Volume 21, 2024

[25] I. M. Sobol’, ”Global sensitivity indices for

nonlinear mathematical models and their Monte

Carlo estimates”, Mathematics and Computers

in Simulation, Vol.55, No.1, 2001, pp.271-280.

[26] N. Chitnis, J. M. Hyman, and J. M. Cushing,

”Determining important parameters in the spread

of malaria through the sensitivity analysis of a

mathematical model”, Bulletin of Mathematical

Biology, Vol.70, 2008, pp.1272-1296.

Contribution of Individual Authors to the

Creation of a Scientific Article (Ghostwriting

Policy)

Clara Horvath carried out the simulations and sensi-

tivity analysis. She was also responsible for the liter-

ature review and research into the physiological con-

text to establish a connection between mathematical

models. Andreas Körner contributed with his exper-

tise in modeling and simulation, supervised the pro-

cess of writing this paper, and was responsible for

prove reading.

Sources of Funding for Research Presented in a

Scientific Article or Scientific Article Itself

No funding was received for conducting this study.

Conflicts of Interest

The authors have no conflicts of interest to

declare that are relevant to the content of this

article.

Creative Commons Attribution License 4.0

(Attribution 4.0 International , CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en

_US

WSEAS TRANSACTIONS on BIOLOGY and BIOMEDICINE

DOI: 10.37394/23208.2024.21.31

Clara Horvath, Andreas Körner

E-ISSN: 2224-2902

322

Volume 21, 2024