Existence and stability analysis of solutions for a ultradian glucocorticoid rhythmicity and acute stress model

TL;DR

This paper analyzes the existence and stability of solutions in mathematical models of the HPA axis, focusing on ultradian glucocorticoid rhythms and stress response, including effects of system delays.

Contribution

It introduces a mathematical framework for analyzing stability and periodic solutions in HPA axis models, incorporating delays and nonlinear dynamics.

Findings

Existence of periodic solutions under specific parameter conditions

Stability of solutions depends on the delay parameter

Linear and nonlinear stability analyses conducted

Abstract

The hypothalamic pituitary adrenal (HPA) axis responds to physical and mental challenge to maintain homeostasis in part by controlling the body's cortisol level. Dysregulation of the HPA axis is implicated in numerous stress-related diseases. For a structured model of the HPA axis that includes the glucocorticoid receptor but does not take into account the system response delay, we analyze linear and non-linear stability of stationary solutions. For a second mathematical model that describes the mechanism of the HPA axis self-regulatory activities and takes into account a delay of system response, we prove existence of periodic solutions under certain assumptions on ranges of parameter values and analyze stability of these solutions with respect to the time delay value.

Figures

Equations

$$\frac{da}{dt}=\frac{CRH}{1+p_{2}or}-p_{3}a=:f_1,$$

$$\frac{dr}{dt}=\frac{(or{)}^2}{p_4+(or{)}^2}+p_5-p_{6}r=:f_2,$$

$$\frac{do}{dt}=a(t-\tau )-o=:f_3$$

$$a(t)=a_{\tau }(t)\forall t\in [-\tau ,0],r(0)=r_0,o(0)=o_0,$$

$$0⩽a_{\tau }(t)\in C^1[-\tau ,0],r_0>0,o_0>0.$$

$$\frac{da}{dt}=\frac{A}{1+p_{2}or}-p_{3}a,$$

$$\frac{dr}{dt}=\frac{(or{)}^2}{p_4+(or{)}^2}+p_5-p_{6}r,$$

$$\frac{do}{dt}=a-o,$$

$$a(0)=a_0⩾0,r(0)=r_0>0,o(0)=o_0>0,$$

$$o=a,\frac{A}{1+p_{2}or}-p_{3}a=0,\frac{(or{)}^2}{p_4+(or{)}^2}+p_5-p_{6}r=0.$$

$$D:=\{(a,r,o)\in R_{+}^3:a=o,0⩽a⩽\frac{A}{p_3},\frac{p_5}{p_6}⩽r⩽\frac{p_5+1}{p_6}\}.$$

$$o=a=\tfrac{1}{2p_{2}r}\biggl{[}\sqrt{1+\tfrac{4p_{2}A}{p_{3}}\,r}-1\biggr{]},$$

$$\tfrac{1}{4p_{2}^{2}}\biggl{[}\sqrt{1+\tfrac{4p_{2}A}{p_{3}}\,r}-1\biggr{]}^{2}=\tfrac{p_{4}(p_{6}r-p_{5})}{1+p_{5}-p_{6}r},\quad\mbox{i.\,e.}\quad A=\tfrac{p_{3}}{r}\sqrt{\tfrac{p_{4}(p_{6}r-p_{5})}{1+p_{5}-p_{6}r}}\Bigl{(}1+p_{2}\sqrt{\tfrac{p_{4}(p_{6}r-p_{5})}{1+p_{5}-p_{6}r}}\Bigr{)}.$$

$$z:=\sqrt{1+\frac{4p_{2}A}{p_3}r}-1⩾0\Rightarrow r=\frac{p_3}{4p_{2}A}z(z+2).$$

$$z^4+2z^3+C_1{z}^2+C_{2}z-C_3=0,$$

$$C_1:=\frac{4p_2(p_2{p}_3{p}_4{p}_6-A(p_5+1))}{p_3{p}_6},C_2:=8p_2^2{p}_4>0,C_3:=\frac{16A{p}_2^3{p}_4{p}_5}{p_3{p}_6}>0.$$

$$a^{*}=o^{*}=\tfrac{1}{2p_{2}r^{*}}\biggl{[}\sqrt{1+\tfrac{4p_{2}A}{p_{3}}\,r^{*}}-1\biggr{]}=\tfrac{1}{r^{*}}\sqrt{\tfrac{p_{4}(p_{6}r^{*}-p_{5})}{1+p_{5}-p_{6}r^{*}}},$$

$$J^{*}=\left(\begin{array}[]{ccc}-p_{3}&-K_{1}&-K_{3}\\ 0&-p_{6}+K_{2}&K_{4}\\ 1&0&-1\\ \end{array}\right),$$

$$K_{1}=\tfrac{Ap_{2}a^{*}}{(1+p_{2}a^{*}r^{*})^{2}}=\tfrac{2A(\sqrt{1+\frac{4p_{2}A}{p_{3}}\,r^{*}}-1)}{r^{*}(1+\sqrt{1+\frac{4p_{2}A}{p_{3}}\,r^{*}})^{2}}=\tfrac{Ap_{2}\sqrt{p_{4}(p_{6}r^{*}-p_{5})(1+p_{5}-p_{6}r^{*})}}{r^{*}(p_{2}\sqrt{p_{4}(p_{6}r^{*}-p_{5})}+\sqrt{1+p_{5}-p_{6}r^{*}})^{2}}=\\ \tfrac{p_{2}p_{3}p_{4}(p_{6}r^{*}-p_{5})}{(r^{*})^{2}(p_{2}\sqrt{p_{4}(p_{6}r^{*}-p_{5})}+\sqrt{1+p_{5}-p_{6}r^{*}})^{2}}\Bigl{(}1+p_{2}\sqrt{\tfrac{p_{4}(p_{6}r^{*}-p_{5})}{1+p_{5}-p_{6}r^{*}}}\Bigr{)}\geqslant 0,$$

$$K_{2}=\tfrac{2p_{4}r^{*}(a^{*})^{2}}{(p_{4}+(a^{*}r^{*})^{2})^{2}}=\tfrac{1}{2p_{2}^{2}p_{4}r^{*}}\bigl{[}\sqrt{1+\tfrac{4p_{2}A}{p_{3}}\,r^{*}}-1\bigr{]}^{2}(1+p_{5}-p_{6}r^{*})^{2}=\\ \tfrac{2}{r^{*}}(p_{6}r^{*}-p_{5})(1+p_{5}-p_{6}r^{*})\geqslant 0\text{ and }\ 0\leqslant K_{2}\leqslant 2p_{6}(\sqrt{1+p_{5}}-\sqrt{p_{5}})^{2},$$

$$K_3=\frac{r^{\ast }}{a^{\ast }}{K}_1=\frac{p_2{p}_3\sqrt{p_4(p_6{r}^{\ast }-p_5)}}{\sqrt{1+p_5-p_6{r}^{\ast }}+p_2\sqrt{p_4(p_6{r}^{\ast }-p_5)}}\text{ and }0⩽K_3⩽p_3,$$

$$K_4=\frac{r^{\ast }}{a^{\ast }}{K}_2=\frac{2r^{\ast }}{\sqrt{p_4}}(1+p_5-p_6{r}^{\ast }{)}^{\frac{3}{2}}\sqrt{(p_6{r}^{\ast }-p_5)}⩾0.$$

$$|J^{*}-\lambda I|=\left|\begin{array}[]{ccc}-p_{3}-\lambda&-K_{1}&-K_{3}\\ 0&-p_{6}+K_{2}-\lambda&K_{4}\\ 1&0&-1-\lambda\\ \end{array}\right|=0,$$

$$(\lambda +1)(\lambda +p_3)(\lambda +p_6-K_2)+K_3(\lambda +p_6)=0,$$

$${\lambda }^3+{\alpha }_1{\lambda }^2+{\alpha }_2\lambda +{\alpha }_3=0,$$

$${\alpha }_1=p_3+p_6-K_2+1,{\alpha }_2=p_3+p_6-K_2+p_3(p_6-K_2)+K_3,{\alpha }_3=p_3(p_6-K_2)+p_6{K}_3.$$

$$\Delta :=18{\alpha }_1{\alpha }_2{\alpha }_3-4{\alpha }_1^3{\alpha }_3+{\alpha }_1^2{\alpha }_2^2-4{\alpha }_2^3-27{\alpha }_3^2.$$

$${\alpha }_1>0\iff 0⩽K_2<p_6+p_3+1\iff 2p_6^2{x}^2+(p_3+1-p_6(1+4p_5))x+2p_5(p_5+1)>0$$

$${\alpha }_3>0\iff 0⩽K_2<p_6(1+\frac{K_3}{p_3})\iff (x-x_1)(x-x_2)<\frac{x}{2p_6}[1+\frac{p_2\sqrt{p_4(x-x_1)}}{\sqrt{x_2-x}+p_2\sqrt{p_4(x-x_1)}}]$$

$$\displaystyle\alpha_{1}\alpha_{2}>\alpha_{3}\Leftrightarrow 0\leqslant K_{2}<p_{6}+\tfrac{1}{2}\biggl{[}p_{3}+1+\tfrac{K_{3}}{1+p_{3}}$$

Full text

1 Introduction

Hormones control a vast array of bodily functions, including sexual reproduction and sexual development, whole-body metabolism, blood glucose levels and so on (see [22]). Hormones are produced, in main, and released from diverse places including the hypothalamus, pituitary, and the adrenal gland. Hormones are capable of a diffusion whole-body effect, as well as a localized effect, depending on the distance between the production site and the site of action. In many ways, the endocrine system is similar to the nervous system, in that it is an intercellular signaling system in which cells communicate via cellular secretions. Further, the distance between the sites of hormone production and action, and the complexities inherent in the mode of transport, make it extraordinarily difficult to construct quantitative models of hormonal control.

The hypothalamus pituitary adrenal axis is a central neuroendocrine system, which consists of the hypothalamus, pituitary, and adrenal glands. The paraventricular nucleus of the hypothalamus secrets corticotropin releasing hormone (CRH), which is transferred to the pituitary and stimulates the synthesis and release of adrenocorticotropic hormone (ACTH). ACTH moves through the bloodstream and reaches the adrenal gland in which it stimulates the secretion of cortisol. In response to stress, the concentrations of the HPA axis hormones are increased.

Disruption of HPA axis regulation is known to contribute to a number of stress-related disorders. For example, increased cortisol has been shown in patients with major depressive disorder (see [11, 6]), and decreased cortisol has been observed in people with post-traumatic stress disorder (see [18]).

Multiple models of the hypothalamus-pituitary-adrenal (HPA) axis have been developed to characterize the oscillations seen in the hormone concentrations and to examine HPA axis dysfunction. Most of these models have been constructed using deterministic coupled ordinary differential equations (see [7]). A major inconsistency among different existing HPA models that was mentioned in [9] is related to their treatment of the circadian and ultradian oscillations. For example, the authors of [20] and [2] assumed that both oscillations can be generated inside the HPA axis system by interaction of its elements; the authors of [21], [1], and [10] treat the circadian and ultradian oscillations differently assuming that only ultradian oscillations are HPA axis based but at the same time circadian rhythms are due to external input. Only one model made no explicit assumption about the origin of the oscillations and was developed to replicate the HPA axis response to CRH injection (see [4]). It has also been suggested that the ultradian rhythm arises from the introduction of a time delay (see [2]). Other models based on delay-differential equations include [14] and [8].

To determine if delay-differential equations could predict the general features of cortisol production, the experimental data was compared to a simulated cortisol curve in [8]. Experimental fitting of ACTH for the model was not possible since hypothalamic derived CRH cannot be measured.

Inclusion of the glucocorticoid receptor in a hypothalamic pituitary adrenal axis model reveals ’bi-stability’ (see [8]). To be more concrete, there arises a nonlinear Gauss type function with compact support, which is characterized by the parameter $p_{4}$. This (Hill function) arises as a result of ’inner’ nonlinearity in the physiological system which is produced by the stress impulse, which is activated by the outer impulse that is called by an acute stress. This situation is provided formally by the two parameters $p_{4}$ and $CRH$. The amplitude of the Hill function determines glucocorticoid receptor $GR$ density in the pituitary, which is coupled nonlinearly in reaction with regulated levels of $CORT$, which in turn mediate a wide range of physiological processes, including metabolic, immunological and cognitive function (see [17, 16]).

The stress response is subserved by the stress system which is located both in the central nervous system and the periphery. The principal effects to the stress system include the corticotropin-releasing hormone $CRH$. The secretion of $CRH$ causes the anterior pituitary to synthesize adrenocorticotropin $ACTH$ which then stimulates the adrenal glands to release cortisol that regulate the blood concentration of $CRH$ and $ACTH$ via different negative feedback mechanisms. The $HPA$ axis is the subject of intensive research in endocrinology. This model is based on the feed-forward and feedback interactions between the anterior pituitary and adrenal glands. Because responsiveness of the stress system to stressors is crucial for life, it is important to consider the simpler case when distributions of hormones in the system become unstable by action on stress, and further to consider influence on the delay time as response of the physiological system on action on stress.

Mathematically, it means that we can consider two mathematical models: the first one is described by a system of ordinary differential equations with initial distributions of hormones at a point $t=0$, and the second one is based upon a system of differential equations with initial distributions of hormones on the interval $[-\tau,0)$, where $\tau$ is a time delay. It turns out that bi-stability is present in both models, i.e limit distributions of hormones may be stable or unstable depending on parameter values. In the model with a time delay, periodic solutions arise for special distributions of hormones when there is a connection between the concentrations of hormones at end points. Thus, the initial distributions of hormones must be coupled in a special manner. This condition may be considered as a ’normal’ reaction of the organism on action of stress.

In this paper, we study a system of delay differential equations (see [8, 17]):

[TABLE]

[TABLE]

[TABLE]

with initial conditions

[TABLE]

where

[TABLE]

Based on the principles of mass action kinetics, these equations describe the production and degradation of the hormones $ACTH$ ($a$), i. e. adrenocorticotropin, and $CORT$ ($o$), i. e. cortisol in humans and corticosterone in rodents, as well as glucocorticoid receptor $GR$ density ($r$) in the pituitary. Here, $CRH$ is corticotrophin-releasing hormone, the parameters $p_{2-6}$ represent dimensionless forms of rate constants of the system, and the dimensionless parameter $\tau$ represents a discrete delay, which accounts for the delayed response of the adrenal gland to $ACTH$. The dimensionless time $t=0$ corresponds to the maximal value of an $ACTH$ pulse.

The paper is organized as follows. In Section 2 we study the stability of the system (1.1)–(1.3) without delay for the given initial distributions of hormones. In Section 3 we analyze the system (1.1)–(1.3) with delay in terms of stability, solvability, and existence of periodic solutions. Moreover, the existence of $\tau$-periodic solutions is proved. The system without delay is always asymptotically stable for strictly positive parameter values. When delay is considered, the system remains unstable after a certain delay value. The theoretical results are then compared with numerical simulations.

2 Stability analysis of the model without time delay

We consider the following nonlinear ODEs without delay:

[TABLE]

[TABLE]

[TABLE]

with initial conditions

[TABLE]

where $A:=CRH\geqslant 0$, and $p_{i}\geqslant 0$. Using Picard’s iteration method, we can show existence of a unique global in time non-negative solution. The proof is similar to the one of Theorem 1.

Lemma 2.1**.**

Assume that $A>0$ and $p_{i}>0$. Then the system (2.1)–(2.3) has a unique fixed point and this point is asymptotically stable.

Proof of Lemma 2.1.

By (2.1)–(2.3) we obtain the following equations for the nullclines:

[TABLE]

The algebraic system (2.5) has a nonnegative solution in the following domain:

[TABLE]

From (2.5) we have

[TABLE]

[TABLE]

Note that the equation (2.8) has a unique solution $r^{*}\in(\frac{p_{5}}{p_{6}},\frac{p_{5}+1}{p_{6}})$ in $D$. Really, the function f_{1}(r):=\frac{1}{4p_{2}^{2}}\biggl{[}\sqrt{1+\frac{4p_{2}A}{p_{3}}\,r}-1\biggr{]}^{2} is nonnegative and monotone increasing on $[0,+\infty)$ such that $f_{1}(0)=0$, $f_{1}(+\infty)=+\infty$, but the function $f_{2}(r):=\frac{p_{4}(p_{6}r-p_{5})}{1+p_{5}-p_{6}r}$ is nonnegative and monotone increasing on $[\frac{p_{5}}{p_{6}},\frac{p_{5}+1}{p_{6}}]$ such that $f^{\prime\prime}_{2}(r)>0$, $f_{2}(\frac{p_{5}}{p_{6}})=0$, $f_{2}(\frac{p_{5}+1}{p_{6}})=+\infty$. Therefore, there is only one intersection of $f_{1}(r)$ and $f_{2}(r)$ on the interval $[\frac{p_{5}}{p_{6}},\frac{p_{5}+1}{p_{6}}]$. On the other hand, let us denote by

[TABLE]

Then (2.8) can be rewritten in the following form

[TABLE]

where

[TABLE]

So, we can find the explicit value of $r^{*}\in\left(\frac{p_{5}}{p_{6}},\frac{p_{5}+1}{p_{6}}\right)$ as a solution of (2.9).

As a result, the system (2.1)–(2.3) has only one fixed point $(a^{*},r^{*},o^{*})$ in $D$. Here,

[TABLE]

and $r^{*}$ is the solution of (2.8) or (2.9).

Next, we find the Jacobian matrix $J^{*}$ for (2.1)–(2.3) at the fixed point $(a^{*},r^{*},o^{*})$.

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Next, we will analyze the stability of the fixed point. First, we look for eigenvalues for $J^{*}$. So,

[TABLE]

whence we obtain the characteristic equation:

[TABLE]

i. e.

[TABLE]

where

[TABLE]

Let us denote by

[TABLE]

If $\Delta>0$, then (2.10) has three distinct real roots. If $\Delta=0$, then (2.10) has a multiple root and all of its roots are real. If $\Delta<0$, then (2.10) has one real root and two complex roots.

To analyze stability, we will use Lemma 5.1 in the Appendix. Let $x_{1}:=\frac{p_{5}}{p_{6}}\leqslant x:=r^{*}\leqslant x_{2}:=\frac{p_{5}+1}{p_{6}}$. Then in our case,

[TABLE]

which is true for all $p_{i}>0$;

[TABLE]

which is true for all $p_{i}>0$;

[TABLE]

provided

[TABLE]

but if \biggl{(}\frac{K_{3}}{1+p_{3}}\biggr{)}^{2}-2\biggl{[}1-\frac{2p_{6}}{1+p_{3}}\biggr{]}K_{3}+(p_{3}-1)^{2}\leqslant 0 then these are true for all $p_{i}>0$. Thus the fixed point is asymptotically stable. ∎

2.1 Stability for different parameter values

Case 1: If $A=0$ then the system (2.1)–(2.3) has the fixed point $\left(0,\frac{p_{5}}{p_{6}},0\right)$. The corresponding characteristic equation is

[TABLE]

whence $\lambda_{i}=-1,\,-p_{3},\,-p_{6}<0$. As a result, $\left(0,\frac{p_{5}}{p_{6}},0\right)$ is stable node.

Case 2: If $p_{2}=0$ then the system (2.1)–(2.3) has the fixed point $\left(\frac{A}{p_{3}},r^{*},\frac{A}{p_{3}}\right)$, where $r^{*}$ is a solution of the equation:

[TABLE]

This equation can also be written as

[TABLE]

which yields the following cases:

• •

if $p_{4}>\left[\frac{A(1+p_{5})}{p_{3}p_{6}}\right]^{2}$ then there is one real root;

• •

if $p_{4}=\left[\frac{A(1+p_{5})}{p_{3}p_{6}}\right]^{2}$ then $r\left(r-\frac{1+p_{5}}{2p_{6}}\right)^{2}=\frac{p_{4}p_{5}}{p_{6}}$, whence

• if $\frac{p_{5}A^{2}}{p_{3}^{2}p_{6}}<\frac{1}{54}$ then we have three real roots,
• if $\frac{p_{5}A^{2}}{p_{3}^{2}p_{6}}=\frac{1}{54}$ then we have two real roots,
• if $\frac{p_{5}A^{2}}{p_{3}^{2}p_{6}}>\frac{1}{54}$ then we have one real root;
• •

if $p_{4}<\left[\frac{A(1+p_{5})}{p_{3}p_{6}}\right]^{2}$ then $r(r-r_{1})(r-r_{2})=\frac{p_{4}p_{5}}{p_{6}}\left(\frac{p_{3}}{A}\right)^{2}$, where

$r_{1,2}=\frac{1}{2p_{6}}\left[1+p_{5}\pm\sqrt{(p_{5}+1)^{2}-p_{4}p_{6}^{2}\left(\frac{p_{3}}{A}\right)^{2}}\right]$, whence

• if $\frac{p_{4}p_{5}}{p_{6}}(\frac{p_{3}}{A})^{2}<\frac{(r_{1}+r_{2}-K)(r_{2}-2r_{1}-K)(r_{1}-2r_{2}-K)}{27}$ then we have three real roots,
• if $\frac{p_{4}p_{5}}{p_{6}}(\frac{p_{3}}{A})^{2}=\frac{(r_{1}+r_{2}-K)(r_{2}-2r_{1}-K)(r_{1}-2r_{2}-K)}{27}$ then we have two real roots,
• if $\frac{p_{4}p_{5}}{p_{6}}(\frac{p_{3}}{A})^{2}>\frac{(r_{1}+r_{2}-K)(r_{2}-2r_{1}-K)(r_{1}-2r_{2}-K)}{27}$ then we have one real root,

where $K^{2}=r_{1}^{2}-r_{1}r_{2}+r_{2}^{2}$.

The corresponding characteristic equation is

[TABLE]

whence $\lambda_{i}=-1,\,-p_{3},\,-p_{6}+K_{2}$. If $K_{2}<p_{6}$ then $(\frac{A}{p_{3}},r^{*},\frac{A}{p_{3}})$ is stable node. If $K_{2}>p_{6}$ then $(\frac{A}{p_{3}},r^{*},\frac{A}{p_{3}})$ is saddle. If $K_{2}=p_{6}$ then it is a non-hyperbolic fixed point.

Case 3: If $p_{3}=0$ then the system (2.1)–(2.3) has the fixed point $(+\infty,\frac{p_{5}+1}{p_{6}},+\infty)$. The corresponding characteristic equation is

[TABLE]

whence $\lambda_{i}=-1,\,0,\,-p_{6}$. As a result, $(+\infty,\frac{p_{5}+1}{p_{6}},+\infty)$ is non-hyperbolic fixed point.

Case 4: If $p_{4}=0$ then the system (2.1)–(2.3) has the fixed point $\left(a^{*},\frac{p_{5}+1}{p_{6}},a^{*}\right)$, where

[TABLE]

In this case, we have that $K_{2}=K_{3}=0$ and $(\lambda+1)(\lambda+p_{3})(\lambda+p_{6})=0,$ whence $\lambda_{i}=-1,\,-p_{3},\,-p_{6}$. Hence, the fixed point is a stable node.

Case 5: If $p_{2}=p_{4}=0$ then we obtain the explicit solution

[TABLE]

Case 6: If $p_{5}=0$ then the one of fixed points is $\left(\frac{A}{p_{3}},0,\frac{A}{p_{3}}\right)$ and

[TABLE]

whence $\lambda_{i}=-1,\,-p_{3},\,-p_{6}$. Hence, this fixed point is stable node. In this case, by (2.5) we obtain that

[TABLE]

whence

[TABLE]

implies that

[TABLE]

Hence, we find that

[TABLE]

where

[TABLE]

and $a_{1}<0<a_{2}$. Note that $a_{2}\leqslant\frac{A}{p_{3}}$ provided $p_{2}^{2}p_{4}\geqslant 1$. As a result,

• •

if $p_{2}^{2}p_{4}\geqslant 1$ then

• if $f_{min}>-\frac{p_{6}}{p_{2}}(\frac{A}{p_{3}})^{2}$ then no real roots;
• if $f_{min}=-\frac{p_{6}}{p_{2}}(\frac{A}{p_{3}})^{2}$ then one positive real root;
• if $f_{min}<-\frac{p_{6}}{p_{2}}(\frac{A}{p_{3}})^{2}$ then two positive real roots;
• •

if $p_{2}^{2}p_{4}<1$ then

• if $f_{min}\geqslant-\frac{p_{6}}{p_{2}}(\frac{A}{p_{3}})^{2}$ then no real roots;
• if $f_{min}<-\frac{p_{6}}{p_{2}}(\frac{A}{p_{3}})^{2}$ then one positive real root.

Case 7: If $p_{6}=p_{5}=0$ then we have the following system

[TABLE]

If $r_{0}=0$ then we find the explicit solution

[TABLE]

If $r_{0}\neq 0$ then we approximately have

[TABLE]

whence

[TABLE]

If $p_{6}=0$ but $p_{5}\neq 0$ then $r(t)$ blows up in a finite time too.

Also, note that if $p_{6}=0$ and $p_{5}=0$ then the system (2.1)–(2.3) has the fixed point $(\frac{A}{p_{3}},0,\frac{A}{p_{3}})$. The corresponding characteristic equation is

[TABLE]

whence $\lambda_{i}=-1,\,-p_{3},\,0$. As a result, $(\frac{A}{p_{3}},0,\frac{A}{p_{3}})$ is non-hyperbolic fixed point.

Example 2.1**.**

Let $A=1$, $p_{2}=15$, $p_{3}=7.2$, $p_{4}=0.05$, $p_{5}=0.11$, and $p_{6}=2.9$. Then $r^{*}\approx 0.03$, $a^{*}=o^{*}\approx 0.12$, $\alpha_{1}=11.1-K_{2}\approx 11.07$, $\alpha_{2}=30.98-8.2K_{2}+K_{3}\approx 30.78$, $\alpha_{3}=20.88-7.2K_{2}+2.9K_{3}\approx 20.75$, whence we find that $\Delta\approx 2509.05>0$, $\alpha_{1}\alpha_{2}>\alpha_{3}$. As a result, all characteristic roots are negative real numbers and the fixed point is stable node.

A visual representation of this stable node can be found in Figure 1. This plot was created using the Matlab ode45 solver ([15]) using various starting values and the parameter values given above. The starting values were selected so that $a_{0}>0$, $r_{0}>0$ and $o_{0}>0$ to imitate real initial hormone levels.

Example 2.2**.**

Let $A=0.106$, $p_{2}=0$, $p_{3}=0.222$, $p_{4}=0.464$, $p_{5}=0.094$, and $p_{6}=0.418$. Then $r^{*}\approx 0.39,0.83,1.38$ and $a^{*}=o^{*}\approx 0.47$. Using similar calculations as above according to the defined values. If $r^{*}\approx 0.39$ then $\alpha_{1}\approx 1.30$, $\alpha_{2}\approx 0.32$, $\alpha_{3}\approx 0.01$, and $K_{2}\approx 0.33<p_{6}$ which means this is a stable node. If $r^{*}\approx 0.83$ then $\alpha_{1}\approx 1.18$, $\alpha_{2}\approx 0.17$, and $\alpha_{3}\approx-0.008$, and $K_{2}\approx 0.45>p_{6}$ which means this is a saddle. If $r^{*}\approx 1.38$ then $\alpha_{1}\approx 1.28$, $\alpha_{2}\approx 0.29$, $\alpha_{3}\approx 0.01$, and $K_{2}\approx 0.35<p_{6}$ which means this is a stable node.

This is illustrated in Figure 2 using the stated above parameter values. The starting values were selected so that $a_{0}>0$, $r_{0}>0$ and $o_{0}>0$ to imitate real initial hormone levels.

2.2 Lyapunov stability analysis

In this section, we show the stability of the fixed point by using the Lyapunov function approach. We consider the system (2.1)-(2.3) and denote

[TABLE]

where $(a^{*},r^{*},o^{*})$ is the fixed point ($a^{*}=o^{*}$).

Lemma 2.2** (Stability).**

Assume that

[TABLE]

and

[TABLE]

Then there exist $W^{*}>0$, $A_{0}>0$ and $p^{*}_{4}>0$ such that

[TABLE]

provided $W(0)<W^{*},\ 0\leqslant A<A_{0},\ p_{4}>p_{4}^{*},$ hence, the fixed point $(a^{*},r^{*},o^{*})$ is globally stable. If $p_{4}\leqslant p_{4}^{*}$ then there exist $A_{0}\leqslant A_{1}<A_{2}$ such that (2.11) holds provided $W(0)<W^{*},\ A_{1}<A<A_{2}.$

Proof of Lemma 2.2.

Using the system (2.1)–(2.3), we have

[TABLE]

As

[TABLE]

[TABLE]

[TABLE]

then

[TABLE]

i. e.

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

provided

[TABLE]

As $0\leqslant a^{*}=o^{*}\leqslant\frac{A}{p_{3}}$, $\frac{p_{5}}{p_{6}}\leqslant r^{*}\leqslant\frac{p_{5}+1}{p_{6}}$ and $o^{*}r^{*}\geqslant\frac{1}{2p_{2}}[\sqrt{1+\frac{4p_{2}p_{5}A}{p_{6}}}-1]$ then by (2.13) we get

[TABLE]

Hence,

[TABLE]

whence

[TABLE]

or

[TABLE]

As the function $F(A)$ has a unique minimum for positive $A$, denote by $A_{\min}$, then there exist $0<A_{1}<A_{\min}<A_{2}$ such that $F(A)<B$ provided $F(A_{\min})<B$.

So, if $W(0)<[\frac{\sqrt{\beta^{2}+4\alpha\gamma}-\beta}{2\gamma}]^{2}$ then by (2.12) we deduce that

[TABLE]

∎

3 Analysis of the model with time

delay

3.1 Stability analysis with respect to time delay

Note the fixed point $(a^{*},r^{*},o^{*})$ for (2.1)–(2.3) coincides with the one for (1.1)–(1.3). Let us denote by

[TABLE]

where $a_{\tau}=a(t-\tau)$, $r_{\tau}=r(t-\tau)$, and $o_{\tau}=o(t-\tau)$. Then $J_{\tau}$ at the point $(a^{*},r^{*},a^{*})$ is equal

[TABLE]

Now we will look for eigenvalues for the matrix $J^{*}+e^{-\lambda\tau}J^{*}_{\tau}$. So,

[TABLE]

whence we obtain the characteristic equation:

[TABLE]

Time delays are known to affect the stability of a fixed point. They can induce stability switches in which the zeros of the characteristic equation may cross the imaginary axis as the delay, $\tau$, increases. Looking at the characteristic equation as a function of $\tau$, and analyzing the location of the roots and the direction of motion as they cross the imaginary axis (see [5]). Destabilization will happen at critical values $\tau_{c}$ which is when there is a pair of purely imaginary characteristic values. Following the ideas of papers [5] and [3], let’s rewrite the characteristic equation as

[TABLE]

Then define

[TABLE]

We want to use Theorem 3 from the Appendix, so we need to check the following conditions:

$P(\lambda)=(\lambda+1)(\lambda+p_{3})(\lambda+p_{6}-K_{2})+K_{3}(\lambda+p_{6})$ and $Q(\lambda)=K_{3}(\lambda+p_{6})$ have no common imaginary zeros since each $p_{i}$ are real values. 2. 2.

It is quick to see that $\overline{P(i\lambda)}=P(i\lambda)$ and $\overline{Q(i\lambda)}=Q(i\lambda)$ for real $\lambda$. 3. 3.

$P(0)+Q(0)=p_{3}(p_{6}-K_{2})+2K_{3}p_{6}\neq 0$ so this is an important restriction in order to use the Theorem 3. 4. 4.

Referring back to (2.10), we see that there are at most 3 roots of (3.1) if $\tau=0$. 5. 5.

(3.2) has at most 6 real zeros for real $y$.

Therefore, by Theorem 3 from the Appendix, if $F(y)$ has no positive roots, the system is stable for all $\tau\geqslant 0$. If $F(y)$ has a simple positive root $y_{0}$, then there exists a pair of purely imaginary roots $\pm iv_{0}$ such that $v_{0}=\sqrt{y_{0}}$. For this $v_{0}$, there is a countable sequence of $\{\tau_{0}^{n}\}$ of delays for which stability switches can occur. Also, there exists a positive $\tau_{c}$ such that the system is unstable for all $\tau>\tau_{c}$. Investigating this further, let $x=y^{2}$

[TABLE]

where $b_{1}=p_{6}^{2}-2K_{2}p_{6}+p_{3}^{2}-2K_{3}+K_{2}^{2}+1$, $b_{2}=K_{2}^{2}+-2K_{2}K_{3}+2K_{3}p_{3}-2K_{2}K_{3}p_{3}+p_{3}^{2}+K_{2}^{2}p_{3}^{2}-2K_{2}p_{6}+2K_{2}K_{3}p_{6}-2K_{2}p_{3}^{2}p_{6}+p_{6}^{2}-2K_{3}p_{6}^{2}+p_{3}^{2}p_{6}^{2}$, and $b_{3}=K_{2}^{2}p_{3}^{2}-2K_{2}K_{3}p_{3}p_{6}-2K_{2}p_{3}^{2}p_{6}+2K_{3}p_{3}p_{6}^{2}+p_{3}^{2}p_{6}^{2}$. Note that

[TABLE]

and

[TABLE]

Now analyzing the roots of (3.3),

• •

If $\Delta_{0}\leqslant 0$, then $F^{\prime}(0)\geqslant 0$ and $F(x)$ is monotonically non-decreasing. Further,

• –

if $F(0)>0$, then $F$ has no positive roots and all the roots of the characteristic will remain to the left of the imaginary axis for all $\tau>0$.

• –

if $F(0)<0$, then since $\lim_{x\rightarrow\infty}F(x)=\infty$, there is at least one positive root of $F$ and thus the roots of the characteristic equation can cross the imaginary axis.

• •

If $\Delta_{0}>0$ then $F$ has critical points

[TABLE]

and if $x_{c_{1}}>0$ and $F(x_{c_{1}})<0$, then $F$ has positive roots (see [5]).

Stability switches are possible for each positive simple root $x_{j}$ of (3.3) and the cross is from left to right if $F^{\prime}(v_{0})>0$, and from right to left if $F^{\prime}(v_{0})<0$ according to Theorem 1 (see [5]). Now let’s analyze the characteristic quasi-polynomial (3.1) for $\lambda=iv$:

[TABLE]

where

[TABLE]

So $x_{j}$ ($j=1,2,3$) is a positive root of $F(x)=0$ and $v_{j}=\sqrt{x_{j}}$. Then $v_{j}$ satisfies (3.1) if its a solution to the system

[TABLE]

This yields

[TABLE]

provided $\max\{|A_{1}(v)A_{3}(v)-A_{2}A_{4}(v)|,|A_{1}(v)A_{2}-A_{3}(v)A_{4}(v)|\}\leqslant A_{2}^{2}+A_{3}^{2}(v)$,

Therefore, for every positive root $v_{j}$, it yields the following sequence of delays $\{\tau_{j}^{n}\}$ for which there are pure imaginary roots (3.1):

[TABLE]

As a result the following statement holds.

Lemma 3.1**.**

The system (2.1)–(2.3) with delay and $p_{3}(p_{6}-K_{2})+2K_{3}p_{6}\neq 0$ is stable for all $\tau\geqslant 0$ if $F(0)>0$ and $\Delta_{0}\leqslant 0$ (where $\Delta_{0}$ is from 3.4). The system has stability switches at some $\{\tau_{j}^{n}\}$ for every positive root $v_{j}$ of (3.1). Furthermore, if $A=0$, $p_{2}=0$ or $p_{5}=p_{6}=0$ then the delay has no affect on the stability of the system.

Example 3.1**.**

This example illustrates the dynamics of eigenvalues with respect to the time delay for the following set of parameters $p_{3}=0.41$, $p_{6}=0.91$, $K_{2}=0.81$, and $K_{3}=0.41$ in the equation (3.1). Taking the real and imaginary parts, we rewrite the equation as a system

[TABLE]

The red lines in Figure 3 represent the solution curves for the first equation and the blue lines in Figure 3 represent the solution curves for the second equation for different values of delay $\tau$. The eigenvalues, which are roots of (3.1), correspond to intersections between the red and blue lines.

When there is no delay, i.e. $\tau=0$, we only have three eigenvalues $\lambda\approx-0.9,-0.2\pm 0.8i$ (see Figure 3). When delay is is non-zero, countably many eigenvalues originate from $-\infty$ and move toward the imaginary axis as $\tau$ increases (see Figure 3). The eigenvalues can cross the imaginary axis only at the points $y_{1}\approx\pm 0.7$ and $y_{2}\approx\pm 0.25$ which are real roots of the equation (3.2) (see Figure 4). The density of complex eigenvalues around these crossing points $y_{1},y_{2}$ is increasing as the $\tau$ gets larger (see Figure 3).

When the delay $\tau<\tau^{*}\approx 2$ (where $\tau^{*}$ is a critical value found as a solution of (3.6) with $v_{1}=\sqrt{|y_{1}|}$) all eigenvalues are stable. The first stability switch happens at $\tau^{*}\approx 2$ when two complex conjugate eigenvalues cross the imaginary axis at $y_{1}\approx\pm 0.7$ changing the sign of the real part from negative to positive. At a later time $\tau^{*}\approx 11$ (where this $\tau^{*}$ is a critical value found as a solution of (3.6) with $v_{2}=\sqrt{|y_{2}|}$) this complex pair will cross the imaginary axis back changing the sign of the real part from positive to negative (see Figure 4).

Solving (3.6) and taking into account the periodicity of the arctangent function one can obtained the infinite sequences of delays associated with $v_{1}$ and another infinite sequence associated with $v_{2}$ at which stability switches may happen. At time delays associated with $v_{1}$ a complex conjugate pair of eigenvalues may cross the imaginary axis from left to right and for time delays associated with $v_{2}$ the pair may cross the imaginary axis from right to left. If the derivative of $F(y)$ (see (3.1)) does not change sign at the corresponding $\tau^{*}$ from either of the two sequences above, then the crossing of the imaginary axis does not happen.

3.2 Global in time existence of solutions

In this section, by Picard’s method we prove the existence of solutions to the problem (1.1)–(1.5).

Theorem 1**.**

If $A>0$, $p_{i}>0$, and

[TABLE]

then the problem (1.1)–(1.5) has a unique non-negative solution $(a(t)$, $r(t)$, $o(t))$ in $C^{2}$ for all $t>0$. Moreover, there exists a time $T^{*}>0$ such that

[TABLE]

For example, if $a_{\tau}(t)=a_{0}+\Lambda\,t^{2}e^{-t}$ then we get

[TABLE]

Proof of Theorem 1.

We will construct a solution to (1.1)–(1.5) by the iterative process. First of all, we will look for a solution on the interval $[0,\tau]$. From (1.3) we obtain that

[TABLE]

So, by (1.5) and (3.8) we have $o(t)\in C^{2}[0,\tau]$, and

[TABLE]

for all $t\in[0,\tau]$. Integrating (1.1) and (1.2) on the interval $[0,\tau]$, taking into account (3.8), we arrive at

[TABLE]

[TABLE]

for all $t\in[0,\tau]$. By (3.10), (3.11) we find that

[TABLE]

[TABLE]

for all $t\in[0,\tau]$. As a result, estimates (3.9), (3.12), (3.13) imply positivity of $o(t),\,a(t),\,r(t)$ on $[0,\tau]$ provided $o_{0}>0$, $r_{0}>0$, and $a_{\tau}(0)>0$. As the right-hand side of (3.11) is Lipschitz continuous on $r$ then there exists a unique solution of (3.11) on whole interval $[0,\tau]$ and, as a result, the one of (3.10). Moreover, obviously the solution $o(t),\,a(t),\,r(t)\in C^{2}[0,\tau]$ if the following fitting condition is true:

[TABLE]

Let us denote the corresponding solution to (3.8), (3.12), (3.13) on $[0,\tau]$ by $(o_{1}(t),a_{1}(t),r_{1}(t))$. Now we will find a solution on the interval $[\tau,2\tau]$. By (1.1)–(1.3) we get

[TABLE]

[TABLE]

[TABLE]

for all $t\in[\tau,2\tau]$. The system (3.15)–(3.17) has a unique solution $o_{2}(t),\,a_{2}(t),\,r_{2}(t)\in C^{2}[\tau,2\tau]$. Moreover,

[TABLE]

[TABLE]

[TABLE]

for all $t\in[\tau,2\tau]$. Continuing this iteration procedure, we derive

[TABLE]

[TABLE]

[TABLE]

for all $t\in[(k-1)\tau,k\,\tau]$, $k\in\mathbb{N}$, where $a_{0}(t)=a_{\tau}(t)$. This system has a unique solution $o_{k}(t)$, $a_{k}(t)$, $r_{k}(t)$ $\in C^{k+1}[(k-1)\tau,k\,\tau]$. Moreover,

[TABLE]

[TABLE]

[TABLE]

for all $t\in[(k-1)\tau,k\,\tau]$. As a result, the problem (1.1)–(1.4) has unique global in time solution, and, letting $k\to+\infty$, we get

[TABLE]

[TABLE]

∎

3.3 Existence of periodic solutions

Theorem 2**.**

Under the conditions of Theorem 1, the system (1.1)–(1.5) has at least one $C^{2}$-smooth $T$-periodic solution, where $T\neq\tau$.

Proof of Theorem 2.

The main line of proof follows (see [12, pp. 278–280]) (see also [13, Theorem 5]). Rewrite the system of (1.1)–(1.3) in the following form

[TABLE]

where

[TABLE]

[TABLE]

Obviously, the right-hand side of (3.27) is $T$-periodic with respect to $t$ as it does not depend on time explicitly. Without loss of generality, we may assume that

[TABLE]

This is true because otherwise we could represent the $\tau$ in the form

[TABLE]

Then shift to the auxilliary equation

[TABLE]

the $T$-periodic solutions of which coincide with the $T$-periodic solutions of (3.27).

On the set of all vector-valued functions $\textbf{x}(t)$ defined on $[0,T]$, let us define an operator $S_{\tau}$ by

[TABLE]

Note that the $T$-periodic solutions of (3.27) coincides with the solutions of the following integral equations:

[TABLE]

The operator $T(\tau,\textbf{x})$ maps every continuous vector-valued function $\textbf{x}(t)$ into a continuous vector-valued function for $0\leqslant t\leqslant T$, therefore $T(\tau,\textbf{x})$ is compact in $C$. Next, we will show that for all $T$-periodic solutions $\textbf{x}_{p}(t)$ there exists $R>0$ such that

[TABLE]

Really, from (3.28) we deduce that

[TABLE]

From here, using Grönwall’s lemma, we arrive at

[TABLE]

where $a=\frac{(A^{2}+(p_{5}+1)^{2})^{\frac{1}{2}}}{(p_{3}^{2}+p_{6}^{2}+1)^{\frac{1}{2}}+1}$, $b=(p_{3}^{2}+p_{6}^{2}+1)^{\frac{1}{2}}+1$. Hence, (3.29) holds with $R=(|\textbf{x}_{p}(0)|+a)e^{bT}-a$. As a result, by the fixed point theorem the integral equation (3.28) has at least one solution, and consequently the equation (3.27) has at least one $T$-periodic solution. ∎

3.4 Periodic solutions with the period $T=\tau$

Lemma 3.2**.**

If

[TABLE]

then the problem (1.1)–(1.4) has at least one $\tau$-periodic solution.

Example 3.2**.**

Let $A=1$, $p_{2}=11$, $p_{3}=1.2$, $p_{4}=0.05$, $p_{5}=0.11$, and $p_{6}=2.9$. Then the initial conditions are $a_{0}=\frac{A}{p_{3}}\left(1+p_{2}\sqrt{\frac{p_{4}(p_{6}r_{0}-p_{5})}{p_{5}+1-p_{6}r_{0}}}\right)^{-1}$, $r_{0}=\frac{1}{2}\left(\frac{p_{5}}{p_{6}}+\frac{p_{5}+1}{p_{6}}\right)$, and $o_{0}=a_{0}$. With these parameter values we solve (1.1)–(1.4) numerically using the Matlab solver dde23 [15]. The resulting periodic solutions can be seen in Figure 5.

If we perturb the parameters by a bit, the periodicity changes. We illustrate a periodicity change by using the parameters $A=1$, $p_{2}=7$, $p_{3}=1.2$, $p_{4}=0.05$, $p_{5}=0.51$, and $p_{6}=3.1$. Then the initial conditions are $a_{0}=\frac{A}{p_{3}}\left(1+p_{2}\sqrt{\frac{p_{4}(p_{6}r_{0}-p_{5})}{p_{5}+1-p_{6}r_{0}}}\right)^{-1}$, $r_{0}=\frac{1}{2}\left(\frac{p_{5}}{p_{6}}+\frac{p_{5}+1}{p_{6}}\right)$, and $o_{0}=a_{0}$. The resulting periodic solutions can be seen in Figure 6.

Periodicity of solutions can also be illustrated by plotting delayed function versus no delay function or function versus derivative as seen in Figure 7.

Proof of Lemma 3.2.

Note that if a function $\Phi(t)\in C^{2}[a,b]$ is periodic with a period $T>0$ then $\Phi(t+T)=\Phi(t)$ and $\Phi^{\prime}(t+T)=\Phi^{\prime}(t)$. Let $(a(t),r(t),o(t))$ be a $T$-periodic solution of (1.1)–(1.3) from Theorem 2. Then for this solution we have

[TABLE]

[TABLE]

[TABLE]

whence

[TABLE]

Let $T=\tau>0$. Then from (3.33) for $t\in[0,\tau]$ we have

[TABLE]

By (3.34) at $t=0$ we get

[TABLE]

[TABLE]

[TABLE]

whence

[TABLE]

[TABLE]

Hence

[TABLE]

[TABLE]

[TABLE]

As $f_{i}(0)=0$ and

[TABLE]

[TABLE]

[TABLE]

then $f_{i}(\tau)=0$ for all $\tau\geqslant 0$ provided

[TABLE]

∎

4 Discussion

Existence of non-negative solutions, uniqueness of a steady state and its stability were analyzed for the minimal model of the HPA in [21]. The analytical difference between the model without delay that we studied and the minimal model is the type of the nonlinearities in the equation that describes production and degradation of the adrenocorticotropic hormone $a(t)$ and in the equation for the density $r(t)$ of glucocorticoid receptor. The nonlinear functions in the model that we analyze depend on the product $o(t)r(t)$ but in the minimal model of the HPA they depend on $o(t)$ only. We obtained similar results but for more complicated non-linearity terms. We also analyzed stability for all possible cases of parameter ranges. We added Lyapunov stability analysis to show the non-linear stability result and we believe that non-linear stability analysis has never done before for this type of model.

For the model with delay we obtained a critical time delay value when the originally stable system becomes unstable as a pair of complex eigenvalues crosses the imaginary axis. To illustrate the dynamics of complex eigenvalues with respect to time delay we used intersections of zero level sets between real and imaginary parts of the characteristic equation as a first approximation of a complex eigenvalue and after that we used Newton iterations to improve the accuracy. We believe that stability analysis for HPA model with respect to time delay was not done before and our results are new in this area.

For certain parameter values we rigorously proved existence of non-negative periodic solutions for a given period $T$ (under an assumption that a given period does not coincide with the value of time delay) and illustrated different periodic solutions numerically. Our numerical simulations revealed that the period of the solution is very sensitive to small perturbations of parameter values.

5 Appendix

Lemma 5.1** (Routh Hurwitz Criteria for a Nonlinear System).**

Suppose

[TABLE]

Suppose $\textbf{x}_{s}$ is a fixed point of (5.1) and the characteristic polynomial at the fixed point is

[TABLE]

If $\alpha_{1}>0,\,\alpha_{3}>0$ and $\alpha_{1}\alpha_{2}>\alpha_{3}$, then the fixed point is asymptotically stable. If $\alpha_{1}<0,\,\alpha_{3}<0$ or $\alpha_{1}\alpha_{2}<\alpha_{3}$, then the fixed point is unstable.

Theorem 3**.**

(see [5]) Consider equation (3.1), where $P(z)$ and $Q(z)$ are analytic functions in a right half-plane $Re(z)>-\delta$, $\delta>0$, which satisfy the following conditions:

1. i)

$P(z)$* and $Q(z)$ have no common imaginary zeros;* 2. ii)

$\overline{P(-iy)}=P(iy)$* and $\overline{Q(-iy)}=Q(iy)$ for real $y$;* 3. iii)

$P(0)+Q(0)\neq 0$; 4. iv)

There are at most a finite number of roots of (3.1) in the right half-plane when $\tau=0$; 5. v)

$F(y)\equiv|P(iy)|^{2}-|Q(iy)|^{2}$* for real $y$ has at most a finite number of real zeros.*

Under these conditions, the following statements are true.

1. a)

Suppose that the equation $F(y)=0$ has no positive roots. Then if (3.1) is stable at $\tau=0$ it remains stable for all $\tau\geq 0$, whereas if it is unstable at $\tau=0$ it remains unstable for all $\tau\geq 0$. 2. b)

Suppose that the equation $F(y)=0$ has at least one positive root and that each positive root is simple. As $\tau$ increases, stability switches may occur. There exists a positive number $\tau_{c}$ such that the equation (3.1) is unstable for all $\tau>\tau_{c}$. As $\tau$ varies from [math] to $\tau_{c}$, at most a finite number of stability switches may occur.