# Dynamical Modelling of Bone Formation and Resorption under Impulsive Estrogen Supplement: Effects of Parathyroid Hormone and Prolactin

Authors: Supassorn Aekthong; Chontita Rattanakul
Journal: International Journal of Mathematics and Mathematical Sciences
 (2021)

---

## Abstract

Osteoporosis, a bone metabolic disease, is one of the major diseases occurring in aging population especially in postmenopausal women. A system of impulsive differential equations is developed in this paper in order to investigate the effects of parathyroid hormone and prolactin on bone-forming cells, namely, osteoblasts, and bone-resorbing cells, namely, osteoclasts, under the impulsive estrogen supplement. The theoretical analysis of the developed model is carried out so that we obtain the conditions on the system parameters in which the stability and permanence of the model can occur. Computer simulations are also provided to illustrate the theoretical predictions.


---

## Body

## 1. Introduction

Osteoporosis is a disease which shows no symptoms until there is a bone fracture incidence [1]. It is characterized by an impairment of bone mass, measured by a Bone Mineral Density (BMD) test as the T-score is less than or equal to 2.5, and microarchitectural deterioration of bone tissue [2–5]. This disease causes a major public health issue for all races and both sexes [1]. Since BMD decreases as the population ages, osteoporosis occurs prevalently in menopause woman and elderly men between 75 and 80 years of age [6]. More than 200 million people worldwide were approximated to be osteoporotic patients with hip fractures [7]. In Europe and the United States, a report revealed that 30% of women are diagnosed with osteoporosis [8]. It was reported in Thailand that 19.8% of postmenopausal women are prevalently osteoporotic [9]. Furthermore, approximately 40% of women who experience menopause and 30% of men might have an osteoporosis-related fracture in the remaining lifetime [10].Bone is a living tissue [11]. There is a process in its place to resorb old bone and form new bone to maintain its strength and health, called bone-remodeling process, and such a process occurs throughout a person’s life [10, 12, 13]. In osteoporosis patients, the loss of bone mass occurs because the removal of older bone is more than the replacement of new bone which is an imbalance between bone resorption and bone formation [10]. There are two families of cells mainly involved in the process of bone remodeling, namely, osteoclastic cells, breaking down the bone cells, and osteoblastic cells, forming the bone tissue [12, 13]. Bone-remodeling process is composed of four phases sequentially, which are activation, resorption, reversal, and formation [12]. In this process, osteoblast lineage cells play a crucial role in the activation stage by taking action on blood cell precursors or hematopoietic cells to establish osteoclastic resorption of bones. Under a layer of lining cells, the removal of bone is done by osteoclasts in the resorption step. The mineral here is dissolved, and the bone matrix is broken down. In the reversal phase, the surface of the resorbed bone creates a thin layer of protein to make it ready for the following stage. Finally, in the formation step, osteoblasts acting in successive waves start to form up new bones by laying down multiple layers of matrix in an orderly arrangement. Then, some of these osteoblasts stay inside the bone and become osteocytes, remaining in contact to adjacent osteocytes and osteoblasts on the bone surface. This process is successfully finished in 3-4 months, while the first three stages take 2 to 3 weeks to finish [12]. There are many factors involved in bone-remodeling process such as parathyroid hormone (PTH), prolactin (PRL), calcitonin, and estrogen.PTH plays an important role in the bone-remodeling process [14]. It is released from the parathyroid gland, and both stimulating and inhibiting effects have been reported on the development of osteoblasts from its progenitors [14–18]. On the one hand, PTH indirectly enhances the development of osteoclasts from its progenitors by operating through osteoblasts since osteoclasts lack PTH receptors while preosteoblast precursors and preosteoblasts possess them [14–18]. On the other hand, PRL is a polypeptide hormone synthesized in and secreted from specialized cells of the anterior pituitary gland, the lactotrophs [19]. According to [20], PRL receptors have been reported on osteoblasts, and hence, PRL then has effects on bone-remodeling process as well. PRL enhances bone resorption by increasing receptor activator of NF-ligand (RANKL) and decreasing osteoprotegerin (OPG) expressions by osteoblastic cells [21].Bone mass reduces significantly due to the declining secretion of estrogen during and after the menopause years in women [22]. Estrogen deficiency in postmenopausal women leads to accelerating resorption of bone by osteoclasts, and osteoporosis will then occur [23]. The direct and indirect studies indicated that estrogen also plays a key role on skeletal health in men [24]. Several hormones play an important role in the process of bone remodeling and also in the treatments of osteoporotic patients. Estrogen replacement therapy has a beneficial effect on postmenopausal women in preventing the loss of bone as shown in clinical studies [25, 26]. In a long-term study of following up postmenopausal women who take estrogen replacement therapy, in vivo study results demonstrated that replacement estrogen therapy can prevent or slow down osteoporosis in postmenopausal patients compared to an unsupplemented control group [27]. A recent research using cancellous rats reported that decreased estrogen levels led to bone changes which affect flow of interstitial around osteocytes [28]. Many studies have shown that estrogen impacts positively in improving bone mineral density, and taking lower doses of estrogen is shown to be effective and cause fewer side effects. At the cellular level, estrogen increases the level of OPG which is able to bind with RANKL to block the differentiation, the activity, and the survival of osteoclasts.Since bone is an alive and dynamic tissue [11], many researchers have been trying to describe the process of bone remodeling inside our bodies using several types of mathematical models. Chudtong et al. [29] introduced an impulsive system to describe bone-remodeling process involving PTH supplements by extending Rattanakul et al.’s model in [30]. Rattanakul et al. [31] constructed a mathematical model as a system of nonlinear differential equations to investigate bone formation and resorption process based on the effect of calcitonin. Motivated by [31], Panitsupakamon and Rattanakul [32] modified the model proposed in [31] by incorporating the time delay observed in the bone-remodeling process. The influence of estrogen supplement is also considered by adding a term to the dynamics of the active osteoclast population in the modified model. Chaiya and Rattanakul [33] proposed an impulsive mathematical model of bone-remodeling process incorporating the effects of prolactin and the impulsive control strategies of parathyroid hormone supplement on osteoblasts and osteoclasts. The dosage of parathyroid hormone can be added appropriately to the system to ensure the suitable levels of active osteoblasts and active osteoclasts. They also proposed another model accounting for the effects of PTH and calcitonin on bone-remodeling process together with the effects of impulsive treatments of estrogen in [34]. However, a mathematical model of bone resorption and bone formation consisting of osteoblastic cells, osteoclastic cells, parathyroid hormone, and prolactin with the effect of impulsive treatment of estrogen has never been established.


## 2. Model Development

In this section, we propose the following impulsive differential equation model to investigate the dynamics of bone-forming cells and bone-resorbing cells based on the effects of PTH, PRL, and estrogen supplements:(1)dxdt=a1k1+z+a2y−b1x,dydt=a3k2+z+a4x−b2y,dzdt=a5xk3+x2+a6yk4+y2zw−b3z,dwdt=a7x−a8xwk5+x−a9ywk6+y−b4w,t≠nT,with(2)Δzt=−ρzt,Δwt=μ,t=nT,where Δzt=zt+−ztandΔwt=wt+−wt. In what follows, xt,yt,zt, and wt account for the level of PTH in blood at time t, the level of PRL in blood at time t, the number of active osteoclasts at time t, and the number of active osteoblasts at time t, respectively. T accounts for the period between each impulsive estrogen treatment, n∈Z+=1,2,3,…,ρ accounts for the inhibiting effect of estrogen treatment on osteoclasts, 0<ρ<1,andμ accounts for the stimulating effect on osteoblasts of estrogen treatment, μ>0. In addition, all parameters in (1) and (2) are assumed to be positive.PTH secretion from the parathyroid gland is controlled by the calcium level in blood. Once the number of active osteoclasts increases, blood calcium levels then increase due to the increase in bone resorption resulting in the decrease of PTH secretion in order to maintain calcium level within the normal range [35, 36]. Hence, we then assumed that the rate of change of the concentration of PTH in blood at time t denoted by dx/dt in equation (1) varies inversely with the number of active osteoclasts as shown in the first term on the right-hand side of the first equation in equation (1). According to [37], an increase in the level of plasma PRL also enhances the release of PTH from the parathyroid gland and hence the second term on the right-hand side denotes the stimulating effect of PRL on PTH. The last term stands for the removal rate of PTH from the system.However, PRL stimulates osteoclastic differentiation by decreasing OPG and increasing RANKL expressions by osteoblastic cells [21]. RANKL then binds to RANK, its receptor on preosteoclast, and stimulates the differentiation of osteoclasts. To balance the bone-remodeling cycle, when the number of osteoclasts increases, the secretion of PRL is then decreased. In the second equation of (1), the rate of change of the level of PRL in blood denoted by dy/dt is then assumed to vary inversely with the number of active osteoclasts as shown in the first term on the right-hand side. On the other hand, the increase in the number of osteoclasts leads to the increase in the calcium level in blood and then the secretion of PTH from the parathyroid gland, as well as the secretion of PRL from the anterior pituitary gland, which also enhances the secretion of PTH will be decreased in order to maintain the calcium level in blood within the normal range [35, 36]. Then, the second term on the right-hand side stands for the stimulating effect of PTH on the release of PRL from the anterior pituitary gland. The last term stands for the removal rate of PRL from the system.Osteoclasts lack PTH and PRL receptors, whereas osteoblasts possess them. Therefore, the effects of PTH and PRL on the differentiation of osteoclasts are both indirect effects. However, it has been reported in [17] that the differentiation of active osteoclasts also requires the production of osteoclast differentiation factor (ODF) and its receptor on osteoclasts as well and hence the differentiation of active osteoclasts requires the presence of both osteoblastic and osteoclastic cells. Even though PTH has the stimulating effects on the differentiation of active osteoclasts as reported in [18, 38], it has also been observed clinically in [17] that PTH inhibits the differentiation of active osteoclasts when the level of PTH increases further as well. On the other hand, PRL enhances bone resorption by decreasing OPG and increasing RANKL expressions by osteoblastic cells [21]. RANKL then binds to RANK, its receptor on preosteoclast, and stimulates the differentiation of osteoclasts [21]. Thus, the rate of change of the number of active osteoclasts denoted by dz/dt is then assumed to depend on the stimulating and inhibiting effects of PTH and PRL as presented on the first term on the right-hand side of the third equation in (1). The last term stands for the removal rate of active osteoclasts.Osteoblasts possess both PTH receptors and PRL receptors, and hence, the effects of PTH and PRL on the differentiation of osteoblasts are both direct effects. PTH prevents a suicidal process called apoptosis of osteoblasts [39, 40]. In addition, both stimulating and inhibiting effects of PTH on osteoblastic differentiation process have been clinically observed depending on the differentiation stages [14]. As for the stimulating effect of PTH on the differentiation of osteoblasts, we assumed that the rate of change of the number of active osteoblasts denoted by dw/dt varies directly with the level of PTH as presented on the first term on the right-hand side of the fourth equation of (1), whereas the inhibiting effect of PTH on the differentiation of osteoblasts is presented on the second term on the right-hand side. On the other hand, it has been reported in [41] that the osteoblast number (DNA content) is declined significantly due to PRL, resulting in a decrease in osteoblast proliferation, and hence, the third term on the right-hand side is assumed to represent the inhibiting effect of PRL on osteoblasts. The last term represents the removal rate of active osteoblasts.The first equation of (2) accounts for the inhibitory effect of impulsive treatment with estrogen supplement on the number of active osteoclasts. The second equation accounts for the simulative effect of impulsive treatment with estrogen supplement on the number of active osteoblasts.Several clinical observations indicate that the dynamics of PTH and PRL are so much faster than the dynamics of osteoclasts and osteoblasts [35, 36, 42, 43]. In what follows, PTH and PRL are then assumed to equilibrate rapidly to their equilibrium for which dx/dt=0 and dy/dt=0, respectively. That is,(3)x=a1b2Ak1+z+a2a3Ak2+z≡f1zand(4)y=a3b2k2+z+a4f1zb2=a3b2k2+z+a1a4A0k1+z+a2a3a4Ab2k2+z≡f2z,where A=b1b2−a2a4 with b1b2>a2a4.Thus, the reduced system of (1) and (2) can be obtained as follows:(5)dzdt=a5f1zk3+f12z+a6f2zk4+f22zzw−b3z,dwdt=a7f1z−a8f1zwk5+f1z−a9f2zwk6+f2z−b4w,t≠nT,with(6)Δzt=−ρzt,Δwt=μ,t=nT.


## 3. Stability Analysis

Now, we let(7)U:R+×R+2⟶R+,where R+=0,∞,R+2=S∈R2:S=z,w,z,w∈R+. Let us denote the map defined by the right-hand side of (5) and (6) by F=F1,F2.Definition 1.
LetU be a continuous function in nT,n+1T×R+2⟶R+. If(8)limt,Y⟶nT+,SUt,Y=UnT+,Sexists and is finite for each S∈R+2,n∈Z+, and V being locally Lipschitzian in S, then the upper right derivative of Ut,S with respect to (5) and (6) is defined as(9)D+Ut,S=limh⟶0+sup1hUt+h,S+hFt,S−Ut,S,where F=F1,F2.
The solutionSt=zt,wt to systems (5) and (6) is assumed to be a piecewise continuous function which means that St is continuous on nT,n+1T for n∈Z+ and limt⟶nT+St=SnT+ exists. According to the smoothness properties of F [44], the solution to systems (5) and (6) exists and is unique. For t≠nT, we can see that dz/dt=0 if zt=0 and dw/dt>0 if wt=0. Furthermore, for t=nT,ZnT+=1−ρznT, where 0<ρ<1 and wnT+=wnT+, where μ>0.Lemma 1.
For allt≥0, the solution of systems (5) and (6), St=zt,wt, is nonnegative provided that S0+≥0. In addition, St is positive for all t≥0 provided that S0+>0.Lemma 2.
Suppose thatzt,wt is a solution to (5) and (6) and(10)b4>a5M3k3+a6M4k4.
Then, for sufficiently larget, there exists a positive constant M such that zt≤M and wt≤M.Proof.
Letut=zt+wt,supf1z=M1,supzf1z=1/Aa1b2+a2a3=M3, supzf2z=a3/b2+a1a4/A+a2a3a4/Ab2=M4, and c=minb3,b4−a5M3/k3−a6M4/k4.
Fort≠nT, it follows that(11)D+u+cu=dzdt+dwdt+cz+cw=a5f1zk3+f12z+a6f2zk4+f22zzw−b3z+a7f1z−a8f1zwk5+f1z−a9f2zwk6+f2z−b4w+cz+cw≤a5zf1zk3+f12z+a6zf2zk3+f12zw+c−b3z+c−b4w+a7M1≤a5M3k3+a6M4k4+c−b4w+c−b3z+a7M1≤a7M1≡M0.
Thus,D+u≤−cu+M0. For t=nT, we can see that(12)unT+=znT++wnT+=1−ρznT+μ+wnT=znT+wnT+μ−ρznT≤unT+μ.
By applying Lemma 2.2 in [44], it can be derived that(13)ut≤u0e∫0t−cds+∫0te∫st−cdτM0ds+∑0<tn<te∫tnt−cdτ=u0e−ct+M0∫0te−ct−sds+∑0<tn<tμe−ct−tn=u0e−ct+M01c−e−ctc+μe−ct−T−e−ct−tn+11−ecT<M0c+μe−ct−T1−ecT+μecTecT−1≤M0c+μecTecT−1≡M.
Therefore,ut≤M as t⟶∞ and ut is uniformly ultimately bounded, which means that when t is large enough with zt≤M and wt≤M for some positive constant M.
Let us consider the following reduced impulsive system of (5) and (6) when osteoclasts are absent z=0:(14)dwdt=B−Cw,t≠nT,wnT+=wnT+μ,t=nT,w0+=w0,where B=a7f10,C=a8f10/k5+f10+a9f20/k6+f20+b4,f10=a1b2/Ak1+a2a3/Ak2 and f20=a3/b2k2+a4f10/b2.
It is obvious thatB and C are positive provided that b1b2>a2a4. By solving system (14), we obtain its periodic solution as follows:(15)w˜t=μe−Ct−nT1−e−CT+BC,t∈nT,n+1T,such that w˜0+=μ/1−e−CT+B/C>0.
Therefore, the positive solution to system (14) can be written as(16)wt=w0−BC−μ1−e−CTe−Ct+w˜t,t∈nT,n+1T.Lemma 3.
System (14) has a positive periodic solution w˜t. In addition, the solution wt of (14) converges to w˜t as t⟶∞.
Hence, the positive periodic solution to systems (5) and (6) in the absence of osteoclasts is(17)0,w˜t=0,μe−Ct−nT1−e−CT+BC,for t∈nT,n+1T with w˜nT+=w˜0+=μ/1−e−CT+B/C,n∈Z+.Theorem 1.
The solution0,w˜t of systems (5) and (6) is locally asymptotically stable if(18)0<T<Tmax,(19)b3<DBC,(20)DμC<ln11−ρ,where D=a5f10/k3+f120+a6f20/k4+f220 and Tmax≡1/DB/C−b3ln1/1−ρ−Dμ/C.Proof.
The local stability of the solution0,w˜t may be determined by considering the behavior of small amplitude perturbations of the solution.
Here, we let(21)zt=u1t,wt=w˜t+u2t.
It follows that(22)u1tu2t=Φtu10u20,0<t<T,where Φtis the fundamental solution matrix, in which(23)dΦtdt=Dw˜t−b30∗−CΦt,and Φ0is the identity matrix I.
Thus, we obtain the fundamental solution matrix as follows:(24)Φt=e∫0tDw˜s−b3ds0∗∗e∫0t−Cds.
Note that the terms∗ and ∗∗ are not required in further analysis and their exact expressions are not necessary.
Linearization of (6) provides(25)u1nT+u2nT+=1−ρ001u1nTu2nT.
Next, the Floquet theory is then applied to guarantee the local stability of the solution0,w˜t.
Let us consider(26)M=1−ρ001ΦT.
The solution0,w˜t is locally asymptotically stable if the magnitudes of λ1 and λ2 and the eigenvalue of M are both less than 1.
We can see that(27)λ1=1−ρe∫0TDw˜s−b3ds=1−ρeDμ/C+DBT/C−b3T,λ2=e∫0T−Cds=e−CT.
Since0<ρ<1andC>0, it follows that λ1<1 provided that (18)–(20) are satisfied. Also, it is obvious that λ2<1, and hence, the proof is complete.


## 4. Permanence of the System

Definition 2.
Systems (5) and (6) are said to be permanent if there exist positive constants m,M, and t0 such that for all solutions zt,wt with positive initial values z0+ and w0+ the following conditions hold for all t>t0>0:(28)m≤zt≤M,(29)m≤wt≤M.Theorem 2.
Suppose that(30)T>T∗,where T∗≡1/DB/C−b3ln1/1−ρ. Systems (5) and (6) are permanent provided that (10), (19), and (30) hold.Proof.
LetSt=zt,wt be any solution of (5) and (6) with z0+>0 and w0+>0. Suppose that (10) holds; it follows that when t is large enough, zt≤M and wt≤M for some M>0 according to Lemma 2.
Consider the second equation of (5); we can see that(31)dwdt≥−a8f1zwk5+f1z−a9f2zwk6+f2z−b4w≥−Cw,t≠nT,wnT+=wnT+μ,t=nT.
It follows that for sufficiently larget, there exists ε>0 such that(32)wt>−BC+w˜t−ε.
Therefore, we obtain(33)wt>μe−CT1−e−CT−ε≡m1,when t is large enough.
Next, we will show that there is a positive constantm2 for which zt>m2. Firstly, for some m3>0, we let(34)M1^=a5f1m3k3+f12m3+a6f2m3k4+f22m3,M2^=a1a7b2k2+a2a3a7k1CAk1+m3k2+m3.
Secondly, we will do the following two steps.Step 1.
We will show thatzt1≥m3 for some t1>0 by contradiction.
Suppose thatzt<m3 for all positive values of t. From the second equation of (5) and (6), we can see that(35)dwdt=a7f1z−a8f1zwk5+f1z−a9f2zwk6+f2z−b4w,t≠nT,≥a1a7b2k2+zAk1+zk2+z+a2a3a7k1+zAk1+zk2+z−Cw,≥a1a7b2k2+a2a3a7k1Ak1+m3k2+m3−Cw,wt+=wt+μ,t=nT.
Let us consider the following comparison system:(36)dPdt=a1a7b2k2+a2a3a7k1Ak1+m3k2+m3−CP,t≠nT,Pt+=Pt+μ,t=nT,P0+=w0+.
A periodic solution of (29) can then be obtained as follows:(37)P˜t≡μe−Ct−nT1−e−CT+a1a7b2k2+a2a3a7k1CAk1+m3k2+m3,t∈nT,n+1T,with P˜0+≡μ/1−e−CT+B/C. Therefore, the positive solution of (29) is(38)Pt=P0+−a1a7b2k2+a2a3a7k1CAk1+m3k2+m3−μ1−e−CTe−CT+P˜t.
Note thatf1m3=a1b2/Ak1+m3+a2a3/Ak2+m3 and f2m3=a3/b2k2+m3+a4/b2f1m3.
From the first equation of reduced system (5), we have(39)dzdt=a5f1zk3+f12z+a6f2zk4+f22zzw−b3z≥a5f1m3k3+f12m3+a6f2m3k4+f22m3w−b3z=M1^w−b3z.
According to the comparison theorem [45], we obtain wt≥Pt. It follows that(40)P˜t−ε1<Pt≤wt,t≠nT,t≥T1,for some T1>0 and for ε1>0 which is small enough.
Hence,(41)dzdt≥M1^P˜t−ε1−b3z,t≠nT,t≥T1,zt+=1−ρzt,t=nT,t≥T1.
Suppose thatn∈Z+ and NT≥T1. Then, by integration over nT,n+1T,n≥N, we obtain(42)zn+1T≥znT1−ρe∫nTn+1TM1^P˜t−ε1−b3dt=znT1−ρeM1^M2^−M1^ε1−b3T+M1^μ/C=znTη,where η≡1−ρeM1^M2^−M1^ε1−b3T+M1^μ/C.
We can see that(43)lnη=ln1−ρ+M1^M2^−M1^ε1−b3T+M1^μC.
Then, for a positive constantε1 which is sufficiently small, we have(44)lnη≈ln1−ρ+M1^M2^−M1^ε1−b3T+M1^μC>ln1−ρ+M1^M2^−b3T=M1^M2^−b3T−ln11−ρ.
Therefore, if (19) and (30) hold, then a small positive constant m3 can be chosen so that lnη>0, and hence η>1. It follows that zn+kT≥znTηk⟶∞ when k⟶∞ which contradicts the boundedness of zt. Hence, there exists t1>0 such that zt1≥m3.Step 2.
Ifzt≥m3 for all t≥t1, then the proof is complete; otherwise, there exists t>t1 such that zt<m3 and we then let t∗=inft>t1zt<m3, so that the following two possible cases are obtained.Case 1.
There existsn1∈Z+, such that t∗=n1T. It follows that for t∈t1,t∗,zt≥m3, and we obtain zt=m3 by the continuity of zt. When t is large enough, zt and wt are both bounded above by a positive constant M and wt is also bounded below by a positive constant m1 which imply that we can choose positive constants M′ and m′ for which zt<M′ and m1′<wt<M′ such that(45)m1′<b3M1^,with(46)wt∗+−a1a7b2k2+a2a3a7k1CAk1+m3k2+m3−μ1−e−CT−μ<M1′.
We also choosen2,n3∈Z+ that satisfy the following conditions:(47)n2T>1BlnM′+με1and(48)1−ρn2eη1n2+1Tηn3>1,where(49)η1≡M1^m1′−b3<0.
Here, we letT′=n2T+n3T. We claim that there exists a constant t2∈t∗,t∗+T′ in which zt2>m3. Otherwise, by considering (38) with Pt∗+=wt∗+, for t∈nT,n+1T and n1≤n≤n1+n2+n3, we obtain(50)Pt=Pt∗+−a1a7b2k2+a2a3a7k1CAk1+m3k2+m3−μ1−e−CTe−Ct−t∗+P˜t.
Forn2T≤t−t∗≤T′,(51)Pt−P˜t=Pt∗+−a1a7b2k2+a2a3a7k1CAk1+m3k2+m3−μ1−e−CTe−Ct−t∗=wt∗+−a1a7b2k2+a2a3a7k1CAk1+m3k2+m3−μ1−e−CTe−Ct−t∗<M′+μe−Ct−t∗<M′+μe−Bn2T<ε1.
It follows that(52)P˜t−ε1<Pt<wt.
Similar to Step1, we then obtain(53)zt∗+T′=zn1T+n2T+n3T≥zt∗+n2Tηn3.
According to the first equation of (5),(54)dzdt=a5f1zk3+f12z+a6f2zk4+f22zzw−b3z,t≠nT≥M1^w−b3z=M1^m1′−b3z=η1z,zt∗=1−ρzt,t=nT.
By integrating (54) over t∗,t∗+n2T, it follows that(55)zt∗+n2T≥zt∗eη1n2T≥zt∗1−ρeη1n2T≥m31−ρn2eη1n2T≥m31−ρn2eη1n2+1T.
Therefore, we obtain(56)zt∗+T′≥zt∗+n2Tηn3≥m31−ρn2eη1n2+1Tηn3>m3,which is a contradiction, and hence, zt2>m3 for some t2∈t∗,t∗+T′.
Next, we lett˜=inft>t∗zt>m3. It follows that zt<m3 for t∈t∗,t˜, and we can obtain zt˜=m3 by the continuity of zt. Then, l∈Z+ is chosen where l≤n2+n3 and t∗+lT≥t˜. Suppose t∈t∗+l−1T,t∗+lT. From (54), we obtain(57)zt≥zt∗+l−1T1−ρl−1eη1t−t∗+l−1T≥m31−ρleη1t−t∗≥m31−ρn2+n3eη1lT≥m31−ρn2+n3eη1n2+n3T.
Sinceη1 is negative and l<n2+n3, we let(58)m2¯=m31−ρn2+n3eη1n2+n3T,so that zt≥m2¯ for t∗,t˜. We can use t˜ instead of t∗, the proof can be proceeded in the same manner and consequently, we will obtain zt≥m2¯ for sufficiently large t.Case 2.
t∗≠nT for all n∈Z+. Then, for t∈t1,t∗, we have zt≥m3 and zt∗=m3. We assume that there exists n1′∈Z+ such that t∗∈n1′T,n1′+1T. This leads to the following two possible subcases:Case 2.1: for allt∈t∗,n1′+1T,zt≤m3. If there exists a constant t2′∈n1′+1T,n1′+1T+T′, then we can claim that zt2′>m3. Otherwise, let us consider (38) with Pn1′+1T+=wn1′+1T+. For t∈nT,n+1T,n1′+1≤n≤n1′+1+n2+n3, it follows that(59)Pt=Pn1′+1T+−a1a7b2k2+a2a3a7k1CAk1+m3k2+m3−μ1−e−CTe−Ct−n1′+1T+P˜t.By a similar argument in Case1, we have(60)Pt−P˜t<ε1,forn2T≤t−t∗. It follows that(61)P˜t−ε1<Pt≤wt.Sincen2T≤n1′+1+n2T−t∗, we then have(62)zn1′+1+n2T≥zt∗eη1n1′+1+n2T−t∗≥zt∗1−ρn2eη1n1′+1+n2T−t∗≥m31−ρn2eη1n2+1T.Thus,(63)zn1′+1+n2+n3T≥zn1′+1+n2Tηn3≥m31−ρn2eη1n2+1Tηn3>m3,which is a contradiction. Therefore, we obtainzt2′>m3 for some t2′∈n1′+1T,n1′+1T+T′.Next, we assume thatt¯=inft>t∗zt>m3. Then, we have zt≤m3 for t∈t∗,t¯ and zt¯=m3. Suppose that t∈n1′T+l′−1T,n1′T+l′T and a positive number l′ is chosen for which l′≤n2+n3+1. From (54), we obtain(64)zt≥zn1′T+l′−1Teη1t−n1′T+l′−1T,zn1′T+l′−1T≥zt∗eη1n1′T+l′−1T−t∗.Therefore,(65)zt≥zn1′T+l′−1T1−ρl−1eη1t−n1′T+l′−1T≥zt∗1−ρl−1eη1t−t∗≥m31−ρl′−1eη1t−t∗.Sincen1′T+l′−1T<t≤n1′T+l′T and n1′T<t∗<n1′+1T, we then obtain t−t∗≤l′T<n2+n3+1T. Thus,(66)zt≥m31−ρn2+n3eη1n2+n3+1T.We let(67)m2=m31−ρn2+n3eη1n2+n3+1T.It follows thatzt≥m2 for t∈t∗,t¯. We can use t¯ instead of t∗, the proof can be proceeded in the same manner, and consequently, we will obtain zt≥m2 for sufficiently large t.Case 2.2: there existst″∈t∗,n1′+1T such that zt″>m3. We assume that t¯=inft>t∗zt>m3. Then, we have zt<m3 for t∈t∗,t¯, and zt¯=m3. Since t<n1′+1T<t∗+T, for t∈t∗,t¯, we obtain the following:(68)zt≥zt∗eη1t−t∗≥m3eη1T=m2.
Sincezt¯≥m3, we can proceed the proof in the same manner for t>t¯. We then have zt≥m2 for t≥t1 because m2<m2¯<m3 and the proof is complete.


## 5. Periodic Solution: Existence and Stability

To investigate the possibility of a bifurcation of the nontrivial periodic solution to systems (5) and (6) near 0,w˜t, we first interchange the variables of z and w for convenience on computation. In what follows, let us consider the following system instead:(69)dzdt=a7f1w−a8f1wzk5+f1w−a9f2wzk6+f2w−b4z,dwdt=a5f1wk3+f12w+a6f2wk4+f22wzw−b3w,t≠nT,with(70)Δzt=μ,Δwt=−ρwt,t=nT.Let(71)G1z,w=a7f1w−a8f1wzk5+f1w−a9f2wzk6+f2w−b4z,G2z,w=a5f1wk3+f12w+a6f2wk4+f22wzw−b3w.Note that(72)f1w=a1b2Ak1+w+a2a3Ak2+w,f1′w=−a1b2Ak1+w2+a2a3Ak2+w2,f2w=a3b2k2+w+a1a4Ak1+w+a2a3a4Ab2k2+w,f2′w=−a3b2k2+w2+a1a4Ak1+w2+a2a3a4Ab2k2+w2.According to Lakmeche and Ario [45], the following notations are used:(73)∂Φ1∂w=∂Φ2∂z=0,∂Φ1∂z=1,∂Φ2∂w=1−ρ,∂2Φ1∂z∂w=∂2Φ2∂z∂w=∂2Φ2∂w2,Θ1z,w=z+μ,Θ2z,w=1−ρw,ζt=w˜t,0T,S0=w˜τ0,0T,τ0=Tmax.And, we have(74)∂G1ζr∂z=−a8f10k5+f10−a9f20k6+f20−b4=−C,∂G2ζr∂w=a5f10k3+f120+a6f20k4+f220w˜r−b3=Dw˜r−b3,∂G1ζr∂w=a7f1′0−a8w˜rf1′0k5+f10−f1′0f10k5+f102−a9w˜rf2′0k6+f20−f2′0f20k6+f202,∂G2ζr∂z∂w=a5f10k3+f120+a6f20k4+f220=D,∂G2ζr∂w2=2a5w˜rf1′0k3+f120−2f1′0f120k3+f1202+2a6w˜rf2′0k4+f220−2f2′0f220k4+f2202,∂Φ1τ0,S0∂τ=∂w˜τ0,S0∂t=−Cμexp−Cτ01−exp−Cτ0<0,∂Φ1τ0,S0∂z=exp∫0τ0∂G1ζr∂zdr=exp−Cτ0>0,∂Φ1τ0,S0∂w=∫0τ0exp∂G1ζr∂zdr∂G1ζν∂wexp∫0ν∂G2ζr∂wdrdν=∫0τ0exp−Cτ0−ν×a7f1′0−a8w˜νf1′0k5+f10−f1′0f10k5+f102−a9w˜νf2′0k6+f20−f2′0f20k6+f202×exp∫0vDw˜r−b3drdν,∂Φ2τ0,S0∂w=exp∫0τ0∂G2ζr∂wdr=exp∫0τ0Dw˜r−b3dr>0,∂2Φ2τ0,S0∂z∂w=∫0τ0exp∫vτ0∂G2ζr∂wdr∂2G2ζν∂z∂wexp∫0v∂G2ζr∂wdrdv=∫0τ0exp∫ντ0Dw˜r−b3drDexp∫0vDw˜r−b3drdν>0,∂2Φ2τ0,S0∂w2=∫0τ0exp∫vτ0∂G2ζr∂wdr∂2G2ζυ∂w2exp∫0v∂G2ζr∂wdrdν+∫0τ0exp∫vτ0∂G2ζr∂wdr∂2G2ζr∂z∂w×∫0vexp∫θv∂G1ζr∂zdr∂G1ζθ∂wexp∫0θ∂G2ζr∂wdrdθdν,∂2Φ2τ0,S0∂w∂τ=∂G2ζτ0∂wexp∫0τ0∂G2ζr∂wdr=Dw˜τ0−b3exp∫0τ0Dw˜r−b3dr.Now, we can compute(75)d0′=1−∂Θ2∂w∂Φ2∂wτ0,S0=1−1−ρexp∫0τ0Dw˜r−b3dr,where τ0 is the root of d0′=0 and(76)a0′=1−∂Θ1∂z∂Φ1∂zτ0,S0=1−exp−Cτ0>0,b0′=−∂Θ1∂z∂Φ1∂w+∂Θ1∂w∂Φ2∂wτ0,S0=−∂Θ1τ0,S0∂w,I=−∂2Θ2∂z∂w∂Φ1τ0,S0∂τ+∂Φ1τ0,S0∂z1a0′∂Θ1∂z∂Φ1τ0,S0∂τ∂Φ2τ0,S0∂w−∂Θ2∂w∂2Φ2τ0,S0∂z∂w1a0′∂Θ1∂z∂Φ1τ0,S0∂τ+∂Φ2τ0,S0∂τ∂w=−1−ρ∂2Φ2τ0,S0∂z∂w11−exp−Cτ0∂Φ1τ0,S0∂τ+Dw˜τ0−b3exp∫0τ0Dw˜r−b3dr=−1−ρexpDμC+DBτ0C−b3τ0−CDτ0μexp−Cτ01−exp−Cτ02+Dw˜τ0−b3,J=−2∂2Θ2∂z∂w−b0′a0′∂Φ1τ0,S0∂z+∂Φ1τ0,S0∂w∂Φ2τ0,S0∂w−∂2Θ2∂w2∂Φ2τ0,S0∂w2+2∂Φ2∂wb0′a0′∂2Φ2τ0,S0∂w∂z−∂Θ2∂w∂2Φ2τ0,S0∂w2=21−ρb0′a0′∂2Φ2τ0,S0∂w∂z−1−ρ∂2Φ2τ0,S0∂w2.We can see thatI>0 if(77)T>1Cln1+μCDb3C−BD,(78)μ<BC−b3D,while J<0 provided that (19) and(79)BCK1+K2>K0,(80)k3<f120,(81)k4<f220,hold where K0=a7f1′0,K1=−a8k5f1′0/k5+f102 and K2=−a9k6f2′0/k6+f202.This leads to the following theorem according to Lakmeche and Ario [45].Theorem 3.
If (19), (77), and (78) hold, then systems (5) and (6) have a positive periodic solution.


## 6. Numerical Simulations

In this section, numerical simulations are performed to illustrate the theoretical predictions as proved in the previous section.To illustrate the theoretical prediction in Theorem1, the parameters in systems (5) and (6) are chosen to satisfy the conditions in (18) and (20) as follows: a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1, a6=0.34,a7=0.91,a8=0.7,a9=0.03, b1=0.98,b2=0.68,b3=0.25,b4=0.15,k1=0.56,k2=0.95,k3=0.96, k4=0.99,k5=0.84,k6=0.65,μ=0.01,ρ=0.2,T=1,z0=0.1, and w0=0.1. The computer simulation is shown in Figure 1. We can see that the solution of systems (5) and (6) converges asymptotically to the oscillating solution 0,w˜t as time passes for which the active osteoblasts are vanished as predicted in Theorem 1.Figure 1
A computer simulation of systems (5) and (6). The solution trajectory approaches oscillatory solution 0,w˜t when time passes. Here, a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1, a6=0.34,a7=0.91,a8=0.7,a9=0.03, b1=0.98,b2=0.68,b3=0.25,b4=0.15,k1=0.56,k2=0.95, k3=0.96,k4=0.99,k5=0.84,k6=0.65, μ=0.01,ρ=0.2,T=1,z0=0.1, and w0=0.1 in which all conditions in Theorem 1 are satisfied.To illustrate the theoretical prediction in Theorem2, the parameters in systems (5) and (6) are chosen to satisfy the conditions in (10), (19), and (30) as follows: a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1, a6=0.34,a7=0.91,a8=0.7,a9=0.03, b1=0.98,b2=0.68,b3=0.01,b4=0.5,k1=0.56,k2=0.95,k3=0.96, k4=0.99,k5=0.84,k6=0.65,μ=0.1,ρ=0.1,T=5,z0=0.1, and w0=0.1. The computer simulation is shown in Figure 2. We can see that the solution of systems (5) and (6) are bounded above and bounded below, and hence, the system shows permanence as predicted in Theorem 2.Figure 2
A computer simulation of systems (5) and (6) showing the permanence of the system. Here, a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1,a6=0.34, a7=0.91,a8=0.7,a9=0.03,b1=0.98,b2=0.68,b3=0.01, b4=0.5,k1=0.56,k2=0.95,k3=0.96,k4=0.99, k5=0.84,k6=0.65,μ=0.1,ρ=0.1,T=5,z0=0.1, and w0=0.1 in which all conditions in Theorem 2 are satisfied.To illustrate the theoretical prediction in Theorem3, the parameters in systems (5) and (6) are chosen to satisfy the conditions in (29), (77), (78), and (79)-(81) as follows: a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1,a6=0.34,a7=0.91, a8=0.7,a9=0.03,b1=0.98,b2=0.68,b3=0.01, b4=0.5,k1=0.56,k2=0.95,k3=0.96,k4=0.2,k5=0.84,k6=0.65, μ=0.1,ρ=0.1,T=10,z0=0.1, and w0=0.1. The computer simulation is shown in Figure 3. We can see that the solution of systems (5) and (6) exhibits sustained oscillations as predicted in Theorem 3.Figure 3
A computer simulation of systems (5) and (6). The solution exhibits sustained oscillations. Here, a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1,a6=0.34,a7=0.91, a8=0.7,a9=0.03,b1=0.98,b2=0.68,b3=0.01,b4=0.5,k1=0.56, k2=0.95,k3=0.96,k4=0.2,k5=0.84,k6=0.65,μ=0.1, ρ=0.1,T=10,z0=0.1, and w0=0.1, in which all conditions in Theorem 3 are satisfied.In addition, to illustrate the effect of different periods between each estrogen administration, computer simulations of systems (5) and (6) for T=1,5,10,15 when a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1,a6=0.34,a7=0.91, a8=0.7,a9=0.03,b1=0.98,b2=0.68,b3=0.25,b4=0.15,k1=0.56, k2=0.95,k3=0.96,k4=0.99,k5=0.84,k6=0.65, μ=0.01,ρ=0.2,z0=0.5,andw0=0.5 are as shown in Figure 4.Figure 4
Computer simulations of systems (5) and (6) for T=1,5,10,15 when a1=0.32,a2=0.58,a3=0.19,a4=0.39,a5=0.1,a6=0.34,a7=0.91, a8=0.7,a9=0.03,b1=0.98,b2=0.68,b3=0.25,b4=0.15,k1=0.56, k2=0.95,k3=0.96,k4=0.99,k5=0.84,k6=0.65,μ=0.01, ρ=0.2,z0=0.5,andw0=0.5.As shown in Figure4, we can see that the period between each estrogen administration T is significant. The smaller value of T yields the better result on the number of osteoclasts and osteoblasts.


## 7. Conclusion

In this paper, an impulsive mathematical model of the process of bone remodeling accounted for bone-resorbing cells, osteoclasts, and bone-forming cells, osteoblasts, is developed in order to investigate the effect of impulsive estrogen supplement. The effects of parathyroid hormone and prolactin are also taken into account. We then apply the Floquet theory and the comparison theorem to derive the conditions in which the periodic solution is locally asymptotically stable. Moreover, the permanence of the system is also investigated as well so that we arrived at the conditions for which the sustained oscillation of the solution is guaranteed. In addition, computer simulations are presented to illustrate the theoretical predictions. The results indicate that the dosage of estrogen supplement indicated byμ and ρ and the frequency of estrogen supplement indicated by 1/T play important roles in the treatment of osteoporosis patients. Even though the smaller value of T yields the better result on the number of osteoclasts and osteoblasts, the side effects of estrogen administration are also needed to take into account. Therefore, the appropriate dosage and frequency of estrogen supplement that could control the number of bone-forming cells and bone-resorbing cells to lie within the desirable range might lead to an efficient treatment in osteoporosis patients.



---

*Source: 5435876-2021-12-15.xml*