Modelling the Effect of the Exposure to Perfluorooctanoate (PFOA) and
Perfluorooctane Sulfonate (PFOS) on Thyroid
Homeostasis
CHONTITA RATTANAKUL
Department of Mathematics,
Faculty of Science,
Mahidol University,
Bangkok,
THAILAND
also with
Centre of Excellence in Mathematics,
MHESI,
Bangkok,
THAILAND
Abstract: - Thyroid homeostasis is crucial for the human body. The imbalance of thyroid homeostasis might
cause diseases such as hypothyroidism. Humans are exposed to PFOS/PFOA frequently since they have been
used in various industrial products. As reported that PFOS/PFOA increase the metabolic clearance rate of
thyroid hormones, we then develop a mathematical model in terms of a system of differential equations to
investigate thyroid homeostasis based on the changes in the levels of thyrotropin-releasing hormone, thyroid-
stimulating hormone and thyroid hormones when the effect of the exposure to PFOS and PFOA is also
incorporated as well. The geometric singular perturbation technique is then employed to identify the possible
dynamic behaviours obtained from the model. Numerical investigations are also presented to illustrate the
results from theoretical analysis. Both theoretical and numerical results imply that a periodic behaviour that has
been observed clinically in the pulsatile secretions of thyroid hormones, thyroid-stimulating hormone and
thyrotropin-releasing hormone could be obtained from our model. In addition, the numerical experiment also
shows that the levels of thyroid hormones and thyroid-stimulating hormone for the case when there is the effect
of exposure to PFOS and PROA are lower than those of the case when there is no effect of the exposure to
PFOS and PFOA which might lead to the imbalance of thyroid homeostasis.
Key-Words: - Thyroid hormones, thyrotropin-releasing hormone, thyroid-stimulating hormone, perfluorooctane
sulfonate, perfluorooctanoate, mathematical model
Received: November 19, 2022. Revised: July 19, 2023. Accepted: August 17, 2023. Published: September 12, 2023.
1 Introduction
Thyroid hormones (THs) secreted by the follicular
cells of the thyroid gland are necessary for the
differentiation of cells, and growth, and also
regulate significant metabolisms, [1]. The release of
thyroid hormones is regulated by the negative
feedback control on the hypothalamus and anterior
pituitary, [2]. When the circulating levels of thyroid
hormones (Triiodothyronine (T3) and Thyroxine
(T4)) are low, the secretion of thyrotropin-releasing
hormone (TRH) from the hypothalamus will be
increased, [2]. Then, the increase in TRH level
stimulates the anterior pituitary gland to produce
thyroid-stimulating hormone (TSH). The thyroid
gland is then stimulated by the increase of TSH to
produce thyroid hormones until levels of thyroid
hormones in the blood return to normal levels, [2].
Thyroid homeostasis is crucial for the human body.
The imbalance of thyroid homeostasis might cause
diseases such as hypothyroidism. Hence, the balance
of thyroid homeostasis must be kept.
Perfluorooctane sulfonate (PFOS) and
perfluorooctanoate (PFOA) are synthetic
compounds in the group of per- and polyfluoroalkyl
substances (PFASs) used in many industrial
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
272
Volume 18, 2023
applications such as firefighting foam, floor wax,
paints, food packaging, cookware, cosmetics and
textiles for decades, [3], [4]. Humans are exposed to
PFAS because they accumulate in soil and waters,
persisting for years and are difficult to eliminate,
[3]. The contamination of PFASs in surface water,
groundwater and soil has been reported frequently,
[5], [6], [7]. Therefore, if the use of PFASs is
continued, the contamination of PFASs will also
continue and accumulate in drinking water, food and
the environment which could lead to a global
environmental pollutants crisis that could also
threaten human health, [5]. Exposure to PFOA is
related to the dysfunctions of the thyroid and cancer
related to the kidney and testicular, [3]. On the other
hand, exposure to PFOS is also related to a decrease
in fertility and adverse effects on the development
of the fetus, [3]. In children, an impaired immune
response is also reported to be related to exposure to
PFOS and PFOA whereas increased levels of
cholesterol and obesity are also observed in adults,
[3]. In this paper, we investigate the effect of
exposure to PFOS/PFOA on thyroid homeostasis
using mathematical modelling which has never been
studied before. As it has been reported in, [8], PFOS
and PFOA increase the metabolic clearance rate of
thyroid hormones which affects thyroid homeostasis
by reducing the circulating levels of thyroid
hormones in diet-exposed animals, a system of
differential equations is then developed in the next
section based upon the changes in the levels of
TRH, TSH, THs when the effect of PFOS and
PFOA is incorporated.
2 Model Development and Analysis
Here, we let X, Y and Z represent the concentrations
of TRH, TSH and THs in blood at time T,
respectively. We propose the following
mathematical model in terms of a system of
differential equations to investigate the effect of
PFOS and PFOA on thyroid homeostasis.
12 1
2
12
a a XZ
dX bX
dT h X hZ
(1)
3 4 5 62
2
3 4 5
()a XY a a Y
dY a Y b Y
dT h X h Z h Y
(2)
783
6
a
dZ a Y b Z
dT h Z
(3)
where we assume the positive values for all
parameters in the system of equations (1)-(3).
Equation (1) accounts for the rate of change of
TRH’s level in blood at time T. The feedback
control of TRH on its secretion from the
hypothalamus is presented by the first term on the
right of (1) while the increase in the secretion of
TRH due to the stimulation of THs is presented by
the second term and the rate at which TRH is
removed from the system is presented by the third
term.
Equation (2) accounts for the rate of change of
TSH’s level in blood at time T. The secretion of
TSH from the anterior pituitary gland due to the
stimulating effect of TRH, the feedback control of
TSH on its secretion from the anterior pituitary
gland and the secretion of TSH due to the level of
this is presented by the first term on the right of
equation (2). The metabolic clearance rate of TSH
due to the exposure to PFOS and PFOA is presented
in the second term while the rate at which TSH is
removed from the system is presented in the third
term.
Equation (3) accounts for the rate of change of
THs’s level in blood at time T. The feedback control
of THs on its secretion from the thyroid gland is
presented by the first term on the right of equation
(3). The secretion of THs from the thyroid gland due
to the level of TSH is presented by the second term
while the rate at which THs are removed from the
system is presented by the third term.
It has been reported that the half-life of TRH and
TSH are approximately 6 minutes, [9], and 60
minutes, [10], respectively. The half-life for THs is
approximately 1 day, [11], and 6-7 days, [10], for T3
and T4, respectively. Therefore, TRH, TSH and THs
have the fastest, intermediate and slowest dynamic
behaviour, respectively, and hence, the geometric
singular perturbation technique which has been
widely used to analyze the system with different
speeds of dynamic behaviours, [12], [13], [14], [15],
[16], can be applied to analyze the system
theoretically. By scaling the variables and
parameters of the system by small positive values
and
as follows:
**
1
1 2 2
2
* * * * *
, , , , , ,
Ta
X Y Z T
x y z t c c T a
X Y Z T X
*
**6
3 4 4 5 5 6
2
**
, , , ,
Ta
T
c c a c Y a c
ZY
* * *
7 8 3
12
7 8 1 2 3
2* * * *
*, , , , ,
T a T Y a h
hh
c c k k k
Z X Z X
Z


WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
273
Volume 18, 2023
*
*
56
42
4 5 6 1 1 2
2
**
*
, , , , ,
hh
h T b
k k k d T b d
ZZ
Y
*3
3Tb
d

, the system of equations (1)-(3) becomes
121
2
12
,,
c c xz
dx d x I x y z
dt k x kz
(4)
3 4 5 62
2
3 4 5
()
,,
c xy c c y
dy c y d y
dt k x k z k y
J x y z




(5)
783
6
,,
c
dz c y d z K x y z
dt k z




(6)
This means that
x
possesses the fastest
dynamics, y possesses the intermediate dynamics
and z possesses the slowest dynamics. Next, let us
investigate the manifold
I =0
,
J =0
and
K =0
in detail.
The manifold
I = 0
By setting equation (4) equals to zero, we obtain
12
12
12
c c z
d
k x x kz

Note that this equation does not depend on y and
hence, the manifold
0I
is parallel to the y-axis.
The intersection of the manifold
0I
and the x-
axis on the
,xz
-plane occurs at the point for
which
2
1 1 1 1 1 1
1
1
40
2
d k d k c d
xx
d
and z = 0.
It also attains the relative maximum along the line
,
MM
x x z z
where
2
1 1 1 1 1 1 22
11
12
4,
22
M
Ak Ak Ac ck
x A d
Ak
and
2M
zk
. Note that
0
M
x
if and only if
10,A
that is
22
12
2
ck
dk
(7)
The manifold
J = 0
By setting equation (5) equals to zero, we obtain
y = 0 and
3 4 5 4 6 2
2
6235
1,
c x c c y
z k c d r x y
cd k x k y





Therefore, the manifold
0J
consists of two
parts. The nontrivial manifold
,z r x y
intersects
the trivial manifold y = 0 along the curve that is
asymptotic to the line
1
zz
where
3 4 4 5 6 2
1
5 6 2
c c k k c d
zk c d

.
Note that
10z
if
3 4 4 5 6 2
c c k k c d
(8)
Moreover, on the
,xz
-plane, the nontrivial
manifold
,z r x y
intersects the x-axis at the
point for which
2
xx
and z = 0 where
3 4 5 6 2
2
3 4 4 5 6 2
k k k c d
xc c k k c d

Note that
20x
if the inequality (8) holds. In
addition, the intersection of the nontrivial manifold
,z r x y
and the
,xy
-plane occurs along the
curve
2
3 4 6 2 5
2
3 4 5 4 6 2 5
k k c d k y
x w y
c c c y k c d k y


wy
has a relative minimum at the point for which
22
4 4 5 5
5
0
m
c c c k
yy c
and
2
3 4 6 2 5
2
3 4 5 4 6 2 5
m
m
mm
k k c d k y
xx c c c y k c d k y


.
Since
2
3 5 4 5 5
2
2
6 2 3 5
2c x c y c y c k
r
y c d k x ky




then,
0
r
y
at the point where
m
yy
. Moreover,
23 5 4
2
22
6 2 3 5
22
0
m
m
yy m
c x c y c
r
y c d k x ky






WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
274
Volume 18, 2023
whenever
0x
. Therefore, the nontrivial manifold
,z r x y
obtains its relative maximum at the
points for which
,m
x y y

and
3 4 5 4 6 2
2
62 35
1m
m
m
c c c y
z z k c d
cd k k y





while
is any positive constant.
The manifold
K = 0
By setting equation (6) equals to zero, we obtain
2
3 3 6 7
8 8 6
d z d k z c
y v z
c z c k


(9)
We can see that this equation does not depend on x
and hence, the manifold
y v z
is parallel to the
x-axis. In addition,
y v z
intersects the z-axis at
the point for which y = 0, and
2
zz
where
2
3 6 3 6 7 3
2
3
40
2
d k d k c d
zd

Moreover,
22
8 3 8 3 6 7 8 8 3 6
2
8 8 6
20
c d z c d k z c c c d k
vz c z c k

whenever
0z
and hence
vz
is an increasing
function and has no relative maximum or relative
minimum for
0z
.
The different dynamic behaviours of
solutions of the system of equations (4)-(6) could be
expected when the locations of the manifolds
0 , 0IJ
and
0K
are different.
Theorem 1 For sufficiently small
and
, suppose
that
21
,xx
(10)
,
Sm
yy
(11)
21
zz
(12)
and the inequality (8) holds then a limit cycle exists
and a periodic solution occurs for the system of
equations (4)-(6).
Given that the inequalities (8) and (10)-(12)
identified in Theorem 1 are satisfied, the manifolds
0 , 0 , 0I J K
will be positioned as
shown in Figure 1 (Appendix). Starting from the
point O as located in Figure 1 (Appendix), a
transition in the direction of decreasing x with the
fast speed, since
0I
here, will bring the solution
trajectory to the fast manifold
0I
where a point
P is reached. Since
0J
here, an intermediate
transition develops along the manifold
0I
in
the direction of decreasing y until the point Q on the
stable branch of
0IJ
is reached. A transition
in the direction of decreasing z at a slow speed then
develops along
0IJ
until it reaches the point
R for which the manifold loses its stability. A jump
with an intermediate speed in the direction of
increasing y from the point R will then bring the
solution trajectory to the other stable branch of
0IJ
where the point T is reached. A slow
transition will develop in the direction of increasing
z, since
0K
here, from the point T to the point U.
A transition of the solution trajectory then brings the
system to the other stable branch of
0IJ
where the point V is reached. This is followed by a
slow transition which brings the system from the
point V to the point R resulting in a closed cycle
RTUV and hence, the system of equations (4) - (6)
exhibits a limit cycle and a periodic solution then
occurs.
Theorem 2 For sufficiently small
and
, suppose
that
14m
x x x
, (13)
then the steady state
1 4 2
,0,S x z
of the system of
equations (4)-(6) is stable.
Given that the inequality (13) identified in
Theorem 2 are satisfied, the manifolds
0 , 0 , 0I J K
will be positioned as
shown in Figure 2 (Appendix). In Figure 2
(Appendix), starting from the in front of the
manifold
0I
at a point O which is above the
manifold
0J
. A fast transition in the direction
of decreasing x will bring the solution trajectory to a
point P on the fast manifold
0I
since
0I
here. Since
0J
here, the transition in the direction
of decreasing y at intermediate speed develops along
the manifold
0I
until it reaches the point Q on
stable branch of
0IJ
. Here
0K
, a
transition at a slow speed in the direction of
decreasing z then develops along this curve until it
reaches the equilibrium point
1 4 2
,0,S x z
where
0I J K
.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
275
Volume 18, 2023
Therefore, the washout equilibrium point
1 4 2
,0,S x z
is stable in this case and the solution
trajectory tends toward
1
S
as time passes.
Theorem 3 For sufficiently small
and
, suppose
that
21
xx
(14)
mS
yy
, (15)
and the inequality (8) holds then the system of
equations (4)-(6) has a stable steady state
2,,
S S S
S x y z
.
Given that the inequalities (8), (14)-(15)
identified in Theorem 3 are satisfied, the manifolds
0 , 0 , 0I J K
will be positioned as
shown in Figure 3 (Appendix). In Figure 3
(Appendix), starting from the in front of the
manifold
0I
at a point O which is above the
manifold
0J
. A fast transition in the direction
of decreasing x will bring the solution trajectory to a
point P on the fast manifold
0I
since
0I
here. A transition in the direction of decreasing y at
intermediate speed develops along the manifold
0I
to the point Q on stable branch of
0IJ
. A transition at a slow speed in the
direction of decreasing z then brings the solution
trajectory to point R in which the stability of the
manifold is lost. A jump in the direction of
increasing y at intermediate speed will then bring
the system to the point T on the other stable branch
of
0IJ
. since
0K
here, a transition in
the direction of increasing z at slow speed then
brings the solution trajectory to the point
2,,
S S S
S x y z
where
0.I J K
Therefore,
the equilibrium point
2,,
S S S
S x y z
is stable and
the solution trajectory tends toward
2
S
as time
passes.
3 Numerical Simulations
In this section, we illustrate the theoretical results by
carrying out numerical simulations for each case.
All simulations are generated by MATLAB using
the Runge-Kutta 4th-order method.
The computer simulation shown in Figure 4
(Appendix) is presented to illustrate the theoretical
prediction in Theorem 1 in which a limit cycle is
expected provided that the inequalities stated in
Theorem 1 are all satisfied.
We can see that the solution of the system of
equations (4)-(6) shown in Figure 4 (Appendix)
tends toward a limit cycle and a periodic solution
occurs as theoretically predicted in Theorem 1.
The computer simulation shown in Figure 5
(Appendix) is presented to illustrate the theoretical
prediction in Theorem 2 in which the equilibrium
point
1
S
is stable provided that the inequality stated
in Theorem 2 is satisfied.
We can see that the solution of the system of
equations (4)-(6) shown in Figure 5 (Appendix)
tends toward the equilibrium point
1
S
as
theoretically predicted in Theorem 2.
The computer simulation shown in Figure 6
(Appendix) is presented to illustrate the theoretical
prediction in Theorem 3 in which the equilibrium
point
2
S
is stable provided that the inequalities
stated in Theorem 3 are all satisfied.
We can see that the solution of the system of
equations (4)-(6) shown in Figure 6 (Appendix)
tends toward the equilibrium point
2
S
as
theoretically predicted in Theorem 3.
4 Discussion and Conclusion
We develop a mathematical model to investigate the
effect of exposure to PFOS/PFOA on thyroid
homeostasis. Theoretically, the geometric singular
perturbation technique is utilized so that we obtain
the conditions on the system parameters that
differentiate various behaviours of the solutions of
our model. Numerical simulations are provided to
illustrate the theoretical results. Both theoretical and
numerical results indicate that a periodic behaviour
that has been observed in the pulsatile secretions of
TRH, TSH and THs, [17], [18], [19], could be
exhibited by our model when the parametric values
are appropriated. In addition, let us consider an
example of the computer simulations in Figure 7
(Appendix) for the case when there is no effect of
PFOS/PFOA (
60c
) and the case when
there is the effect of PFOS/PFOA (
60.08c
) while
the other parametric values are the same that is
12
0.15, 0.30,cc
3 4 5 7 8 1
0.55, 0.80, 0.7, 0.21, 0.7, 0.95,c c c c c k
2 3 4 5 6 1
0.75, 0.9, 0.5, 0.3, 0.7, 0.55,k k k k k d
23
0.1, 0.03, 0.95, 0.035dd

with
0 0.1, 0 0.5,xy
and
02z
. We can see
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
276
Volume 18, 2023
that when there is the effect of PFOS/PFOA, the
levels of TSH and THs are lower than those of the
case when there is no effect of PFOS/PFOA which
might lead to the imbalance of thyroid homeostasis
and hence, the hypothyroidism might be expected.
The developed model might be useful to investigate
further the treatment of hypothyroidism due to
exposure to PFOS and PFOA.
Acknowledgement:
This work is supported by the Centre of Excellence
in Mathematics, Thailand and Mahidol University,
Thailand.
References:
[1] Leko, M. B., Gunjaca, I., Pleic, N., Zemunik,
T. Environmental factors affecting thyroid-
stimulating hormone and thyroid hormone
levels. Int. J. Mol. Sci., 2021. 22: 6521.
https://doi.org/10.3390/ijms22126521
[2] Dietrich, J.W., Landgrafe, G., Fotiadou, E.H.
TSH and thyrotropic agonists: key actors in
thyroid homeostasis. J. Thyroid Res., 2012.
2012: 351864. doi: 10.1155/2012/351864
[3] Coperchini, F., Croce, L., Ricci, G., Magri, F.,
Rotondi, M., Imbriani, M., Chiovato, L.
Thyroid disrupting effects of old and new
generation PFAS. Front. Endocrinol., 2021.
11: 612320. doi: 10.3389/fendo.2020.612320
[4] Buck, R.C., Franklin, J., Berger, U., Conder,
J.M., Cousins, I.T., De Voogt, P., Jensen, A.A.,
Kannan, K., Mabury, S.A., van Leeuwen,
S.P.J. Perfluoroalkyl and polyfluoroalkyl
substances in the environment: Terminology,
classification, and origins. Integr. Environ.
Assess. Manag., 2011. 7: 513541,
doi:10.1002/ieam.258.
[5] Jian, J.M., Chen, D., Han, F.J., Guo, Y., Zeng,
L., Lu, X., et al. A short review on human
exposure to and tissue distribution of per- and
polyfluoroalkyl substances (PFASs). Sci. Total
Environ., 2018. 636:105869. doi: 10.1016/j.
scitotenv.2018.04.380.
[6] Stubleski, J., Salihovic, S., Lind, L., Lind,
P.M., van Bavel, B., Kärrman, A. Changes in
serum levels of perfluoroalkyl substances
during a 10-year follow-up period in a large
population-based cohort. Environ. Int., 2016.
95:8692. doi: 10.1016/j.envint.2016. 08.002
[7] Hong, S., Khim, J.S., Park, J., Kim, M., Kim,
W.K., Jung, J., et al. In situ fate and
partitioning of waterborne perfluoroalkyl acids
(PFAAs) in the Youngsan and Nakdong River
Estuaries of South Korea. Sci. Total Environ.,
2013. 445-446:13645. doi: 10.1016/
j.scitotenv.2012.12. 040
[8] Coperchini, F., Awwad, O., Rotondi, M.,
Santini, F., Imbriani, M., Chiovato, L. Thyroid
disruption by perfluorooctane sulfonate
(PFOS) and perfluorooctanoate (PFOA). J.
Endocrinol. Invest., 2017. 40:105-121.
[9] Duntas, L., Keck, F.S., Rosenthal, J., Wolf,
C.H., Loos, U., Pfeiffer, E.F., Wochenshrift, K.
Single-compartment model analysis of
thyrotropin-releasing hormone kinetics in
hyper- and hypothyroid patients. Klin.
Wochenschr., 1990. 68:1013-1019.
[10] Barrett, K., Brooks, H., Boitano, H., Barman,
S. The Thyroid Gland. Review of Medical
Physiology, 23rd ed., The McGraw-Hill
Companies, Inc., 2010. 305-307.
[11] Jonklaas, J., Burman, K.D., Wang, H., Latham,
K.R. Single dose T3 administration: kinetics
and effects on biochemical and physiologic
parameters. Ther. Drug. Monit., 2015. 37(1):
110-118.
[12] Muratori, S., Rinaldi, S. Low-and high-
frequency oscillations in three dimensional
food chain systems. Siam J. Appl. Math., 1992.
52:1688-1706.
[13] Manorod, S., Rattanakul, C. Modelling the
population dynamics of brown planthopper,
Cyrtorhinus Lividipennis and Lycosa
Pseudoannulata. Advances in Difference
Equations, 2019. 2019:265. doi: 10.1186/
s13662-019-2217-y
[14] Rattanakul, C., Lenbury, Y., Krishnamara, N.,
Wollkind, D.J. Mathematical modelling of
bone formation and resorption mediated by
parathyroid hormone: Responses to
estrogen/PTH therapy. BioSystems., 2003. 70:
55-72.
[15] Lenbury, Y., Ruktamatakul, S.,
Amornsamankul, S. Modeling insulin kinetics:
responses to a single oral gluclose
administration or ambulatory-fed conditions.
BioSystems, 2001. 59:15-25.
[16] Rattanakul, C., Rattanamongkonkul, S. Effect
of Calcitonin on Bone Formation and
Resorption: Mathematical Modeling Approach.
Int. J. Mathl. Mod. Meth. Appl. Sci., 2011.
5(8):1363-1371.
[17] Villares, S., Knoepfelmacher, M., Salgado, L.
et al. Pulsatile Release and Circadian Rhythms
of Thyrotropin and Prolactin in Children with
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
277
Volume 18, 2023
Growth Hormone Deficiency. Pediatr. Res.,
1996. 39:10061011.
[18] Roelfsema, F., Veldhuis, J.D. Thyrotropin
Secretion Patterns in Health and Disease.
Endocrine Reviews, 2013. 34(5):619657.
[19] Brabant, G., Prank, K., Ranft, U.,
Schuermeyer, T.H., Wagner, T.O.F., Hauser,
H., Kummer, B., Feistner, H., Hesch, R.D.,
MüHlen, A.V.Z. Physiological Regulation of
Circadian and Pulsatile Thyrotropin Secretion
in Normal Man and Woman. The Journal of
Clinical Endocrinology & Metabolism, 1990.
70(2):403409.
Appendix
I = K = 0
y0
x
z
I = J = 0
J = 0
J = 0
I = 0
2
S
1
S
1
x
s
y
m
y
O
P
Q
RT
U
V
2
S
1
S
2
x
2
z
Fig. 1: The positions of the manifolds
0 , 0IJ
and
0K
where all inequalities stated in
Theorem 1 are satisfied. A solution trajectory tends toward a limit cycle where the fast, intermediate and slow
transitions are identified by the three, two and one arrows, respectively.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
278
Volume 18, 2023
Fig. 2: The positions of the manifolds
0 , 0IJ
and
0K
where the inequality stated in Theorem 2
are satisfied. A solution trajectory tends toward a stable equilibrium point
1
S
where the fast, intermediate and
slow transitions are identified by the three, two and one arrows, respectively.
y0
x
z
J = 0
J = 0
Q
I = 0
O
P
I = K = 0
1
S
1
x
4
x
m
x
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
279
Volume 18, 2023
I = K = 0
y0
x
z
I = J = 0
J = 0
J = 0
I = 0
2
S
1
S
1
x
2
x
s
y
m
y
O
P
Q
RT
N
2
S
1
S
Fig. 3: The positions of the manifolds
0 , 0IJ
and
0K
where all inequalities stated in
Theorem 3 are satisfied. A solution trajectory tends toward a stable equilibrium point
2
S
where the fast,
intermediate and slow transitions are identified by the three, two and one arrows, respectively.
Fig. 4: A simulation result of the system of equations (4)-(6) where
1 2 3
0.15, 0.30, 0.55,c c c
4 5 6 7 8 1 2 3 4 5 6
0.80, 0.70, 0.08, 0.21, 0.7, 0.95, 0.75, 0.90, 0.50, 0.30, 0.70,c c c c c k k k k k k
1 2 3
0.55, 0.02, 0.03, 0.95, 0.035,d d d

0 0.1, 0 0.5,xy
and
02z
.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
280
Volume 18, 2023
Fig. 5: A simulation result of the system of equations (4)-(6) where
1 2 3
0.015, 0.3, 0.55,c c c
4 5 6 7 8 1 2 3 4 5 6 1
0.8, 0.7, 0.08, 0.21, 0.7, 0.95, 0.75, 0.9, 0.5, 0.3, 0.7, 0.55,c c c c c k k k k k k d
23
0.02, 0.03, 0.5, 0.35, 0 0.1,d d x

0 0.1,y
and
05z
.
Fig. 6: A simulation result of the system of equations (4)-(6) where
1 2 3
0.05, 0.30, 0.55,c c c
4 5 6 7 8 1 2 3 4 5 6 1
0.80, 0.4, 0.1, 0.21, 0.7, 0.95, 0.75, 0.10, 0.5, 0.3, 0.7, 0.55,c c c c c k k k k k k d
23
0.11, 0.03, 0.15, 0.75, 0 0.5,d d x

0 0.5,y
and
01z
.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
281
Volume 18, 2023
Fig. 7: Comparison of the simulation results of the system of equations (4)-(6) for the case
60c
(without the
effect of PFOS/PFOA) and
60.08c
(with the effect of PFOS/PFOA). Here,
12
0.15, 0.30,cc
3 4 5 7 8 1 2 3 4 5 6 1
0.55, 0.80, 0.7, 0.21, 0.7, 0.95, 0.75, 0.9, 0.5, 0.3, 0.7, 0.55,c c c c c k k k k k k d
23
0.1, 0.03, 0.95, 0.035, 0 0.1, 0 0.5,d d x y

and
0 2.z
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The author contributed to the present research at all
stages including the formulation of the problem, the
theoretical analysis of the model, numerical
investigations as well as preparation of the
manuscripts.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
This work is supported by the Centre of Excellence
in Mathematics, Thailand and Mahidol University,
Thailand.
Conflict of Interest
The author has 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 SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.27
Chontita Rattanakul
E-ISSN: 2224-2856
282
Volume 18, 2023