Bone Remodeling Process and Covid-19: A Modelling Approach
CHONTITA RATTANAKUL
Department of Mathematics,
Faculty of Science,
Mahidol University,
Bangkok,
THAILAND
also with
Centre of Excellence in Mathematics,
MHESI,
Bangkok,
THAILAND
Abstract: - Apart from the effects on the lungs, COVID-19 caused by the severe acute respiratory syndrome
coronavirus 2 (SARS-CoV-2) also has effects on bone metabolism including the bone remodelling process. The
bone remodelling process involves bone formation by bone-forming cells (osteoblasts) and bone resorption by
bone-resorbing cells (osteoclasts). The infection with SARS-Cov-2 decreases the inhibiting effects of
Angiotensin-converting enzyme 2 (ACE2) on osteoclastic reproduction and inhibits the osteogenic ability of
osteoblasts which might lead to the imbalance in the bone remodeling process. In this study, we modify the
system of differential equations to investigate the effects of Covid-19 on bone formation and bone resorption.
The geometric singular perturbation method is utilized to analyze the modified model theoretically. To
illustrate the theoretical results, numerical investigations are also demonstrated. The results indicate that the
oscillations in the numbers of osteoclastic cells and osteoblastic cells observed in the clinical evidence could
still be expected when the effects of SARS-Cov-2 are incorporated, however, the oscillations occur at the
higher level of the number of osteoclasts and hence more bone loss might occur when infected with Covid-19.
Key-Words: - Covid-19, osteoclast, osteoblast, mathematical model, geometric singular perturbation
Received: December 17, 2022. Revised: August 24, 2023. Accepted: September 29, 2023. Published: October 17, 2023.
1 Introduction
The bone remodelling process aims to shape and
sculpt the skeleton during growth, to repair micro-
damaged bones that occur from everyday stress and
also to regulate calcium homeostasis, [1], [2], [3]. It
is a significant life-long process, the skeleton is
replaced almost 100% in the first year of life while
the bone remodelling process of the skeleton occurs
approximately 10% per year in adults, [4], [5], [6].
The process involves two types of bone cells, bone
tissue will be removed by bone resorbing cells,
osteoclasts (OCs), and new bone cells will be
formed by bone-forming cells, osteoblasts (OBs),
[7]. If the net bone resorption is over the net bone
formation after the completion of a bone
remodelling cycle, an imbalance will then occur
leading to the increase of bone loss and hence,
osteoporosis might be expected, [4], [5], [6].
COVID-19 is a contagious disease caused by
the SARS-CoV-2. The interaction of SARS-CoV-2
and ACE2 receptors expressed on cellular targets
such as alveolar cells of lungs and bone cells
enhances the production of inflammatory cytokines,
[8]. Also, the binding of SARS-CoV-2 and the
ACE2 receptors on osteoclastic cells induces the
differentiation of osteoclasts, [8]. However, it has
been observed that SARS-CoV-2 also enhances the
levels of angiotensin II by downregulating the
expression of ACE2 and then induces inflammatory
conditions, [9]. Moreover, angiotensin II enhances
the differentiation of OCs by enhancing the
expression of RANKL on OBs and hence,
accelerates the development of osteoporosis in the
rat model, [10]. On the other hand, the binding of
inflammatory cytokines and cytokine receptors
inhibits the osteogenic ability of osteoblasts and
hence, bone loss might be increased, [8].
In, [11], the mathematical model developed to
investigate the bone remodelling process is as
follows.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
316
Volume 18, 2023
1
1
1
r
dP sP
dt h C
23
2
2
2
r r P
dC BC s C
dt h P
5
43
3
r PB
dB r P s B
dt h P
where P is the PTH concentration above the basal
level, C is the number of active osteoclastic cells
and B is the number of active osteoblastic cells, [1].
Many researchers, [12], [13], [14], [15], [16],
[17], [18], [19], also investigated the effects of
involving hormones such as estrogen, parathyroid
hormone, vitamin D and calcitonin using
mathematical modelling. Mathematical models in
various forms were developed such as a system of
differential equations, a system of delay-differential
equations and an impulsive system of differential
equations. However, the effects of the inflection
with SARS-Cov-2 on the bone remodelling process
have never been taken into account. Therefore, the
mathematical model developed by, [11], is modified
in this paper to investigate bone resorption and bone
formation when the effects of the infection with
SARS-Cov-2 are also incorporated and the modified
model is then analyzed theoretically so that we
obtain the conditions for which different dynamics
behaviour can be occurred.
2 A Modified Mathematical Model
and Geometric Singular Perturbation
Analysis
In what follows, X stands for PTH concentration
above the basal level at time T, Y stands for the
number of active osteoclasts at time T and Z stands
for the number of active osteoblasts at time T. The
system of differential equations developed by, [11],
is modified here to investigate the effects of SARS-
Cov-2 on bone resorption and bone formation as
follows.
1
1
1
a
dX bX
dT k Y
(1)
23 4
2
2
23
Ya a X aY
dYZ b Y
dT k X k Y
(2)
6
5 7 3
4
a XZ
dZ a X a Z b Z
dT k X
(3)
Here, all parametric values in the system of
equations (1)-(3) are assumed to be positive.
As described in, [11], the change in PTH
concentration in blood above the basal level is
represented by equation (1). PTH is assumed to be
secreted from the parathyroid gland with the rate
presented by the first term on the right-hand side of
equation (1) which decreases when the number of
active osteoclasts increases so that the calcium
levels in blood is maintained to be within the normal
range, [11]. PTH is assumed to be removed from the
system as represented by the last term of equation
(1), [11].
In equation (2), as described in [11], It has been
observed that the reproduction of osteoclasts is
stimulated by PTH, [20], [21], [22]. However, it has
also been observed in, [19], as well that osteoclastic
reproduction is inhibited if the level of PTH
increases further, the saturation expression
23
2
2
a a X
kX
is then assumed for the stimulating effect
of PTH on the osteoclastic reproduction as
presented in the first term of the right-hand side of
equation (2), [11].
When infected with COVID-19, SARS-Cov-2
will bind with ACE2 receptors on osteoclastic cells
and induce the differentiation of osteoclasts, [8].
Moreover, the number of ACE2 receptors available
for binding with ACE2 is then decreased resulting in
the decrease of the inhibiting effect of ACE2 on
osteoclastic reproduction, [23], [24]. Therefore, the
number of osteoclastic cells will be increased.
Moreover, SARS-CoV-2 also enhances the levels of
angiotensin II by downregulating the expression of
ACE2, [9]. Angiotensin II then enhances the
differentiation of OCs by enhancing the RANKL
expression on OBs, [10]. Therefore, the term
4
3
aY
kY
is assumed on the first term of the right-hand side of
equation (2) to account for the effect of COVID-19
on the reproduction of osteoclastic cells. Noted that
the effects of both PTH and the infection with
COVID-19 on osteoclastic reproduction require the
osteoclast differentiation factor (ODF) and its
receptor on osteoclastic cells, [19], the term YZ is
then also presented. Active osteoclastic cells are
assumed to be removed from the system as
presented by the last term on the right-hand side of
equation (2), [11].
In equation (3), as described in [11], PTH has
been observed to stimulate osteoblastic reproduction
and the first term on the right-hand side is then
assumed to account for PTH stimulating effect on
osteoblastic reproduction, [11]. Moreover, PTH has
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
317
Volume 18, 2023
also been observed to inhibit osteoblastic
reproduction and the second term on the right-hand
side is then assumed to account for PTH inhibiting
effect on osteoblastic reproduction, [11]. On the
other hand, the interaction of SARS-CoV-2 and
ACE2 receptors enhances the production of
inflammatory cytokines, [8]. Since the binding of
inflammatory cytokines and cytokine receptors
inhibits the osteogenic ability of osteoblasts, [8], the
number of osteoblastic cells will then decrease when
infected with COVID-19 and hence the third term
on the right-hand side is presented to account for the
effect of Covid-19 on the reproduction of
osteoblastic cells. Active osteoblastic cells are
assumed to be removed from the system by the last
term on the right-hand side, [11].
The changes in the levels of PTH in the blood
occur within minutes to control the calcium levels in
the blood while the bone resorption process takes
approximately 2 weeks and the bone formation
process takes approximately 3 months, [25], [26].
Hence, we assume that the dynamics of PTH, OCs
and OBs are fastest, intermediate and slowest,
respectively. We then analyze our modified model
by utilizing the geometric singular perturbation
technique. Firstly, we scale the variables of the
system by two small dimensionless positive
parameters
and
. The new parameters are
introduces as follows.
Letting
* * *
12
12
2
* * * * * * *
, , , , , ,
T a T Z a
X Y Z T
x y z t c c
X Y Z T X Y X
* * * * *
**
3 5 6
4
3 4 5 6
**
, , , ,
T Z a T X a T a
T Z a
c c c c
XZ
*
*
73
12
7 1 2 3 1 1
2
**
*
, , , , ,
T a k
kk
c m m m d T b
YY
X
*
2
2,
Tb
d
and
*
3
3
Tb
d
, our model then becomes
1
1
1
,,
c
dx d x f x y z
dt m y
(4)
23 4
2
2
3
2
,,
c c x cy
dy yz d y
dt m y
mx
g x y z
(5)
6
5 7 3
4
,,
c xz
dz c x c z d z h x y z
dt m x
(6)
This means that, PTH possesses a faster time
response than osteoclasts, while, osteoblasts have
slowest dynamics provided that
and
are small
which is consistent with the clinical evidence, [19],
[25], [26].
For our system of equations (4)-(6), we now
show that the manifolds
0f
,
0g
and
0h
are shaped and located as shown in Figure 1
(Appendix) and Figure 2 (Appendix) with some
appropriate parametric values.
Manifold
0f
By setting
0
dx
dt
in equation (4), we obtain the
manifold
0f
which is represented by
1
1
1
0
cdx
my
or
1
11
c
x r y
d m y
Here,
ry
is independent of
z
. As a result,
ry
is also parallel to the z-axis. Moreover, the
intersection of
ry
and the (x,z)-plane occurs
along the straight line
1
1
11
.
c
xx
dm
In addition, we can see that
ry
is a decreasing
function and
1
11
lim lim 0
yy
c
ry d m y
that is,
0x
as
y
along the surface
.x r y
Manifold
0g
By setting
0
dy
dt
in equation (5), we obtain the
manifold
0g
which is represented by
23 4
2
2
23
0
c c x cy yz d y
m x m y
or
23 4
2
2
23
0
c c x cy
y z d
m x m y
composing of two parts, the trivial manifold
0y
,
and the nontrivial manifold
2
23
2
2 3 3 4 2
,
m x m y
z s x y
c c x m y c y m x
We can see that on the
,xz
-plane,
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
318
Volume 18, 2023
2
2
23
.
mx
z u x
c c x
Note that
ux
intersects the
z
-axis at the point in
which
0x
and
2
1
2
.
m
zz
c
Since,
2
3 2 3 2
2
23
2c x c x c m
ux c c x
,
0ux
is then
occurred in the first octant at the point in which
2
22
22
33
.
cc
x m x
cc
We remark that
2
x
is always positive, when all
parameters are assumed to be positive.
Moreover,
2
3 3 2 3 2
32
23
2 3 2 3
22
22 c c x c x c m
c x c
ux c c x c c x
and
20x
, then
2
2 3 2 3 2 2
24
2 3 2
22
0.
c c x c x c
ux c c x
Therefore,
z u x
attains its relative minimum at
the point at which
2
xx
,
2
22
22
2 3 2
.
mx
z u x z
c c x
On the
,yz
-plane,
2 3 2
2 3 2 4 2
.
m m m y
z v y
c m c c m y
We can see that,
vy
intersects the
y
-axis at the
point for which
0z
and
3.ym
Note that
2
*
2 4 2
lim
y
m
v y z
c c m
. The intersection of
vy
and the
z
-axis occurs at the point for which
0y
and
2
1
2
.
m
zz
c
Since
2
4 2 3
2
2 3 2 4 2
0
c m m
vy c m c c m y
then
z v y
is an decreasing function when all
parameters are assumed to be positive.
On the
,xy
-plane,
3.ym
Since
2
23
2
2 3 3 4 2
,
m x m y
z s x y
c c x m y c y m x
and
22
3 3 2 3 2
2
2
2 3 3 4 2
2m y c x c x c m
S
xc c x m y c y m x
then
0
S
x
at
2
xx
.
Note that
2
2
4 3 2
2
2
2 3 3 4 2
0
c m m x
S
yc c x m y c y m x
when all parameters are assumed to be positive.
Manifold
0h
By setting
0
dz
dt
in equation (6), we obtain the
manifold
0h
which is represented by
6
5 7 3
4
0
c xz
c x c z d z
mx
or
54
6 7 3 4
.
c x m x
z w x
c x c d m x
Here,
wx
is parallel to the
y
-axis because
wx
is independent of the variable
y
. Also,
wx
intersects the
,yz
-plane along the
y
-axis.
Since,
2
5 6 7 3
2
6 7 3 4 7 3
2
5 4 7 3 5 4 7 3
2
6 7 3 4 7 3
2
c c c d x
wx c c d x m c d
c m c d x c m c d
c c d x m c d
then
0wx
if
22
6 7 3 4 7 3 4 7 3
2 0.c c d x m c d x m c d
That is
4 7 3
6 7 3
2
22
4 7 3 4 7 3 6 7 3
6 7 3
0
m c d
xc c d
m c d m c d c c d
c c d
Thus,
wx
has no relative minimum or maximum
where all parametric values are assumed to be
positive. Moreover,
0wx
for
0.x
Thus,
wx
is an increasing function in the first octant.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
319
Volume 18, 2023
Theorem 1 Given that
and
are sufficiently
small. A limit cycle exists for the system of
equations (4)-(6) if the inequalities
2 3 1,x x x
(7)
34
,yy
(8)
and
2 3 4 1
z z z z
(9)
hold.
Note that Theorem 1 can be proven by utilizing the
geometric singular perturbation techniques, [27],
[28].
In Figure 1 (Appendix), the curve
0gh
intersects the curve
0fg
at the point
1
S
on
the
,xz
-plane and the point
2
S
is located between
the points
2 2 2
,,x y z
and
11
,0,xz
along the curve
0fg
. Note that the high, intermediate and
slow speed transitions are represented by the three,
two and one arrows, respectively.
Starting from the generic point A =
0 0 0
,,x y z
for which
0 0 0
, , 0f x y z
in Figure 1 (Appendix).
Here, the position of point A is assumed to be as
shown in Figure 1 (Appendix). A transition in the
direction of decreasing x at a fast speed will bring
the solution trajectory from point A to point B
located on the manifold
0f
, since
0f
here.
Since
0g
here, a transition of the solution
trajectory with an intermediate speed in the
direction of decreasing
y
will bring the system to
the curve
0fg
where a point C is reached. A
transition of the solution trajectory in the direction
of increasing z along this curve with a slow speed,
since
0h
here, will then reach some point D on
this curve. The stability of the submanifold will be
lost and an intermediate transition of the solution
trajectory will then bring the system to the other
stable part of
0fg
where a point E is
reached. A transition of the solution trajectory at
slow speed in the direction of decreasing
z
, since
0h
here, will bring the system to the point F in
which the stability of the submanifold is lost again.
An intermediate transition will then bring the
solution trajectory to the other stable part of
0fg
where point G is reached. A transition
in the direction of increasing
z
with a slow speed
will then bring the solution trajectory to point D
again, since
0h
here, resulting in the closed orbit
DEFG and hence, we obtain a limit cycle for our
model of equations (4)-(6).
Theorem 2 Given that
and
are sufficiently
small, if the inequalities
3 2 1,x x x
(10)
43
,yy
(11)
and
5 3 2 4 1
z z z z z
(12)
hold then the equilibrium point
2
S
of the system of
equations (4)-(6) is stable.
Starting from the generic point A =
0 0 0
,,x y z
for which
0 0 0
, , 0f x y z
in Figure 2 (Appendix).
The position of A is assumed to be located as in
Figure 2 (Appendix). Since
0f
here, a transition
of the solution trajectory in the direction of
decreasing x will bring the system to the manifold
0f
with a fast speed where a point B is
reached. An intermediate transition of the solution
trajectory in the direction of decreasing
y
on the
curve
0fg
will bring the system to a point
C, since
0g
here. A transition of the solution
trajectory along this curve at slow speed in the
direction of increasing z, since
0h
here, will then
bring the system to some point D. The stability of
the submanifold will be lost and an intermediate
transition of the solution trajectory will then bring
the system to the other stable part of
0fg
where a point E is reached. A transition of the
solution trajectory at a slow speed in the direction of
decreasing
z
will bring the system to the
equilibrium point
2
S
where
0f g h
. Thus, the
equilibrium point
2
S
is stable.
3 Numerical Investigations
In this section, we demonstrate the theoretical
predictions by carrying out numerical investigations
for each case. All computer simulations are
generated by MATLAB using the Runge-Kutta 4th-
order method.
In Figure 3 (Appendix), numerical simulation is
presented to demonstrate the theoretical result stated
in Theorem 1 for which a limit cycle is expected.
The parametric values in the system of equations
(4)-(6) are chosen to satisfy all inequalities stated in
Theorem 1 as follows:
12
0.05, 0.009,cc
3 4 5 6 7
0.675, 0.1, 0.01, 0.005, 0.001,c c c c c
1 2 3 4 1 2
0.1, 0.5, 1, 0.025, 0.1, 0.3,m m m m d d
30.01, 0.1, 0.9, 0 0.5, 0 0.1,d x y
and
0 0.1.z
In Figure 3 (Appendix), we can see that
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
320
Volume 18, 2023
a limit cycle exists for the system of equations (4)-
(6) and a periodic solution occurs as theoretically
predicted in Theorem 1.
In Figure 4 (Appendix), a numerical simulation is
presented to demonstrate the theoretical result
indicated in Theorem 2. The parametric values in
the system of equations (4)-(6) are chosen to satisfy
all inequalities stated in Theorem 2 as follows:
1 2 3 4 5
0.05, 0.0009, 0.882, 0.1, 0.01,c c c c c
6 7 1 2 3
0.005, 0.0001, 0.1, 0.5, 1,c c m m m
4 1 2 3
0.025, 0.1, 0.3, 0.01, 0.01, 0.5,m d d d
0 0.1,x
0 0.5,y
and
0 0.5.z
We can see
that the equilibrium point
2
S
is stable as
theoretically predicted in Theorem 2.
4 Discussion
In addition to the numerical simulations represented
in Figure 3 (Appendix) and Figure 4 (Appendix), let
us consider an example of the computer simulations
in Figure 5 (Appendix) for which the number of
osteoclasts when there are no effects from SARS-
Cov-2 infection (
47
0cc
) is compared to the
number of osteoclasts when there is the effects from
SARS-Cov-2 infection (
47
0.1, 0.001cc
) when
the other parameters are the same as the values used
in Figure 3 (Appendix). Also, let us consider an
example of the computer simulations in Figure 6
(Appendix) for which the number of osteoclasts
when there is no effects from SARS-Cov-2 infection
(
47
0cc
) is comparing to the number of
osteoclasts when there is the effects from SARS-
Cov-2 infection (
47
0.1, 0.0001cc
) when the
other parameters are the same as the values used in
Figure 4 (Appendix). We can see that the number of
osteoclasts in the cases when there is the effect of
SARS-Cov-2 infection is higher than that of the
cases when there is no effect of SARS-Cov-2
infection in both Figure 5 (Appendix) and Figure 6
(Appendix) whereas the number of osteoblasts in
the cases when there is the effect of SARS-Cov-2
infection is lower than that of the cases when there
is no effect of SARS-Cov-2 infection in Figure 6
(Appendix). Therefore, more bone loss might be
expected when infected with Covid-19.
5 Conclusion
The model developed by, [1], is modified to
incorporate the effects of Covid-19 on the
reproductions of osteoblastic cells and osteoclastic
cells. The modified model is then analyzed
theoretically using the geometric singular
perturbation technique so that we obtain the
conditions for which different behaviours of
solution can be expected. Numerical simulations are
also provided to demonstrate the theoretical results.
The results indicate that a periodic behaviour which
has been observed clinically in the oscillations of
the number of osteoclasts and osteoblasts, [19], as
well as the pulsatile secretions of PTH, [29], could
also be exhibited by our model. In addition, the
simulation results in Figure 5 (Appendix) and
Figure 6 (Appendix) indicate that the number of
osteoclasts in the cases where there is the effect of
SARS-Cov-2 infection is higher than that of the
cases when there is no effect of SARS-Cov-2
infection whereas the number of osteoblasts in the
cases when there is the effect of SARS-Cov-2
infection are lower than that of the cases when there
is no effect of SARS-Cov-2 infection and hence,
more bone loss can be expected when infected with
Covid-19. The modified model in this paper might
be modified further to investigate treatments for
osteoporosis patients when infected with COVID-
19.
Acknowledgement:
This research has received funding support from the
NSRF via the Program Management Unit for
Human Resource and institutional Development,
Research and Innovation (grant number
B05F640231). Also, this research is partially
supported by Mahidol University, Thailand.
References:
[1] Russell, T., Turner, B., Lawrence, R.,
Thomas, C.S. Skeletal effects of estrogen.
Endocr. Rev., 15(3) (1994): p.275-300.
[2] Lobo, R.A., Kelsey, J.L. and Marcus, R.
Menopause: Biology and Pathobiology.
Academic Press, 2000, pp. 287-307.
[3] Albright, J.A. and Sauders, M. The Scientific
Basis of Orthopaedics. Norwalk, Conn.:
Appleton & Lange, 1990.
[4] Lewellen, T.K., Neip, W.B., Murano, R.,
Hinn, G.M. and Chesnut, C.H. Absolute
Measurement of Total-Body Calcium by the
Ar-37 Method- Preliminary Results: Concise
Communication. J. Nucl. Med., 18 (1977):
p.929-932.
[5] Morley, P., Whitfield, J.F. and Willick, G.E.
Parathyroid hormone: an anabolic treatment
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
321
Volume 18, 2023
for osteoporosis. Curr Phar Design, 7 (2001):
p.671-687.
[6] Parfitt, A.M. Targeted and nontargeted bone
remodelling: relationship to basic
multicellular unit origination and progression.
Bone, 30 (2002): p.5-7.
[7] Heersche, J.N.M. and Cherk, S. Metabolic
Bone Disease: Cellular and tissue
mechanisms. Boca Raton, FI: CRC Press,
1989.
[8] Sapra, L., Saini, C., Garg, B., Gupta, R.,
Verma, B., Mishra, P.K., Srivastava, R.K.
Long-term implications of COVID-19 on
bone health: pathophysiology and
therapeutics. Inflammation Research,
71(2022): p.1025-1040.
[9] Kuba, K., Imai, Y., Penninger, J. Angiotensin-
converting enzyme 2 in lung diseases. Curr
Opin Pharmacol. 6 (2006): 271–6.
[10] Shimizu, H., Nakagami, H., Osako, M.K.,
Hanayama, R., Kunugiza, Y., Kizawa, T.,
Tomita, T., Yoshikawa, H., Ogihara, T.,
Morishita, R. Angiotensin II accelerates
osteoporosis by activating osteoclasts. FASEB
J. 22(7) (2008): 2465-75.
[11] Rattanakul, C., Lenbury, Y., Krishnamara, N.
and Wollkind, D.J. Mathematical Modelling
of Bone Formation and Resorption Mediated
by Parathyroid Hormone: Responses to
Estrogen/PTH Therapy. BioSystems,
70(2003): p.55-72.
[12] Chaiya I., Rattanakul C., Rattanamongkonkul
S., Kunpasuruang W., Ruktamatakul S.
Effects of Parathyroid Hormone and
Calcitonin on Bone Formation and
Resorption: Mathematical Modeling
Approach. Int. J. Math. Comp. Simul. 5(6)
(2011): p.510-519.
[13] Ayati, B.P., Edwards, C.M., Webb, G.F. et
al. A mathematical model of bone remodeling
dynamics for normal bone cell populations
and myeloma bone disease. Biol Direct. 5(28)
(2010).
[14] Rattanakul C., Rattanamongkonkul S. Effect
of Calcitonin on Bone Formation and
Resorption: Mathematical Modeling
Approach. Int. J. Mathl. Mod. Meth. Appl. Sci.
5(8) (2011): p.1363-1371.
[15] Chen-Charpentiera B.M., Diakite I. A
mathematical model of bone remodeling with
delays. Journal of Computational and Applied
Mathematics 291(2016):76–84.
[16] Chaiya, I., Rattanakul, C. Effects of Prolactin
on Bone Remodeling Process with
Parathyroid Hormone Supplement: An
Impulsive Mathematical Model. Advances in
Difference Equations (2017) 2017: 147.
[17] Trejo I., Kojouharov H. A Mathematical
Model to Study the Fundamental Functions of
Phagocytes and Inflammatory Cytokines
During the Bone Fracture Healing
Process. Letters in Biomathematics 7(1)
(2020), p.171–189.
[18] Aekthong, S., Rattanakul, C. Dynamical
Modelling of Bone Formation and Resorption
under Impulsive Estrogen Supplement:
Effects of Parathyroid Hormone and Prolactin.
International Journal of Mathematics and
Mathematical Sciences. Article ID 5435876,
(2021) 17 pages. https://doi.org/10.1155/2021/
5435876.
[19] Kroll, M.H. Parathyroid hormone temporal
effects on bone formation and resorption.
Bull. Math. Bio, 62 (2000): p.163-188.
[20] Dempster, D.W., Cosman, F., Parisisen, M.,
Shen, V. and Linsay, R. Anabolic actions of
parathyroid hormone on bone. Endocr Rev. 14
(1993): p.690-709.
[21] Weryha, G. and Leclere, J. Paracrine
regulation of bone remodeling. Horm Res. 43
(1995): p.69-75.
[22] Brown, E.M. Extracellular Ca2+ sensing,
regulation of parathyroid cell function, and
role of Ca2+ and other ions as extracellular
(first) messengers. Physiol Rev. 71 (1991):
p.371-411.
[23] Zheng, K., Zhang, W., Xu, Y., Geng, D.
Covid-19 and the bone: underestimated to
consider. European Review for Medical and
Pharmacological Sciences. 24 (2020):
p.10316-10318.
[24] Tao, H., Bai, J., Zhang, W., Zheng, K., Guan,
P. Ge, G., Li, M., Geng, D. Bone biology and
COVID-19 infection: Is ACE2 a potential
influence factor? Medical Hypotheses. 144
(2020): 110178.
[25] Whitfield, J.F., Morley, P. and Willick, G.E.
The parathyroid hormone: an unexpected
bone builder for treating osteoporosis. Austin,
Tex: Landes Bioscience Company, 1998.
[26] The epidemiology and pathogenesis of
osteoporosis [Online 2002]. Available from
http://www.endotext.org/parathyroid/parathyr
oid11/parathyroid11.htm.
[27] Kaper, T.J. An introduction to geometric
methods and dynamical systems theory for
singular perturbation problems. Analyzing
multiscale phenomena using singular
perturbation methods, Proc. Symposia Appl
Math, 56 (1999).
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
322
Volume 18, 2023
[28] Rinaldi, S. and Muratori, S. A separation
condition for the existence of limit cycle in
slow-fast systems. Appl Math Modelling. 15
(1991): p.312-318.
[29] Prank, K., Harms, H., Brabant, G. and Hesch,
R.D. Nonlinear dynamics in pulsatile
secretion of parathyroid hormone in normal
human subjects. Chaos. 5 (1995): p.76-81.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
323
Volume 18, 2023
Appendix
1
S
2
x
0
2
z
1
z
2
S
z
y
0fg
0h
0g
0f
4
z
1
x
0gh
3
x
3
y
3
z
4
y
0fg
x
A
B
C
DE
F
G
Fig. 1: The manifolds
0 , 0fg
and
0h
where all conditions in Theorem 1 are satisfied. Here, the
fast, intermediate, and slow transitions of the trajectories are represented by three arrows, two arrows, and one
arrow, respectively. The solution trajectory tends towards a limit cycle in this case.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
324
Volume 18, 2023
1
S
2
x
0
2
z
1
z
2
S
z
y
0fg
0h
0g
0f
4
z
1
x
0gh
3
x
3
y
3
z
4
y
0fg
x
5
z
A
B
C
DE
Fig. 2: The manifolds
0 , 0fg
and
0h
where all conditions in Theorem 2 are satisfied. Here, the
fast, intermediate, and slow transitions of the trajectories are represented by three arrows, two arrows, and one
arrow, respectively. The solution trajectory tends towards the stable equilibrium point
2
S
in this case.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
325
Volume 18, 2023
Fig. 3: A numerical simulation of the model system (4)-(6), with the parametric values satisfying all
inequalities identified in Theorem 1. Here,
1 2 3 4 5 6
0.05, 0.009, 0.675, 0.1, 0.01, 0.005,c c c c c c
7 1 2 3 4 1 2 3
0.001, 0.1, 0.5, 1, 0.025, 0.1, 0.3, 0.01, 0.1, 0.9, 0 0.5, 0 0.1,c m m m m d d d x y
and
0 0.1.z
A limit cycle exists as predicted in Theorem 1.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
326
Volume 18, 2023
Fig. 4: A numerical simulation of the model system (4)-(6), with parametric values satisfying all inequalities
stated in Theorem 2. Here,
1 2 3 4 5 6 7 1
0.05, 0.0009, 0.882, 0.1, 0.01, 0.005, 0.0001, 0.1,c c c c c c c m
2 3 4 1 2 3
0.5, 1, 0.025, 0.1, 0.3, 0.01, 0.01, 0.5, 0 0.1, 0 0.5,m m m d d d x y
and
0 0.5.z
The equilibrium point
2
S
is stable as predicted in Theorem 2.
Fig. 5: Comparison of the numerical simulations of the system of equations (4)-(6) for the case when
47
0cc
(without effects of SARS-Cov-2 infection) and
47
0.1, 0.001cc
(with effects of SARS-Cov-2
infection). Here,
1 2 3 5 6 1 2 3 4
0.05, 0.009, 0.675, 0.01, 0.005, 0.1, 0.5, 1, 0.025,c c c c c m m m m
1 2 3
0.1, 0.3, 0.01, 0.1, 0.9, 0 0.5, 0 0.1, 0 0.1.d d d x y z
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2023.18.32
Chontita Rattanakul
E-ISSN: 2224-2856
327
Volume 18, 2023
Fig. 6: Comparison of the numerical simulations of the system of equations (4)-(6) for the case when
47
0cc
(without effects of SARS-Cov-2 infection) and
47
0.1, 0.0001cc
(with effects of SARS-Cov-2
infection). Here,
1 2 3 5 6 1 2 3 4
0.05, 0.0009, 0.882, 0.01, 0.005, 0.1, 0.5, 1, 0.025,c c c c c m m m m
1 2 3
0.1, 0.3, 0.01,d d d
0.01, 0.5, 0 0.1, 0 0.5, 0 0.5.x y z
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The author contributed in 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 research has received funding support from the
NSRF via the Program Management Unit for
Human Resource and institutional Development,
Research and Innovation (grant number
B05F640231). Also, this research is partially
supported by 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.32
Chontita Rattanakul
E-ISSN: 2224-2856
328
Volume 18, 2023
