Effects of prolactin on bone remodeling process with parathyroid hormone supplement: an impulsive mathematical model
Inthira Chaiya1 &
Chontita Rattanakul ORCID: orcid.org/0000-0002-2362-58851,2
In this paper, we investigate the effect of prolactin on bone remodeling process with the impulsive supplement of parathyroid hormone. An impulsive mathematical model is developed and analyzed theoretically in order to obtain the conditions on the parameters of the model for which the net bone formation might be expected from the parathyroid hormone supplement. Numerical investigation is also carried out to confirm our theoretical predictions. The dosage and the frequency of parathyroid hormone supplement turned out to be the keys for the effective parathyroid hormone supplement.
Osteoporosis is a metabolic bone disease detected mostly in postmenopausal women and is recognized as a major health problem in many countries. It is characterized mainly by low bone mass [1–3]. Nowadays, the population of those aged over 60 years is increasing rapidly in Thailand, and hence, osteoporosis is then becoming more prevalent in Thailand. Osteoporosis occurs when bone remodeling process, the process that normally occurs throughout life to control the repair or replacement of bone following injuries, is imbalanced with a net bone resorption greater than bone formation at bone multicellular units (BMUs) [1–3]. Bone remodeling process consists of bone resorption, for which osteoclasts are responsible, and bone formation, for which osteoblasts are responsible [1–3].
There are several hormones involved in bone remodeling process including prolactin (PRL). It has been reported that PRL-receptors have been found on osteoblasts [4]. Moreover, PRL has also been reported to enhance bone resorption in part by increasing receptor activator of NF-ligand (RANKL) and decreasing osteoprotegerin (OPG) expressions by osteoblasts [5]. On the other hand, the increase in the number of osteoclasts results in the increase in calcium level in blood and the secretion of parathyroid hormone (PTH), one of the principal hormones that regulate calcium levels in blood to vary in the normal ranges, from the parathyroid gland is then decreased [1–3]. The decrease in the PTH level then results in the reduction of PRL secretion from the anterior pituitary gland because PTH normally inhibits the reuptake and release of dopamine (DA) which is PRL inhibiting factor [1–3, 6].
Drug treatments in osteoporosis patients will slow down the activity of bone resorbing cells (osteoclasts) or stimulate the activity of bone forming cells (osteoblasts) or do both. PTH supplement is an option for treating osteoporosis patients. PTH has effects on both bone resorbing cells and bone forming cells [1–3]. As reported in [7], PTH at a high level inhibits the production of bone resorbing cells while it has both stimulating and inhibiting effects on the production of bone forming cells depending on the stage of the development of bone forming cells. In addition, the studies of [7–9] showed the paradoxical results of net bone formation and net bone loss when PTH is administered in a different manner. Hence, more insightful study on the effects of PTH supplements on bone remodeling process is needed.
Even though several mathematical models of bone remodeling process were proposed and analyzed theoretically and numerically [7, 10–14], the model that incorporates the effects of PRL and the impulsive treatments of PTH has not been proposed and analyzed yet. In the next section, we develop an impulsive mathematical model to investigate bone remodeling process based on the effects of PRL and the impulsive PTH supplement.
An impulsive mathematical model
Let \(x(t)\) denote the concentration of PRL above the basal level in blood at time t, \(y(t)\) denote the number of active osteoclasts at time t and \(z(t)\) denote the number of active osteoblasts at time t. We then propose an impulsive mathematical model to incorporate the effects of PRL and PTH-supplement on bone remodeling as follows:
where \(\Delta y(t)=y({{t}^{+}})-y(t)\), \(\Delta z(t)=z({{t}^{+}})-z(t)\), T represents the period of impulsive treatment of PTH, \(n\in {{Z}_{+}}\), \({{Z}_{+}}=\{1,2,3,\ldots\}\), ρ represents the inhibiting effect of PTH supplement on osteoclasts, \(0<\rho<1\), and μ represents the stimulating effect of PTH supplement on osteoblasts, \(\mu >0\). All parameters in the model are also assumed to be positive.
Equation (1a) represents the rate of change of PRL concentration in blood. As mentioned above, the increase in the number of active osteoclasts results in the increase in calcium level in blood, and the secretion of PTH from the parathyroid gland is then decreased [1–3]. The decrease in the PTH level then results in the reduction of PRL secretion from the anterior pituitary gland because of the inhibiting effect of PTH on the uptake and release of DA which is PRL inhibiting factor [1–3, 6]. Hence, the first term on the right-hand side accounts for the secretion rate of PRL from the lactotroph cells in the anterior pituitary gland which is inhibited by the increase in the number of active osteoclasts. The last term on the right-hand side represents the removal rate of PRL from the system.
Equation (1b) represents the rate of change of the number of active osteoclasts. Even though osteoclasts do not possess PRL receptors as osteoblasts as reported in [4], PRL has been reported to enhance bone resorption by increasing RANKL and decreasing OPG expressions by osteoblasts [5]. RANKL then binds to RANK, its receptor on preosteoclast, and stimulates proliferation and maturation into osteoclast. The first term on the right-hand side then accounts for the stimulating effect of PRL on the production of active osteoclasts. The last term on the right-hand side represents the removal rate of active osteoclasts from the system.
Equation (1c) represents the rate of change of the number of active osteoblasts. The first term on the right-hand side accounts for the production of active osteoblasts depending on the production and differentiation of active osteoclasts. That is, the more active osteoclasts, the more active osteoblasts are needed to refill the resorption cavities on BMUs created by active osteoclasts [1–3]. The second term on the right-hand side accounts for the inhibiting effect of PRL on the production of active osteoblasts which has been reported in [15]. The last term on the right-hand side accounts for the removal rate of active osteoblasts from the system.
Since the dynamics of PRL is very fast compared to the changes in the number of active osteoclasts and active osteoblasts [1–3], we then assume in what follows that PRL equilibrates quickly to its equilibrium for which \(\frac{dx}{dt}=0\). That is,
$$ x=\frac{{{c}_{1}}}{{{d}_{1}}({{k}_{1}}+y)}\equiv{{f}}(y). $$
System (1a)-(1e) is then reduced to the following system:
In the next section, we state the definitions and lemmas needed to prove the main results.
$$ V:{{R}_{+}}\times R_{+}^{2} \to{{R}_{+}}, $$
where \({{R}_{+}}=[0,\infty)\), \({R}_{+}^{2}=\{S\in{{R}^{2}}:S=(y,z), y\ge 0, z\ge0\}\). The map defined by the right-hand side of (3a)-(3b) is denoted by \(F=({{F}_{1}},{{F}_{2}})\).
Definition 3.1
The function V defined in (4) is said to belong to class \({V_{0}}\) if
V is continuous in \((nT,(n+1)T]\times R_{+}^{2}\to{R_{+}}\) and for each \(S\in R_{+}^{2}\), \(n\in{Z_{+}}\), \(\lim_{(t,Y)\to(nT^{+},S)} V(t,Y)=V(nT^{+},S)\) exists, and
V is locally Lipschitzian in S.
Suppose \(V\in{V_{0}}\). For \((t,S)\in(nT,(n+1)T]\times R_{+}^{2}\), the upper right derivative of \(V(t,S)\) with respect to (4)-(7) is defined by
$${D^{+}}V(t,S)=\limsup_{h \to0^{+}}\frac{1}{h} \bigl[V \bigl(t+h,S+hF(t,S) \bigr)-V(t,S) \bigr], $$
where \(F=({F_{1}},{F_{2}})\).
Assume that the solution of (3a)-(3d), \(S(t)=(y(t),z(t))\), is a piecewise continuous function. That is, \(S(t):{{R}_{+}}\to{R}_{+}^{2}\), \(S(t)\) is continuous on \((nT,(n+1)T]\), \(n\in{{Z}_{+}}\) and \(\lim_{t\to nT^{+}}S(t)=S(nT^{+})\) exists. Then the smoothness properties of F guarantee the global existence and uniqueness of solution to (3a)-(3d).
Since \(\frac{dy}{dt}=0\) whenever \(y(t)=0\), \(t\ne nT\), \(\frac {dz}{dt}>0\) whenever \(z(t)=0\), \(t\ne nT\) and \(y(n{{T}^{+}})=(1-\rho )y(nT)\), \(0<\rho<1\), \(z(n{{T}^{+}})=z(n{{T}})+\mu\), \(\mu>0\), the following lemmas are obtained.
Lemma 3.1
Suppose \(S(t)=(y(t),z(t))\) is a solution of (3a)-(3d) with \(S({{0}^{+}})\ge0\). Then \(S(t)\ge0\) for all \(t\ge0\).
There exists a constant \(M>0\) such that, for sufficiently large t, \(y(t)\le M\) and \(z(t)\le M\) provided that
$$ d_{3}>\frac{{c_{1}}{c_{2}}}{{d_{1}}{k_{2}}}, $$
where \((y(t),z(t))\) is a solution of (3a)-(3d).
Let \(v(t)=y(t)+z(t)\), \({M_{1}}=\sup yf(y)=\frac{c_{1}}{d_{1}}\), and \({M_{2}}= c_{3}\).
For \(t\ne nT\), we choose a positive constant c for which \(c=\min \{d_{2},{d}_{3}-\frac{{c_{1}}{c_{2}}}{{d_{1}}{k_{2}}} \}\).
$$\begin{aligned} D^{+}v+cv =&\frac{dy}{dt}+\frac{dz}{dt}+cy+cz \\ =& \frac{{c_{2}}f(y)yz}{k_{2}+f^{2}(y)}-{d_{2}}y+\frac{{c_{3}}y}{k_{3}+y}- \frac {{c_{4}}f(y)z}{k_{4}+f(y)}-{d_{3}}z+cy+cz \\ \le&(c-{d_{2}})y+ \biggl(\frac{{c_{2}}{M_{1}}}{k_{2}}-{{d}_{3}}+c \biggr)z+c_{3} \\ \le& c_{3} = {M_{2}}. \end{aligned}$$
Therefore, \({{D}^{+}}v\le-cv+{{M}_{2}}\).
For \(t=nT\),
$$\begin{aligned} v \bigl(n{{T}^{+}} \bigr) =&y \bigl(n{{T}^{+}} \bigr)+z \bigl(n{{T}^{+}} \bigr) \\ =&(1-\rho)y(nT)+z(nT)+\mu \\ =&y(nT)+z(nT)+\mu-\rho y(nT) \\ \le& v(nT)+\mu. \end{aligned}$$
According to Lemma 2.2 of Liu et al.[16], for \(t\in (nT,(n+1)T]\), we have
$$\begin{aligned} v(t) \le& v(0){e^{-ct}}+ \int _{0}^{t}{{M_{2}} {e^{-c(t-s)}} \,ds}+\sum_{0< {t_{n}}< t}{\mu} {e^{-c(t-{t_{n}})}} \\ \le& v(0){e^{-ct}}+{M_{2}} \biggl(\frac{1}{c}- \frac{{e^{-ct}}}{c} \biggr)+\mu \biggl(\frac {{e^{-c(t-T)}}-{e^{-c(t-{{t}_{n+1}})}}}{1-{e^{cT}}} \biggr) \\ < &\frac{{M_{2}}}{c}+\mu \biggl(\frac{{e^{cT}}}{{e^{cT}}-1} \biggr)\equiv M \quad \mbox{as } t\to\infty. \end{aligned}$$
Hence, \(v(t)\) is uniformly ultimately bounded and there exists a constant \(M>0\) such that \(y(t)\le M\) and \(z(t)\le M\) for sufficiently large t. □
Stability in the absence of active osteoclasts
Let us consider the impulsive system (3a)-(3d) in the absence of active osteoclasts \(y=0\), system (3a)-(3d) becomes
$$\begin{aligned}& \frac{dz}{dt} = -Az, \quad t\ne nT, \end{aligned}$$
$$\begin{aligned}& z \bigl(nT^{+} \bigr) = z(nT)+\mu, \quad t= nT, \end{aligned}$$
$$\begin{aligned}& z \bigl(0^{+} \bigr) = z_{0}, \end{aligned}$$
where \(A\equiv\frac{{c_{1}}{c_{4}}}{{d_{1}}{k_{1}}{k_{4}}}+d_{3}\). Here \(A>0\), a periodic solution of (6)-(8) is
$$ \tilde{z}(t)=\frac{\mu{e^{-A(t-nT)}}}{1-{e^{-AT}}}, \quad t\in\bigl(nT,(n+1)T\bigr] $$
with \(\tilde{z}(0^{+})=\frac{\mu}{1-{e^{-AT}}}\).
Hence, the positive solution of (6)-(8) is
$$z(t)= \biggl(z_{0}-\frac{\mu}{1-{e^{-AT}}} \biggr){{e}^{-At}}+ \tilde{z}(t), \quad t\in\bigl(nT,(n+1)T\bigr]. $$
This leads to the following result.
System (6)-(8) has a positive periodic solution \(\tilde{z}(t)\), and for every solution \(z(t)\) of (6)-(8), we have \(z(t)\to\tilde{z}(t)\) as \(t\to\infty\).
Therefore, system (3a)-(3d) has a periodic solution in the absence of active osteoclasts:
$$\bigl(0,\tilde{z}(t) \bigr)= \biggl(0,\frac{\mu {e^{-A(t-nT)}}}{1-{e^{-AT}}} \biggr) $$
for \(t\in(nT,(n+1)T]\) and \(\tilde{z}(nT^{+})=\tilde{z}(0^{+})=\frac{\mu }{1-{e^{-AT}}}\), \(n\in{Z_{+}}\). The conditions which guarantee the local stability of a periodic solution \((0,\tilde{z}(t))\) are then given in the following theorem.
Theorem 4.1
The solution \((0,\tilde{z}(t))\) of (3a)-(3d) is locally asymptotically stable if
$$ T>{T}_{\mathrm{max}} $$
$$ \ln \biggl(\frac{1}{1-\rho} \biggr)< \frac{D\mu}{A}, $$
$$ T_{\max}\equiv\frac{1}{d_{2}} \biggl( \frac{D\mu}{A}-\ln \biggl(\frac{1}{1-\rho } \biggr) \biggr) $$
and \(D\equiv\frac {{c_{1}}{c_{2}}{d_{1}}{k_{1}}}{c_{1}^{2}+d_{1}^{2}k_{1}^{2}{{k}_{2}}}\).
Consider a small perturbation from the point \((0,\tilde{z}(t))\)
$$\begin{aligned}& y(t) = {u_{1}}(t), \\& z(t) = \tilde{z}(t)+{u_{2}}(t). \end{aligned}$$
Then we may write
$$ \left ( \begin{matrix} {u_{1}}(t) \\ {u_{2}}(t) \end{matrix} \right )=\Phi(t) \left ( \begin{matrix} {u_{1}}(0) \\ {u_{2}}(0) \end{matrix} \right ), \quad0< t< T, $$
where \({\Phi}(t)\) satisfies
$$\frac{d{\Phi}(t)}{dt}=\left ( \begin{matrix} D\tilde{z}(t)-{d_{2}} & 0 \\ * &-A \end{matrix} \right ){\Phi}(t) $$
and \({\Phi}(0)=I\), the identity matrix.
Therefore, the fundamental solution matrix is
$${\Phi}(t)=\left ( \begin{matrix} \exp\int _{0}^{t}({D\tilde{z}(s)-{d_{2}}})\,ds & 0& \\ * & \exp\int _{0}^{t}{(-A)}\,ds \end{matrix} \right ). $$
Note that it is not necessary to find the exact expression for (*) because the term (*) is not involved in further analysis.
Linearization of (3c)-(3d) yields
$$\left ( \begin{matrix} {u_{1}}(n{{T}^{+}}) \\ {u_{2}}(n{{T}^{+}}) \end{matrix} \right )=\left ( \begin{matrix} 1-\rho& 0 \\ 0 & 1 \end{matrix} \right )\left ( \begin{matrix} {u_{1}}(nT) \\ {u_{2}}(nT) \end{matrix} \right ). $$
Hence, Floquet theory implies that the solution \((0,\tilde{z}(t))\) of (3a)-(3d) is locally stable if \(\vert \lambda_{i} \vert <1\), \(i=1,2\), where \(\lambda_{i}\) is an eigenvalue of
$$M=\left ( \begin{matrix} 1-\rho& 0 \\ 0 & 1 \end{matrix} \right ){\Phi}(T). $$
The eigenvalues of M are
$$\begin{aligned}& \lambda_{1}=(1-\rho)\exp \int _{0}^{T} \bigl({D\tilde {z}(s)-{d_{2}}} \bigr)\,ds=(1-\rho)\exp \biggl(\frac{D\mu}{A}-{d_{2}}T \biggr), \\& \lambda_{2}=\exp \int _{0}^{T}(-A)\,ds=\exp(-AT). \end{aligned}$$
Since \(0<\rho<1\), \(A>0\) and (10) holds, then
$$T>\frac{1}{d_{2}} \biggl(\frac{D\mu}{A}-\ln \biggl(\frac{1}{1-\rho} \biggr) \biggr). $$
$$\vert \lambda_{1} \vert =(1-\rho)\exp \biggl( \frac{D\mu }{A}-{d_{2}}T \biggr)< 1 $$
$$\vert \lambda_{2} \vert =\exp(-AT)< 1. $$
Therefore, Floquet theory implies that the solution \((0,\tilde{z}(t))\) is locally stable and the proof is complete. □
In the next section, the permanence of impulsive system (3a)-(3d) is investigated.
Permanence of the system
System (3a)-(3d) is said to be permanent if there are constants \(m, M>0\) (independent of the initial values) and a finite time \({{t}_{0}}\) such that for all solutions with initial values \(y(0^{+})>0\), and \(z(0^{+})>0\),
$$\begin{aligned}& m\le y(t)\le M, \\& m\le z(t)\le M, \end{aligned}$$
for all \(t>{t_{0}}\). Note that \({t_{0}}\) may depend on the initial values.
System (3a)-(3d) is permanent if (5) and (11) hold and
$$ T< {T}_{\mathrm{max}} . $$
Suppose that \(S(t)=(y(t),z(t))\) is a solution of system (3a)-(3d) with \(y(0^{+})>0\) and \(z(0^{+})>0\). Since (5) holds, Lemma 3.2 implies that there is a constant \(M>0\) such that, for sufficiently large t, \(y(t)\le M\) and \(z(t)\le M\).
Since \(\frac{{c_{3}}y}{{k_{3}}+y}>0\) and \((\frac {{c_{4}}f(y)z}{{k_{4}}+f(y)}+{d_{3}} )\) is a decreasing function when \(y>0\), (3b) implies that
$$\begin{aligned}& \frac{dz}{dt} \ge -Az, \quad t\ne nT, \\& z \bigl(nT^{+} \bigr) = z(nT)+\mu, \quad t= nT \end{aligned}$$
and then we have
$$z(t)>\tilde{z}(t)-\varepsilon $$
for some \(\varepsilon>0\) and for sufficiently large t.
$$z(t)>\frac{\mu{e^{-AT}}}{1-{e^{-AT}}}-\varepsilon\equiv{m_{1}} $$
for sufficiently large t.
Therefore, we only need to show that there exists a constant \({m_{2}}>0\) such that \(y(t)>{m_{2}}\). In order to do so, for some \({m_{3}}>0\), we first let
$$\hat{M}_{1}=\frac{{c_{1}}{c_{2}}{d_{1}}{k_{1}}}{{{c_{1}}^{2}}+{{d_{1}}^{2}}{k_{2}}({k_{1}}+{m_{3}})^{2}}. $$
Next, we do the following two steps.
Step 1 We prove by contradiction that there exists \(t_{1}\) such that \(y(t_{1})\ge{m_{3}}\). Suppose that \(y(t)<{m_{3}}\) for all positive t. From (3b) and (3d)
$$\begin{aligned}& \begin{aligned} \frac{dz}{dt}&=\frac{{c_{3}}y}{{k_{3}}+y}-\frac {{c_{4}}f(y)z}{{k_{4}}+f(y)}-{d_{3}}z, \quad t\ne nT \\ &\ge\frac{{c_{3}}(0)}{{k_{3}}+{m_{3}}}-\frac {{c_{4}}f(0)z}{{k_{4}}+f(0)}-{d_{3}}z \\ &= -Az, \end{aligned} \\& z \bigl(t^{+} \bigr)=z(t)+\mu, \quad t= nT. \end{aligned}$$
Let us consider the comparison system
$$\begin{aligned}& \frac{dP}{dt} = -AP, \quad t\ne nT, \end{aligned}$$
$$\begin{aligned}& P \bigl(t^{+} \bigr) = P(t)+\mu, \quad t= nT \end{aligned}$$
$$ P \bigl(0^{+} \bigr)=z \bigl(0^{+} \bigr). $$
$$ \tilde{P}(t)\equiv\frac{\mu{e^{-A(t-nT)}}}{1-{e^{-AT}}}, \quad t \in (nT,(n+1)T] $$
is a periodic solution of (14)-(16) with \(P(0^{+})\equiv \frac{\mu}{1-{e^{-AT}}}>0\). The positive solution of (14)-(16) is
$$ P(t)= \biggl(P \bigl(0^{+} \bigr)-\frac{\mu}{1-{e^{-AT}}} \biggr){e^{-At}}+ \tilde{P}(t), \quad t \in(nT,(n+1)T] $$
and \(P(t)\to\tilde{P}(t)=\frac{\mu{e^{-A(t-nT)}}}{1-{e^{-AT}}}\), as \(t\to\infty\).
The comparison theorem [17] then implies that \(z(t)\geq P(t)\).
From (3a), we have
$$\begin{aligned} \frac{dy}{dt} =&\frac{{c_{2}}f(y)yz}{{k_{2}}+f^{2}(y)}-{d_{2}}y \\ =& \biggl(\frac{{c_{1}}{c_{2}}{d_{1}} ({k_{1}}+y )z}{c_{1}^{2}+d_{1}^{2}{k_{2}}{{ ({k_{1}}+y )}^{2}}}-{d_{2}} \biggr)y \\ \ge& ({\hat{M}_{1}}z-{d_{2}} )y. \end{aligned}$$
Since \(z(t)\geq P(t)\), there is \({T_{1}}>0\) such that
$$ \tilde{P}(t)-{\varepsilon_{1}}< P(t)\le z(t),\quad t\ne nT, t\ge T_{1} $$
for sufficiently small \({\varepsilon_{1}} > 0\).
$$ \frac{dy}{dt}\ge \bigl({\hat{M}_{1}} \bigl(\tilde{P}(t)-{ \varepsilon_{1}} \bigr)-{d_{2}} \bigr)y, \quad t\ne nT, t \ge{T_{1}} $$
$$ y \bigl(t^{+} \bigr)=(1-\rho)y(t), \quad t= nT, t\ge T_{1}. $$
Letting \(N\in Z_{+}\) and \(NT\ge T_{1}\) and integrating over \((nT,(n+1)T]\), \(n\ge N\), we obtain
$$\begin{aligned} y \bigl((n+1)T \bigr) \ge& y(nT) (1-\rho)\exp \biggl( \int _{nT}^{(n+1)T}{ \bigl({\hat{M}_{1}} \bigl(\tilde{P}(t)-{\varepsilon_{1}} \bigr)-{d_{2}} \bigr) \,dt} \biggr) \\ =&y(nT) (1-\rho)\exp \biggl(\frac{{\hat{M}_{1}}\mu}{A}-({\varepsilon _{1}} { \hat{M}_{1}}+{d_{2}})T \biggr) \\ =&y(nT){\eta}, \end{aligned}$$
where \(\eta\equiv(1-\rho)\exp (\frac{{\hat{M}_{1}}\mu }{A}-({\varepsilon_{1}}{\hat{M}_{1}}+{d_{2}})T )\).
$$\ln\eta=\ln(1-\rho)+\frac{{\hat{M}_{1}}\mu}{A}-({\varepsilon_{1}} {\hat {M}_{1}}+{d_{2}})T. $$
For sufficiently small \({\varepsilon_{1}}>0\),
$$\begin{aligned} \ln\eta \approx&\ln(1-\rho)+\frac{{\hat{M}_{1}}\mu}{A}-{d_{2}}T \\ =&\frac{{\hat{M}_{1}}\mu}{A}-{d_{2}}T-\ln \biggl(\frac{1}{1-\rho} \biggr). \end{aligned}$$
Since \({\hat{M}_{1}}< D\), \(\ln (\frac{1}{1-\rho} )<\frac{D\mu }{A}\) and (11) hold, we can choose a small constant \({m_{3}}>0\) such that \(\eta>0\), and hence,
$$ \eta\equiv(1-\rho)\exp \biggl(\frac{{\hat{M}_{1}}\mu}{A}- ({\varepsilon _{1}} {\hat{M}_{1}}+{d_{2}} )T \biggr)>1. $$
Then \(y((n+k)T)\ge y(nT)\eta^{k}\to\infty\) as \(k\to\infty\), which contradicts the boundedness of \(y(t)\). Hence, there is \(t_{1}>0\) such that \(y(t_{1})\ge{m_{3}}\).
Step 2 If \(y(t)\ge{m_{3}}\) for all \(t>{t_{1}}\), then the proof is complete. Otherwise, \(y(t)<{m_{3}}\) for some \(t>{t_{1}}\). Let \(t^{*}=\inf_{t>t_{1}}\{y(t)< m_{3}\}\). There are two possible subcases as follows.
Case 1: \({t^{*}}={n_{1}}T\) for some \({n_{1}}\in{Z_{+}}\). This means \(y(t)\ge {m_{3}}\) for \(t\in(t_{1},t^{*}]\) and, by the continuity of \(y(t)\), we have \(y(t^{*})={m_{3}}\).
Since there are \(M>0\) and \(m_{1}>0\) such that \(y(t)< M\) and \(m_{1}< z(t)< M\) for sufficiently large t, we choose \(M'>0\) and \(m'_{1}>0\) such that
$$ y(t)< M' \quad\text{and} \quad m'_{1}< z(t)< M' $$
$$ m'_{1}< \frac{d_{2}}{\hat{M}_{1}} $$
such that
$$ \biggl\vert z \bigl(t^{*+} \bigr)-\frac{\mu}{1-e^{-AT}} \biggr\vert - \mu< M'. $$
Then choose \(n_{2}, n_{3}\in Z_{+}\) such that
$$ {n_{2}}T>\frac{1}{A}\ln \biggl(\frac{M'+\mu}{\varepsilon_{1}} \biggr) $$
$$ (1-\rho)^{n_{2}}\exp \bigl(({n_{2}}+1){\eta_{1}}T \bigr)\eta^{n_{3}}>1, $$
$$ \eta_{1}\equiv{\hat{M}_{1}}m'_{1}-d_{2}< 0. $$
Let \(T'={n_{2}}T+{n_{3}}T\). We claim that there is \(t_{2} \in(t^{*},t^{*}+T']\) such that \(y(t_{2})>m_{3}\). Otherwise, considering (18) with \(P(t^{*+})=z(t^{*+})\), we have
$$\begin{aligned} P(t)= \biggl(P \bigl(t^{*+} \bigr)-\frac{\mu}{1-e^{-AT}} \biggr)e^{-A(t-t^{*})}+\tilde{P}(t) \end{aligned}$$
for \(t\in(nT,(n+1)T]\) and \(n_{1}\le n\le n_{1}+n_{2}+n_{3}\).
For \({n_{2}}T\le t-t^{*}\le T'\), we have
$$\begin{aligned} \bigl\vert P(t)-\tilde{P}(t) \bigr\vert =& \biggl\vert P \bigl(t^{*+} \bigr)-\frac {\mu}{1-e^{-AT}} \biggr\vert e^{-A(t-t^{*})} \\ =& \biggl\vert z \bigl(t^{*+} \bigr)-\frac{\mu}{1-e^{-AT}} \biggr\vert e^{-A(t-t^{*})} \\ < & \bigl(M'+\mu \bigr)e^{-A(t-t^{*})} \\ < & \bigl(M'+\mu \bigr)e^{-A{n_{2}}T} \\ < &\varepsilon_{1}. \end{aligned}$$
$$ \tilde{P}(t)-\varepsilon_{1}< P(t)< z(t). $$
As in Step 1, we have
$$\begin{aligned} y \bigl(t^{*}+T' \bigr) =&y({n_{1}}T+{n_{2}}T+{n_{3}}T) \\ \ge& y \bigl(t^{*}+{n_{2}}T \bigr)\eta^{n_{3}}. \end{aligned}$$
$$\begin{aligned}& \begin{gathered} \frac{dy}{dt} = \frac{{c_{2}}f(y)yz}{{k_{2}}+f^{2}(y)}-{d_{2}}y, \quad t\ne nT \\ \hphantom{\frac{dy}{dt}} \ge({\hat{M}_{1}}z-d_{2})y \\ \hphantom{\frac{dy}{dt}} \ge \bigl({\hat{M}_{1}}m'_{1}-d_{2} \bigr)y \\ \hphantom{\frac{dy}{dt}} ={\eta_{1}}y, \\ y \bigl(t^{+} \bigr)=(1-\rho)y(t),\quad t= nT. \end{gathered} \end{aligned}$$
Integrating the above over \([t^{*},t^{*}+{n_{2}}T]\), we obtain
$$\begin{aligned}& \begin{aligned} y \bigl(t^{*}+{n_{2}}T \bigr)&\ge y \bigl(t^{*} \bigr){{(1-\rho)}^{n_{2}}}\exp \biggl( \int _{{n_{1}}T}^{{n_{1}}T+{n_{2}}T}{{\eta_{1}}\,dt} \biggr) \\ &\ge{m_{3}} {{(1-\rho)}^{n_{2}}}\exp({n_{2}} { \eta_{1}}T) \\ &\ge{m_{3}} {{(1-\rho)}^{n_{2}}}\exp \bigl(({n_{2}}+1){ \eta_{1}}T \bigr), \end{aligned} \\& \begin{aligned} y \bigl(t^{*}+T' \bigr)&\ge y \bigl(t^{*}+{n_{2}}T \bigr)\eta^{n_{3}} \\ &\ge{m_{3}}(1-\rho)^{n_{2}}\exp \bigl(({n_{2}}+1){ \eta_{1}}T \bigr)\eta^{n_{3}} \\ &>m_{3}, \end{aligned} \end{aligned}$$
which contradicts the definition of \(m_{3}\). Hence, there is \(t_{2} \in (t^{*},t^{*}+T']\) such that \(y(t_{2})>m_{3}\).
Now, let \(\tilde{t}=\inf_{t>t^{*}}\{y(t)>m_{3}\}\). Then \(y(t)< {m_{3}}\) for \(t\in(t^{*},\tilde{t})\), and by the continuity of \(y(t)\), we have \(y(\tilde{t})={m_{3}}\). We choose \(l\in{Z_{+}}\) such that \(l\le{n_{2}}+{n_{3}}\) and \(t^{*}+lT\ge\tilde{t}\), and suppose \(t\in (t^{*}+(l-1)T,t^{*}+lT]\). From (26), we have
$$\begin{aligned} y(t) \ge& y \bigl(t^{*+} \bigr) (1-\rho)^{l-1}\exp \bigl((l-1){ \eta_{1}}T \bigr)\exp \bigl({\eta _{1}} \bigl(t- \bigl(t^{*}+(l-1)T \bigr) \bigr) \bigr) \\ =&y \bigl(t^{*} \bigr) (1-\rho)^{l}\exp \bigl((l-1){ \eta_{1}}T \bigr)\exp \bigl({\eta_{1}} \bigl(t- \bigl(t^{*}+(l-1)T \bigr) \bigr) \bigr) \\ =&{m_{3}}(1-\rho)^{l}\exp \bigl({\eta_{1}} \bigl(t-t^{*} \bigr) \bigr) \\ \ge&{m_{3}}(1-\rho)^{n_{2}+n_{3}}\exp({\eta_{1}}lT) \\ \ge&{m_{3}}(1-\rho)^{n_{2}+n_{3}}\exp \bigl(({n_{2}+n_{3}}){ \eta_{1}}T \bigr). \end{aligned}$$
Since \(\eta_{1}<0\) and \(l\le n_{2}+n_{3}\), letting
$$ \bar{m}_{2}={m_{3}}(1-\rho)^{n_{2}+n_{3}}\exp \bigl((n_{2}+n_{3}){\eta_{1}}T \bigr), $$
we have \(y(t)\ge\bar{m}_{2}\) for \(t\in(t^{*},\tilde{t})\). We can continue in the same way by using t̃ instead of \(t^{*}\). Then we shall have \(y(t)\ge\bar{m}_{2}\) for all t large enough.
Case 2: \(t^{*}\ne nT\) for all \(n\in{Z_{+}}\). This means \(y(t)\ge{m_{3}}\) for \(t\in[t_{1},t^{*})\) and \(y(t^{*})={m_{3}}\). Suppose \(t^{*} \in ({n'_{1}}T,({n'_{1}}+1)T)\) for some \({n'_{1}}\in{Z_{+}}\). There are two possible subcases as follows.
Case 2.1: \(y(t)\le{m_{3}}\) for all \(t\in(t^{*},({n'_{1}}+1)T]\). We claim that there is \({t'_{2}}\in[({n'_{1}}+1)T,({n'_{1}}+1)T+T']\) such that \(y(t'_{2})>{m_{3}}\). Otherwise, consider (18) with \(P(({n'_{1}}+1){T^{+}})=z(({n'_{1}}+1){T^{+}})\). For \(t\in(nT,(n+1)T]\), \({n'_{1}}+1\le n\le{n'_{1}}+1+{n_{2}}+{n_{3}}\), we obtain
$$\begin{aligned} P(t) =& \biggl(P \bigl( \bigl(n'_{1}+1 \bigr)T^{+} \bigr)- \frac{\mu}{1-e^{-AT}} \biggr)e^{-A(t-(n'_{1}+1)T)}+\tilde{P}(t). \end{aligned}$$
Similar to Case 2.1, for \({n_{2}}T\le t-t^{*}\), we obtain
$$ \bigl\vert P(t)-\tilde{P}(t) \bigr\vert < \varepsilon_{1}. $$
$$ \tilde{P}(t)-\varepsilon_{1}< P(t)\le z(t). $$
Since \({n_{2}}T\le(n'_{1}+1+n_{2})T-t^{*}\), we have
$$\begin{aligned} y \bigl( \bigl(n'_{1}+1+n_{2} \bigr)T \bigr) \ge& y \bigl(t^{*} \bigr) (1-\rho)^{n_{2}}\exp \bigl({\eta _{1}} \bigl( \bigl(n'_{1}+1+n_{2} \bigr)T-t^{*} \bigr) \bigr) \\ \ge&{m_{3}}(1-\rho)^{n_{2}}\exp \bigl({\eta_{1}} \bigl( \bigl(n'_{1}+1+n_{2} \bigr)T-{n'_{1}}T \bigr) \bigr) \\ \ge&{m_{3}}(1-\rho)^{n_{2}}\exp \bigl((n_{2}+1){ \eta_{1}}T \bigr). \end{aligned}$$
$$\begin{aligned} y \bigl( \bigl(n'_{1}+1+n_{2}+n_{3} \bigr)T \bigr) \ge& y \bigl( \bigl(n'_{1}+1+n_{2} \bigr)T \bigr)\eta^{n_{3}} \\ \ge&{m_{3}}(1-\rho)^{n_{2}}\exp \bigl((n_{2}+1){ \eta_{1}}T \bigr)\eta^{n_{3}} \\ >&m_{3}, \end{aligned}$$
which contradicts the definition of \(m_{3}\). Hence, there is \(t'_{2} \in [(n'_{1}+1)T,(n'_{1}+1)T+T']\) such that \(y(t'_{2})>m_{3}\).
Now, let \(\bar{t}=\inf_{t>t^{*}}\{y(t)>m_{3}\}\). Then \(y(t)\le m_{3}\) for \(t\in[t^{*},\bar{t})\), and \(y(\bar{t})=m_{3}\). We choose \(l'\in Z_{+}\) such that \(l'\le n_{2}+n_{3}+1\) and suppose \(t\in ({n'_{1}}T+(l'-1)T,{n'_{1}}T+l'T]\). From (26), we have
$$\begin{aligned} y(t) \ge& y \bigl( \bigl({n'_{1}}T+ \bigl(l'-1 \bigr)T \bigr)^{+} \bigr)\exp \bigl(\eta_{1} \bigl(t- \bigl({n'_{1}}T+ \bigl(l'-1 \bigr)T \bigr) \bigr) \bigr) \\ =& y \bigl({n'_{1}}T+ \bigl(l'-1 \bigr)T \bigr) (1-\rho)\exp \bigl(\eta_{1} \bigl(t- \bigl({n'_{1}}T+ \bigl(l'-1 \bigr)T \bigr) \bigr) \bigr) \\ \ge& y \bigl(t^{*} \bigr) (1-\rho)^{l'-2}\exp \bigl(\eta_{1} \bigl( \bigl({n'_{1}}T+ \bigl(l'-1 \bigr)T \bigr)-t^{*} \bigr) \bigr) (1-\rho ) \\ &{}\times\exp \bigl(\eta_{1} \bigl(t- \bigl({n'_{1}}T+ \bigl(l'-1 \bigr)T \bigr) \bigr) \bigr) \\ =& {m_{3}}(1-\rho)^{l'-1}\exp \bigl(\eta_{1} \bigl(t-t^{*} \bigr) \bigr). \end{aligned}$$
Since \(\eta_{1}<0\) and \(t-t^{*}\le l'T\). Hence,
$$ y(t) \ge{m_{3}}(1-\rho)^{n_{2}+n_{3}}\exp \bigl((n_{2}+n_{3}+1){ \eta_{1}}T \bigr). $$
$$ m_{2} = {m_{3}}(1-\rho)^{n_{2}+n_{3}}\exp \bigl((n_{2}+n_{3}+1){\eta_{1}}T \bigr). $$
We have \(y(t)\ge m_{2}\) for \(t\in(t^{*},\bar{t})\). We can continue in the same manner by using t̄ instead of \(t^{*}\). Then we shall have \(y(t)\ge m_{2}\) for all t large enough.
Case 2.2: There is \(t''\in(t^{*},({n'_{1}}+1)T]\) such that \(y(t'')>{m_{3}}\). Let \(\underline{t}=\inf_{t>t^{*}}\{y(t)>m_{3}\}\). Hence, \(y(t)<{m_{3}}\) for \(t\in[t^{*},\underline{t})\), and \(y(\underline {t})={m_{3}}\). For \(t\in[t^{*},\underline{t})\), (26) holds, we have
$$\begin{aligned} y(t) \ge& y \bigl(t^{*} \bigr)\exp \biggl( \int _{t^{*}}^{t}{\eta_{1}}\,dt \biggr) \\ =& {m_{3}}\exp \bigl(\eta_{1} \bigl(t-t^{*} \bigr) \bigr) \\ \ge& {m_{3}}\exp({\eta_{1}}T) \\ >& m_{2} \end{aligned}$$
since \(t<{n'_{1}}T+T<t^{*}+T\). For \(t>\underline{t}\), we can continue in the same manner since \(y(\underline{t})\ge m_{3}\). Since \(m_{2}<\bar {m}_{2}<m_{3}\), we have \(y(t)\ge m_{2}\) for \(t \ge t_{1}\). The proof is complete. □
Existence and stability of the positive periodic solution
We now investigate the possibility of bifurcation of a positive periodic solution to system (3a)-(3d) near \((0,\tilde {z}(t))\). Let us exchange the state variables for convenience so that system (3a)-(3d) becomes
$$\begin{aligned}& \frac{dy}{dt} = \frac{{c_{3}}z}{k_{3}+z}-\frac {{c_{4}}yf(z)}{k_{4}+f(z)}-{d_{3}}y, \quad t\neq nT, \end{aligned}$$
$$\begin{aligned}& \frac{dz}{dt} = \frac{{c_{2}}yzf(z)}{k_{2}+f^{2}(z)}-{d_{2}}z, \quad t\neq nT \end{aligned}$$
$$\begin{aligned}& \Delta y(t) = \mu, \quad t=nT, \end{aligned}$$
$$\begin{aligned}& \Delta z(t) = -\rho z(t), \quad t=nT. \end{aligned}$$
$$F_{1}(y,z)=\frac{{c_{3}}z}{k_{3}+z}-\frac{{c_{4}}yf(z)}{k_{4}+f(z)}-{d_{3}}y, \qquad F_{2}(y,z)=\frac{{c_{2}}yzf(z)}{k_{2}+f^{2}(z)}-{d_{2}}z. $$
According to Lakmeche and Arini [18],
$$\begin{gathered} \Theta_{1}(y,z)=y+\mu,\qquad \Theta_{2}(y,z)=(1- \rho)z,\qquad \varsigma(t)= \bigl(\tilde {z}(t),0 \bigr)^{T}, \\ S_{0}= \bigl(\tilde{z}(\tau_{0}),0 \bigr)^{T}, \qquad \tau_{0}=T_{\mathrm{max}} \end{gathered} $$
$$\begin{aligned}& \frac{\partial{\Phi_{1}}(\tau_{0},S_{0})}{\partial\tau}=\frac{\partial\tilde {z}(\tau_{0},S_{0})}{\partial t}=\frac{-A\mu\exp(-A\tau_{0})}{1-\exp(-A\tau _{0})}< 0, \\& \frac{\partial{\Phi_{1}}(\tau_{0},S_{0})}{\partial y}=\exp \biggl( \int _{0}^{\tau_{0}}{\frac{{\partial F_{1}}(\varsigma(r))}{\partial y}\,dr} \biggr)= \exp(-A\tau_{0}), \\& \begin{aligned} \frac{\partial{\Phi_{1}}(\tau_{0},S_{0})}{\partial z}&= \int _{0}^{\tau _{0}}{\exp \biggl( \int _{\upsilon}^{\tau_{0}}{\frac{\partial F_{1}\varsigma(r)}{\partial y}\,dr} \biggr) \frac{\partial F_{1}( \varsigma (\upsilon))}{\partial z}}\exp \biggl( \int _{0}^{\upsilon}{\frac {\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr)\,d \upsilon \\ &= \int _{0}^{\tau_{0}}{\exp \bigl(-A(\tau_{0}- \upsilon) \bigr) \biggl( \frac {c_{3}}{k_{3}}-\frac{{c_{4}}{k_{4}}f'(0)\tilde{z}(\upsilon )}{({k_{4}}+f(0))^{2}} \biggr)} \\ &\quad {}\times\exp \biggl( \int _{0}^{\upsilon}{ \bigl(D\tilde {z}(r)-d_{2} \bigr)\,dr} \biggr)\,d\upsilon, \end{aligned} \\& \begin{aligned} \frac{\partial{\Phi_{2}}(\tau_{0},S_{0})}{\partial z}&=\exp \biggl( \int _{0}^{\tau_{0}}{\frac{\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr)=\exp \biggl( \int _{0}^{\tau_{0}}{ \bigl(D\tilde{z}(r)-d_{2} \bigr)\,dr} \biggr) \\ \frac{{\partial^{2}}\Phi_{2}(\tau_{0},S_{0})}{\partial y\,\partial z}&= \int _{0}^{\tau_{0}}{\exp \biggl( \int _{\upsilon}^{\tau_{0}}{\frac {\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr) \frac{{\partial ^{2}}F_{2}(\varsigma(\upsilon))}{\partial y\,\partial z}\exp \biggl( \int _{0}^{\upsilon}{\frac{\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr)\,d \upsilon} \\ &=\frac{D\tau_{0}}{1-\rho}>0, \end{aligned} \\& \begin{aligned} \frac{{\partial^{2}}\Phi_{2}(\tau_{0},S_{0})}{\partial z^{2}}&= \int _{0}^{\tau_{0}}{\exp \biggl( \int _{\upsilon}^{\tau_{0}}{\frac{\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr) \frac{{\partial^{2}}F_{2}(\varsigma (\upsilon))}{\partial z^{2}}\exp \biggl( \int _{0}^{\upsilon}{\frac {\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr)\,d \upsilon} \\ &\quad {}+ \int _{0}^{\tau_{0}}{ \biggl[\exp \biggl( \int _{\upsilon }^{\tau_{0}}{\frac{\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr) \frac {{\partial^{2}}F_{2}(\varsigma(\upsilon))}{\partial y\,\partial z} \biggr]} \\ &\quad {}\times \biggl[ \int _{0}^{\upsilon} \biggl(\exp \biggl( \int _{\theta}^{\upsilon}{\frac{\partial{{F}_{1}} ( \varsigma ( r ) )}{\partial y}\,dr} \biggr) \frac{\partial {{F}_{1}} ( \varsigma (\theta ) )}{\partial z}\\ &\quad {}\times\exp \biggl( \int _{0}^{\theta}{\frac{\partial{{F}_{2}} (\varsigma ( r ) )}{\partial z}\,dr} \biggr) \biggr)\,d\theta \biggr]\,d\upsilon \\ &= \int _{0}^{{{\tau}_{0}}}{\frac{2{c_{2}}f'(0)({k_{2}}-f^{2}(0))\tilde {z}(\upsilon)}{(1-\rho)(k_{2}+f^{2}(0))^{2}}\,d\upsilon} \\ &\quad {}+D \int _{0}^{{\tau}_{0}}{ \biggl[\exp \biggl( \int _{\upsilon }^{{{\tau}_{0}}}{ \bigl(D\tilde{z}(r)-d_{2} \bigr)\,dr} \biggr) \biggr]} \\ &\quad {}\times \biggl[ \int _{0}^{\upsilon} \biggl({\exp \bigl(-A(\upsilon - \theta) \bigr) \biggl(\frac{c_{3}}{k_{3}}-\frac{{c_{4}}{k_{4}}f'(0)\tilde {z}(\theta)}{(k_{4}+f(0))^{2}} \biggr)} \\ &\quad {}\times\exp \biggl( \int _{0}^{\theta}{ \bigl(D\tilde{z}(r)-d_{2} \bigr)\,dr} \biggr) \biggr)\,d\theta \biggr]\,d\upsilon, \end{aligned} \\& \begin{aligned} \frac{\partial^{2}{\Phi_{2}}(\tau_{0},S_{0})}{\partial z\,\partial\tau}&=\frac {\partial F_{2}(\varsigma(\tau_{0}))}{\partial z}\exp \biggl( \int _{0}^{\tau_{0}}{\frac{\partial F_{2}(\varsigma(r))}{\partial z}\,dr} \biggr) \\ &= \bigl(D\tilde{z}(\tau_{0})-d_{2} \bigr)\exp \biggl( \int _{0}^{\tau_{0}}{ \bigl(D\tilde {z}(r)-d_{2} \bigr)\,dr} \biggr) \\ &=\frac{1}{1-\rho} \biggl[\frac{D\mu\exp(-A\tau_{0})}{1-\exp(-A\tau _{0})}-d_{2} \biggr]. \end{aligned} \end{aligned}$$
Here, we can compute
$$\begin{aligned} d'_{0}=1-{ \biggl(\frac{\partial\Theta_{2}}{\partial z}\frac{\partial\Phi _{2}}{\partial z} \biggr)}_{(\tau_{0},S_{0})}=1-(1-\rho)\exp \biggl( \int _{0}^{\tau_{0}}{ \bigl(D\tilde{z}(r)-d_{2} \bigr)\,dr} \biggr), \end{aligned}$$
where \(\tau_{0}\) is the root of \(d'_{0}=0\). Note that \(d'_{0}>0\) if \(T>T_{\mathrm{max}} \) and \(d'_{0}<0\) if \(T< T_{\mathrm{max}} \).
$$\begin{aligned}& a'_{0} = 1-{ \biggl( \frac{\partial\Theta_{1}}{\partial y} \frac{\partial\Phi _{1}}{\partial y} \biggr)}_{(\tau_{0},S_{0})}=1-\exp(-A{\tau_{0}})>0, \\& \begin{aligned} b'_{0} &= -{ \biggl(\frac{\partial\Theta_{1}}{\partial y} \frac{\partial\Phi _{1}}{\partial z}+\frac{\partial\Theta_{1}}{\partial z}\frac{\partial\Phi _{2}}{\partial z} \biggr)}_{(\tau_{0},S_{0})} \\ & = -\frac{\partial\Phi_{1}(\tau_{0},S_{0})}{\partial z} \\ &= - \int _{0}^{\tau_{0}}{\exp \bigl(-A(\tau_{0}- \upsilon) \bigr) \biggl( \frac {c_{3}}{k_{3}}-\frac{{c_{4}}{k_{4}}f'(0)\tilde{z}(\upsilon)}{(k_{4}+f(0))^{2}} \biggr)}\exp \biggl( \int _{0}^{\upsilon}{ \bigl(D\tilde{z}(r)-d_{2} \bigr)\,dr} \biggr)\,d\upsilon. \end{aligned} \end{aligned}$$
Here, \(b'_{0}<0\) provided that
$$ c_{1}^{2}>d_{1}^{2} k_{1}^{2} k_{2}. $$
Since \((1-\frac{A{\tau_{0}}}{1-\exp(-A{\tau_{0}})} )\) is always negative for \({\tau_{0}}>0\), then
$$\begin{aligned}& \begin{aligned} C^{*}&=2\frac{\partial\Theta_{2}}{\partial z}\frac{b'_{0}}{a'_{0}}\frac {\partial^{2}\Phi_{2}}{\partial y\,\partial z}- \frac{\partial\Theta _{2}}{\partial z}\frac{\partial^{2}\Phi_{2}}{\partial z^{2}} \\ &=2( 1-\rho)\frac{b'_{0}}{a'_{0}}\frac{\partial^{2}\Phi_{2}}{\partial y\,\partial z}-(1-\rho)\frac{\partial^{2}\Phi_{2}}{\partial z^{2}}< 0, \end{aligned} \\& \begin{aligned} B^{*}&=-\frac{\partial\Theta_{2}}{\partial z} \biggl(\frac{\partial^{2}\Phi _{2}}{\partial y\,\partial z} \frac{1}{a'_{0}}\frac{\partial\Theta _{1}}{\partial y}\frac{\partial\Phi_{1}}{\partial\tau}+\frac{\partial ^{2}\Phi_{2}}{\partial\tau\partial z} \biggr) \\ &=- \biggl[\frac{D\mu\exp(-A{\tau_{0}})}{1-\exp(-A{\tau_{0}})} \biggl( 1-\frac {A\tau_{0}}{1-\exp(-A{\tau_{0}})} \biggr)-d_{2} \biggr]>0, \end{aligned} \end{aligned}$$
and hence, \(B^{*}C^{*}<0\). According to Lakmeche and Arini [18], the following result is obtained.
System (27)-(30) has a positive periodic solution provided (5), (11), (31) hold, and \(T< T_{\mathrm{max}} \).
In the next section, numerical simulations are given in order to confirm our theoretical predictions.
Numerical results
To support our prediction in Theorem 4.1, we choose all parameters so that all the conditions in Theorem 4.1 are satisfied. The simulation result is as shown in Figure 1. We can see that the solution converges asymptotically to the oscillatory solution \((0,\tilde{z}(t))\) as predicted. The simulation result in Figure 2 is obtained from choosing all parameters to satisfy all the conditions in Theorem 5.1. In this case, we can see that the solution is bounded within a positive range as we expected as well. Finally, we also carry out numerical simulation to confirm our prediction in Theorem 6.1. The result is as shown in Figure 3 in which the solution trajectory tends toward the sustained oscillations when all parameters are chosen to satisfy all the conditions in Theorem 6.1.
Numerical simulation of equations (3a)-(3d). The solution trajectory approaches oscillatory solution \((0,\tilde {z}(t))\) as time passes. Here, all parameters are chosen to satisfy the conditions in Theorem 4.1, i.e., \({c_{1}}=0.15\), \({c_{2}}=0.8\), \({c_{3}}=0.35\), \({c_{4}}=0.9\), \({d_{1}}=0.9\), \({d_{2}}=0.5\), \({d_{3}}=0.1\), \({k_{1}}=1.1\), \({k_{2}}=0.9\), \({k_{3}}=0.9\), \({k_{4}}=3.9\), \(\mu=0.8\), \(\rho =0.5\), \(T=20\), \(y(0)=0.13\), and \(z(0)=5\). (a) The solution trajectory projected on \((y,z)\)-plane. (b) The corresponding time course of the number of active osteoclasts \((y)\) tending towards zero. (c) The corresponding time course of the number of active osteoblasts \((z)\) exhibiting positive oscillation.
Numerical simulation of equations (3a)-(3d). The solution trajectory is bounded within a positive range as time passes. Here, all parameters are chosen to satisfy the conditions in Theorem 5.1. Here, \({c_{1}}=0.5\), \({c_{2}}=0.9\), \({c_{3}}=0.35\), \({c_{4}}=0.9\), \({d_{1}}=0.95\), \({d_{2}}=0.05\), \({d_{3}}=0.5\), \({k_{1}}=1.2\), \({k_{2}}=0.95\), \({k_{3}}=0.9\), \({k_{4}}=3.9\), \(\mu=0.9\), \(\rho=0.1\), \(T=5\), \(y(0)=0.1\), and \(z(0)=2\). \((a)\) The solution trajectory projected on \((y,z)\)-plane. \((b)\) The corresponding time course of the number of active osteoclasts \((y)\) and \((c)\) the corresponding time course of the number of active osteoblasts \((z)\).
Numerical simulation of equations (3a)-(3d). The solution trajectory approaches a limit cycle as time passes. Here, all parameters are chosen to satisfy the conditions in Theorem 6.1. Here, \({c_{1}}=0.7\), \({c_{2}}=0.2\), \({c_{3}}=0.35\), \({c_{4}}=0.9\), \({d_{1}}=0.95\), \({d_{2}}=0.005\), \({d_{3}}=0.7\), \({k_{1}}=1.2\), \({k_{2}}=0.35\), \({k_{3}}=0.9\), \({k_{4}}=3.9\), \(\mu=0.9\), \(\rho=0.1\), \(T=10\), \(y(0)=0.1\), and \(z(0)=2\). (a) The solution trajectory projected on \((y,z)\)-plane. (b) The corresponding time course of the number of active osteoclasts \((y)\) and (c) the corresponding time course of the number of active osteoblasts \((z)\) exhibiting sustained oscillation.
The impulsive mathematical model is proposed to incorporate the effects of PRL and PTH supplement on bone remodeling process. The conditions on system parameters are then derived so that we can expect desirable behaviors of the solution of the system such as the solution will be bounded within a positive range, which means that the number of bone forming cells and bone resorbing cells is able to be controlled and lies within the appropriate levels.
From Theorem 5.1, one of the conditions used to guarantee the boundedness of the number of active osteoblasts and active osteoclasts is \(T< T_{\mathrm{max}} \), where \(T_{\mathrm{max}} \) depends on μ and ρ. We can see that, if the dosage of PTH supplement reflected by μ and ρ is fixed, we can choose an appropriate frequency of PTH supplement \(\frac {1}{T}\) for which \(T< T_{\mathrm{max}} \) as required by Theorem 5.1 to guarantee the appropriate levels of active osteoblasts and active osteoclasts. Moreover, Theorem 6.1 implies that if all the conditions required in Theorem 5.1 are satisfied and PRL is secreted from the anterior pituitary gland with the sufficient amount, i.e., \(c_{1}^{2}>d_{1}^{2} k_{1}^{2} k_{2}\), the oscillatory behavior resembling clinical observation can be expected. Hence, the frequency and dosage of PTH supplement are the keys to the effectiveness of PTH supplement in osteoporosis patients.
Goodman, HM: Basic Medical Endocrinology, 3rd edn. Academic Press, San Diego (2003)
Albright, JA, Sauders, M: The Scientific Basis of Orthopaedics. Appleton and Lange, Norwalk (1990)
Marcus, R: Osteoporosis. Blackwell Scientific Publication, Oxford (1994)
Lacroix, PC, Ormandy, C, Lepescheux, L, Ammann, P, Damotte, D, Goffin, V, Bouchard, B, Amling, M, Kelly, MG, Binart, N, Baron, R, Kelly, PA: Osteoblasts are a new target for prolactin: analysis of bone formation in prolactin. Endocrinology 140, 96-105 (1999)
Charoenphandhu, N, Tudpor, K, Thongchote, K, Saengamnart, W, Puntheeranurak, S, Krishnamara, N: High-calcium diet modulates effects of long-term prolactin exposure on the cortical bone calcium content in ovariectomized rats. Am. J. Physiol: Endocrinol. Metab. 292, E443-E452 (2007)
Momsen, G, Schwarz, P: A mathematical/physiological model of parathyroid hormone secretion in response to blood-ionized calcium lowering in vivo. Scand. J. Clin. Lab. Invest. 57, 381-394 (1997)
Kroll, MH: Parathyroid hormone temporal effects on bone formation and resorption. Bull. Math. Biol. 62, 163-188 (2000)
Tam, CS, Heersche, JNM, Murray, TM, Parsons, JA: Parathyroid hormone stimulates the bone apposition rate independently of its resorptive action: differential effects of intermittent and continuous administration. Endocrinology 110(2), 506-512 (1982)
Hock, JM, Gera, I: Effects of continuous and intermittent administration and inhibition of resorption on the anabolic response of bone to parathyroid hormone. J. Bone Miner. Res. 7(1), 65-72 (1992)
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. Comput. Simul. 6(5), 510-519 (2011)
Rattanakul, C, Lenbury, Y, Krishnamara, N, Wollkind, DJ: Mathematical modelling of bone formation and resorption mediated by parathyroid hormone: responses to estrogen/PTH therapy. Biosystems 70, 55-72 (2003)
Rattanamongkonkul, S, Kunpasuruang, W, Ruktamatakul, S, Rattanakul, C: A mathematical model of bone remodeling process: effect of vitamin D. Int. J. Math. Comput. Simul. 6(5), 489-498 (2011)
Rattanakul, C, Rattanamongkonkul, S: Effect of calcitonin on bone formation and resorption: mathematical modeling approach. Int. J. Math. Models Methods Appl. Sci. 5(8), 1363-1371 (2011)
Rattanakul, C, Rattanamongkonkul, S, Kunpasuruang, W, Ruktamatakul, S, Srisuk, S: A mathematical model of bone remodeling process: effects of parathyroid hormone and vitamin D. Int. J. Math. Models Methods Appl. Sci. 5(8), 1388-1397 (2011)
Coss, D, Yang, L, Kuo, CB, Xu, X, Luben, RA, Walker, AM: Effects of prolactin on osteoblast alkaline phosphatase and bone formation in the developing rat. Am. J. Physiol: Endocrinol. Metab. 279, E1216-E1225 (2000)
Liu, B, Zhi, Y, Chen, LS: The dynamics of a predator-prey model with Ivlev's functional response concerning integrated pest management. Acta Math. Appl. Sin. 20(1), 133-146 (2004)
Lakshmikantham, V, Bainov, DD, Simeonov, PS: Theory of Impulsive Differential Equations. World Scientific, Singapore (1989)
Lakmeche, A, Arini, O: Bifurcation of nontrivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment. Dyn. Contin. Discrete Impuls. Syst. 7, 265-287 (2000)
This work was supported by the Royal Golden Jubilee Ph.D. Program (contract number PHD53K0191), the Centre of Excellence in Mathematics, Commission on Higher Education, Thailand and Mahidol University, Thailand.
Department of Mathematics, Faculty of Science, Mahidol University, 272 Rama VI Road, Bangkok, 10400, Thailand
Inthira Chaiya
& Chontita Rattanakul
Centre of Excellence in Mathematics, Commission on Higher Education, 328 Sri Ayuthya Road, Bangkok, 10400, Thailand
Chontita Rattanakul
Search for Inthira Chaiya in:
Search for Chontita Rattanakul in:
Correspondence to Chontita Rattanakul.
The first author analyzed the model theoretically. The second author developed the model and carried out numerical simulations. All authors read and approved the final manuscript.
Chaiya, I., Rattanakul, C. Effects of prolactin on bone remodeling process with parathyroid hormone supplement: an impulsive mathematical model. Adv Differ Equ 2017, 147 (2017) doi:10.1186/s13662-017-1199-x
DOI: https://doi.org/10.1186/s13662-017-1199-x
bone remodeling process
mathematical model
parathyroid hormone supplement
Proceedings of the International Conference in Mathematics and Application 2016