# Stability and Hopf bifurcation analysis for the hypothalamic-pituitary-adrenal axis model with memory

Stability and Hopf bifurcation analysis for the hypothalamic-pituitary-adrenal axis model with... Abstract This article generalizes the existing minimal model of the hypothalamic-pituitary-adrenal (HPA) axis in a realistic way, by including memory terms: distributed time delays, on one hand and fractional-order derivatives, on the other hand. The existence of a unique equilibrium point of the mathematical models is proved and a local stability analysis is undertaken for the system with general distributed delays. A thorough bifurcation analysis for the distributed delay model with several types of delay kernels is provided. Numerical simulations are carried out for the distributed delays models and for the fractional-order model with discrete delays, which substantiate the theoretical findings. It is shown that these models are able to capture the vital mechanisms of the HPA system. 1. Introduction The hypothalamuspituitaryadrenal (HPA) axis is a self-regulated dynamic feedback neuroendocrine system i.e. engaged in the rapid response to stressful stimuli and is responsible for the return to homeostasis through complex feedback mechanisms (Conrad et al., 2009; Swanson, 2000). By regulating the plasma levels of corticosteroids secreted from adrenal glands, it also controls many bodily processes, including mood and emotions, digestion, sexuality, the immune system, energy storage and expenditure. The HPA axis is organized into three distinct regions: the hypothalamus, pituitary gland and adrenal gland, with a complex set of direct influences and feedback interactions among the three endocrine glands. These glands work together by producing and secreting, or responding to common hormones including corticotropin-releasing hormone (CRH), corticotropin (ACTH), and cortisol (CORT) Kyrylov et al. (2005). Both physical stressors (e.g. infection, thermal exposure, dehydration) and psychological stressors (e.g. fear, anticipation) activate the hypothalamus to release CRH, which induces the ACTH production in the pituitary. Then, ACTH is transported by the blood to the adrenal cortex of the adrenal gland, where it stimulates the production of CORT, which in turn suppresses the production of both CRH and ACTH (see Fig. 1). Fig. 1. View largeDownload slide A simple schematic representation of the HPA axis with the negative feedback. Fig. 1. View largeDownload slide A simple schematic representation of the HPA axis with the negative feedback. It is important to keep CORT concentration within a certain physiological range. CORT overproduction, which is often due to a pituitary tumour causing high levels of circulating ACTH, leads to Cushings disease. CORT underproduction, which generates Addisons disease, is most frequently the result of adrenal destruction. In the past few decades, mathematical modelling has started to play an increasingly important role in the study of metabolic and endocrine processes, both in physiology and in clinical medicine (Bairagi et al., 2008). The theory of nonlinear dynamical systems has become a promising research tool for studying rhythmicity in humans. Moreover, it is often mandatory to introduce time delays in the mathematical models describing real life phenomena. Analysing a mathematical model of the HPA axis is helpful in understanding the simultaneous feedback mechanisms and revealing the different ways in which a malfunction can occur (Andersen et al., 2013). It is important to emphasize that time delays of up to 60 min (according to Boscaro et al. (1998); Posener et al. (1997)) unavoidably exist in the HPA axis, due to the transportation of the hormones among the three glands. Several mathematical models of the HPA axis have been developed and analysed by Jelić et al. (2005); Kyrylov et al. (2005); Lenbury & Pornsawad (2005); Savić & Jelić (2005); Savić et al. (2006); Bairagi et al. (2008); Conrad et al. (2009); Markovic et al. (2011); Vinther et al. (2011); Andersen et al. (2013); Pornsawad (2013); Gudmand-Hoeyer et al. (2014). Experimental data show the circadian as well as ultradian rhythm of hormone levels (Carroll et al., 2007; Veldhuis et al., 2008), which should be reflected by the mathematical models of the HPA axis, through the existence of oscillatory solutions. The ultradian rhythm is usually considered an inherent behaviour of the HPA axis, whereas the circadian rhythm is regarded as an external input to the axis (Andersen et al., 2013). A frequently studied model of the HPA axis, called ‘minimal model’ (Vinther et al., 2011), consists of a system of three coupled, non-linear differential equations, with the hormones CRH, ACTH and free CORT as variables. While the mathematical model investigated by Kyrylov et al. (2005) included two more differential equations, accounting for corticosteroid-binding globulin (CBG) bound CORT, and albumin bound CORT, (Vinther et al., 2011) pointed out that only free CORT is capable of interacting with the rest of the HPA-axis, which is the reason for considering only three differential equations in the formulation of the minimal model. If time delays are not incorporated in the mathematical model, no oscillatory behaviour has been detected (Savić & Jelić, 2005; Vinther et al., 2011; Andersen et al., 2013). For systems of ordinary differential equations, sufficiently large time delays are often linked with generating oscillatory solutions. A particular case of the minimal model with discrete time-delays and exponential negative-feedback mechanism has been studied in Lenbury & Pornsawad (2005) and later, in Pornsawad (2013), revealing the occurrence of Hopf bifurcations resulting in the emergence of periodic orbits. Other variants of the minimal model with Hill-type feedback functions and discrete time-delays have been considered in Savić & Jelić (2005); Savić et al. (2006), but no oscillations have been reported. More recently, Vinther et al. (2011) have included discrete time-delays in the minimal model and have observed periodic solutions by computer simulations. Their investigations show that time delays of at least 18 min in the feedback mechanisms are needed for generating oscillations. A slightly modified version of the minimal model with discrete time-delays has been analysed in Bairagi et al. (2008), successfully obtaining the desired Hopf bifurcation and thereby, oscillating solutions for sufficiently large time-delays. In this article, we generalize the existing minimal model of the HPA axis in two realistic ways: firstly, by including distributed time delays and secondly, by considering fractional-order derivatives. On one hand, distributed time delays represent the situation where the delays occur in certain ranges of values with some associated probability distributions, taking into account the whole past history of the variables. In many real world applications, distributed time delays are more realistic and more accurate than discrete time delays (Cushing, 2013). Distributed delay models appear in a wide range of applications such as, population biology (Ruan & Wolkowicz, 1996; Faria & Oliveira, 2008), haematopoiesis (Adimy & Crauste, 2003; Adimy et al., 2005, 2006; Özbay et al., 2008), neural networks (Jessop & Campbell, 2010). On the other hand, the main benefit of fractional-order models in comparison with classical integer-order models is that fractional derivatives provide a good tool for the description of memory and hereditary properties of various processes (Podlubny, 1999; Kilbas et al., 2006; Lakshmikantham et al., 2009). Due to the fact that the whole past history of the variables is accounted for in the formulation of both distributed time-delays as well as fractional-order derivatives, these generalizations of the minimal model are able to reflect biological variability in a better way than other approaches. The article is structured as follows. Section 2 provides the mathematical model of the HPA axis, where instead of considering the transportation of different hormones as an instantaneous process, we introduce distributed time delays to account for the time needed by the hormones to travel from source to destination. In Section 3, the existence of a unique equilibrium point of the system is shown. Local stability analysis of the system with general distributed delays is analysed in Section 4. In Section 5, we undertake a bifurcation analysis for the distributed delay model in the case of several types of delay kernels. The fractional order mathematical model of the HPA axis is presented and shortly analysed in Section 6. Numerical simulations are carried out and discussed in Section 7, followed by concluding remarks in Section 8. 2. Mathematical model of HPA with distributed delays In formulating the mathematical model which describes the variation in time of the concentrations of the three hormones CRH, ACTH and CORT, the following sequence of typical events is considered, according to the schema presented in Fig. 1. CRH is secreted from the hypothalamus and released into the portal blood vessel of the hypophyseal stalk, and then transported to the anterior pituitary where it stimulates the secretion of ACTH, with an average time delay $$\tau_1$$. Then, in the cortex of the adrenal glands, ACTH stimulates the secretion of the stress hormone CORT with the average time delay $$\tau_2$$. CORT has a negative feedback effect on the hypothalamus and the pituitary, expressed by two feedback functions $$f_1$$ and $$f_2$$, affecting the synthesis and release of CRH and ACTH, respectively. On one hand, CORT inhibits the secretion of CRH through glucocorticoid receptors (GRs) situated in the hypothalamus Landsberg et al. (1992), with an average time delay $$\tau_{31}$$. On the other hand, CORT also performs a negative feedback on the secretion of ACTH through GRs situated in the pituitary, with an average time delay $$\tau_{32}$$. The hormone concentrations of CRH, ACTH and CORT are depleted through the rate constants $$w_1$$, $$w_2$$ and $$w_3$$, respectively. The mathematical model of the HPA axis studied in this article is based on the minimal model introduced in Vinther et al. (2011). The main improvement is that we consider distributed time delays to account for the transport of the hormones among the glands. Denoting the hormone concentrations, for simplicity, by $$CRH(t)=x_1(t)$$, $$ACTH(t)=x_2(t)$$, $$CORT(t)=x_3(t)$$, the following system of differential equations with distributed delays is considered:   {x˙1(t)=f1(∫−∞tx3(s)h31(t−s)ds)−w1x1(t),x˙2(t)=f2(∫−∞tx3(s)h32(t−s)ds)∫−∞tx1(s)h1(t−s)ds−w2x2(t),x˙3(t)=k3∫−∞tx2(s)h2(t−s)ds−w3x3(t), (2.1) where all the first terms on the right hand side represent production and all the second terms represent depletion of hormones. The constant $$k_3$$ as well as the elimination constants $$w_1,w_2,w_3$$ are positive. The functions $$f_1,f_2:[0,\infty)\to (0,\infty)$$, which represent the negative feedback from CORT on CRH and ACTH, respectively, are assumed to be strictly decreasing, smooth and bounded on $$[0,\infty)$$. In particular, the results presented in this article are also applicable when Hill functions are being used in the expression of the feedback functions (Vinther et al., 2011; Andersen et al., 2013):   f1(u)=k1(1−ηuαcα+uα),f2(u)=k2(1−μuαcα+uα), (2.2) with $$\alpha\geq 1$$, $$k_1,k_2>0$$, $$\eta,\mu\in(0,1)$$, $$c>0$$. It is easy to verify that functions (2.2) satisfy all the properties mentioned above. However, it may be possible to model the negative feedback using different types of functions $$f_1$$ and $$f_2$$. In this article, our aim is to obtain general results which will also be applicable to other choices of negative feedback functions, besides functions (2.2), often used in the literature. In system (2.1), the delay kernels $$h_1,h_2,h_{31},h_{32}:[0,\infty)\to[0,\infty)$$ are probability density functions representing the probability that a particular time delay occurs. They are assumed to be bounded, piecewise continuous and satisfy   ∫0∞h(s)ds=1. (2.3) The average delay of a delay kernel $$h(t)$$ is given by   τ=∫0∞sh(s)ds<∞. Two important classes of delay kernels, which are often used in the literature, are worth mentioning: Dirac kernels: $$h(s)=\delta(s-\tau)$$, where $$\tau\geq 0$$. In this particular case, the distributed delay is reduced to a discrete time delay:   ∫−∞tx(s)h(t−s)ds=∫0∞x(t−s)δ(s−τ)ds=x(t−τ). Gamma kernels: $$h(s)=\displaystyle\frac{ s^{p-1}e^{-s/\beta}}{\beta^p\Gamma(p)}$$, where $$p,\beta>0$$. The average delay of a Gamma kernel is $$\tau=p\beta$$. The analysis of the mathematical models including particular classes of delay kernels (such as weak Gamma kernels with $$p=1$$ or strong Gamma kernels with $$p=2$$) may shed a light on how distributed delays affect the dynamics differently from discrete delays. However, in the modelling of real world phenomena, one usually does not have access to the exact distribution, and approaches using general kernels may be more useful (Bernard et al., 2001; Campbell & Jessop, 2009; Yuan & Bélair, 2011; Diekmann & Gyllenberg, 2012). Initial conditions associated with system (2.1) are as follows:   xi(s)=φi(s),∀s∈(−∞,0],i=1,2,3, where $$\varphi_i$$ are bounded continuous functions defined on $$(-\infty,0]$$, with values in $$[0,\infty)$$. 3. Existence of a unique equilibrium point An equilibrium point of system (2.1) is a solution of the following algebraic system:   {f1(x3)=w1x1,f2(x3)x1=w2x2,k3x2=w3x3, (3.1) which is equivalent to   {x1=f1(x3)w1,x2=w3k3x3,x3=k3w1w2w3f1(x3)f2(x3). (3.2) Due to the properties of $$f_1$$ and $$f_2$$, the function $$\displaystyle\frac{k_3}{w_1w_2w_3}f_1(x)f_2(x)$$ which appears in the right hand side of the last equation of system (3.2), is strictly positive and strictly decreasing on $$[0,\infty)$$, and therefore, it has a unique fixed point $$x^\star>0$$. It follows that system (2.1) has a unique equilibrium point   E=(f1(x⋆)w1,w3x⋆k3,x⋆). (3.3) In the following, necessary and sufficient conditions will be explored for the local asymptotic stability of the equilibrium point $$E$$ and the occurrence of limit cycles in a neighbourhood of $$E$$ (due to Hopf bifurcations) that can explain the ultradian rhythm. 4. Local stability analysis of system (2.1) In this section, considering general delay kernels, we seek to obtain delay independent sufficient conditions for the local asymptotic stability of the equilibrium point $$E$$. Such results prove to be useful if one is unable to accurately estimate the time delays in system (2.1). Using the transformation $$y_1(t)=x_1(t)-\displaystyle\frac{f_1(x^\star)}{w_1}$$, $$y_2(t)=x_2(t)-\displaystyle\frac{w_3x^\star}{k_3}$$ and $$y_3(t)=x_3(t)-x^\star$$, the linearized system of (2.1) at the equilibrium point $$E$$ is:   {y˙1(t)=f1′(x⋆)∫−∞ty3(s)h31(t−s)ds−w1y1(t),y˙2(t)=f2(x⋆)∫−∞ty1(s)h1(t−s)ds+1w1f1(x⋆)f2′(x⋆)∫−∞ty3(s)h32(t−s)ds−w2y2(t),y˙3(t)=k3∫−∞ty2(s)h2(t−s)ds−w3y3(t). (4.1) The associated characteristic equation of the linearized system (4.1) is:   (z+w1)(z+w2)(z+w3)+a(z+w1)H2(z)H32(z)+bH1(z)H2(z)H31(z)=0, (4.2) where $$H_i( z)=\int_{0}^\infty e^{- z s}h_i(s)ds$$ represent the Laplace transforms of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$ and   a=−k3w1f1(x⋆)f2′(x⋆)=−w2w3x⋆f2′(x⋆)f2(x⋆)>0, (4.3)  b=−k3f1′(x⋆)f2(x⋆)=−w1w2w3x⋆f1′(x⋆)f1(x⋆)>0. (4.4) The following inequalities will be useful for the theoretical analysis:   (I1)8f1(x⋆)+x⋆f1′(x⋆)≥0;(I2)1+x⋆f1′(x⋆)f1(x⋆)+x⋆f2′(x⋆)f2(x⋆)>0;(I2¯)1+x⋆f1′(x⋆)f1(x⋆)+x⋆f2′(x⋆)f2(x⋆)≤0. Remark 4.1 Consider the feedback function   f1(u)=k1(1−ηuαcα+uα) with $$k_1>0$$, $$\eta\in(0,1)$$, $$c>0$$, $$\alpha\geq 1$$. A straightforward computation yields:   8f1(u)+uf1′(u)=k18(1−η)u2α+(16−8η−αη)uαcα+8c2α(cα+uα)2. It can be easily seen that if $$\alpha\leq 8$$, we have   8f1(u)+uf1′(u)≥k18(1−η)u2α+16(1−η)uαcα+8c2α(cα+uα)2≥0, and hence, the inequality $$8f_1(u)+uf_1'(u)\geq 0$$ is satisfied for any $$u\geq 0$$ (regardless of the choice of parameters $$k_1$$, $$\eta$$ or $$c$$). Therefore, if $$\alpha\leq 8$$, inequality $$(I_1)$$ holds as well. It is worth mentioning that according to Murray (2002), $$\alpha>7$$ is considered unphysiological. Remark 4.2 Consider the particular case of feedback functions given by (2.2), with $$k_1,k_2>0$$, $$\eta,\mu\in(0,1)$$, $$\alpha\geq 1$$ and with the constant $$c$$ given by   c=x⋆. This assumption comes from the fact that the constant $$c$$ is chosen to be equal to the mean value of free CORT (Vinther et al., 2011), which, in turn, is equal to the last component of the equilibrium point $$E$$. In this case, the term that appears in the left hand-side of inequalities $$(I_2)$$ and $$(\overline{I_2})$$ becomes   1+x⋆f1′(x⋆)f1(x⋆)+x⋆f2′(x⋆)f2(x⋆)=1−α2[η2−η+μ2−μ]. It is important to note that in this case, inequalities $$(I_2)$$ and $$(\overline{I_2})$$ only depend on the parameters $$\alpha,\eta,\mu$$, and do not depend on the choice of the parameters $$k_1,k_2$$. For instance, if $$\eta=\mu=0.5$$, inequality $$(I_2)$$ is equivalent to $$\alpha<3$$. Theorem 4.1 (Local asymptotic stability) (1) In the non-delayed case, if inequality $$(I_1)$$ is satisfied, then the equilibrium point $$E$$ of system (2.1) is locally asymptotically stable. (2) For any delay kernels $$h_i(t)$$, $$i\in\{1,2,31,32\}$$, if inequality $$(I_2)$$ holds, then the equilibrium point $$E$$ of system (2.1) is locally asymptotically stable. Proof 1. In the non-delayed case, the characteristic equation (4.2) becomes:   λ3+c1λ2+c2λ+c3=0, (4.5) where   c1=w1+w2+w3>0,c2=w1w2+w2w3+w1w3+a>0,c3=w1w2w3+aw1+b>0. We can easily compute   c1c2−c3=a(w2+w3)+γ(w1,w2,w3)+2w1w2w3−b, where   γ(w1,w2,w3)=w12w2+w1w22+w12w3+w1w32+w2w32+w22w3. By the inequality of arithmetic and geometric means, we deduce   γ(w1,w2,w3)≥6w1w2w3,for any w1,w2,w3>0, and hence, based on (4.4) and $$(I_1)$$, we obtain:   c1c2−c3≥a(w2+w3)+8w1w2w3−b==a(w2+w3)+w1w2w3(8+x⋆f1′(x⋆)f1(x⋆))>0. By the Routh–Hurwitz stability test, the equilibrium point $$E$$ is asymptotically stable. 2. From (4.3) and (4.4) it can be easily seen that inequality $$(I_2)$$ is equivalent to   a+bw1<w2w3. The characteristic equation (4.2) can be written as   φ(z)=ψ(z), where the functions $$\varphi$$ and $$\psi$$ are given by   φ(z)=−(z+w1)(z+w2)(z+w3),ψ(z)=a(z+w1)H2(z)H32(z)+bH1(z)H2(z)H31(z). These functions are holomorphic in the right half-plane. Let $$z\in \mathbb{C}$$ with $$\Re(z)\geq 0$$. For any $$i\in\{1,2,31,32\}$$, from (2.3) we obtain:   |Hi(z)|=|∫0∞e−zshi(s)ds|≤∫0∞|e−zs|hi(s)ds=∫0∞e−ℜ(z)shi(s)ds≤∫0∞hi(s)ds=1, and hence, we have:   |ψ(z)|≤a|z+w1||H2(z)||H32(z)|+b|H1(z)||H2(z)||H31(z)|≤a|z+w1|+b=|z+w1|(a+b|z+w1|)=|z+w1|(a+b|z|2+w12+2ℜ(z)w1)≤|z+w1|(a+bw1)<|z+w1|w2w3≤|z+w1|[(|z|2+w22+2ℜ(z)w2)(|z|2+w32+2ℜ(z)w3)]1/2=|z+w1||z+w2||z+w3|=|φ(z)|. Hence, the inequality $$|\psi(z)|<|\varphi(z)|$$ holds for any $$z\in \mathbb{C}$$, $$\Re(z)\geq 0$$. Therefore, the characteristic equation $$\varphi(z)=\psi(z)$$ does not have any root in the right half-plane (or the imaginary axis). This means that all the roots of the characteristic equation (4.2) have strictly negative real part, and the equilibrium $$E$$ is asymptotically stable. ■ Corollary 4.1 For any delay kernels $$h_i(t)$$, $$i\in\{1,2,31,32\}$$, if the equilibrium point $$E$$ of system (2.1) is unstable, then inequality $$(\overline{I_2})$$ holds. In other words, inequality $$(\overline{I_2})$$ is a necessary condition for the occurrence of bifurcations in system (2.1). Remark 4.3 If $$f_1$$ is given by (2.2) with $$\alpha\leq 8$$, it follows from Remark 4.1 that inequality $$(I_1)$$ holds. Based on Theorem 4.1, the equilibrium point $$E$$ is asymptotically stable in the non-delayed case. This improves the sufficient condition $$\alpha\leq 7.46$$ presented in Vinther et al. (2011). As $$\alpha>7$$ is considered unphysiological (Murray, 2002), the equilibrium point $$E$$ is locally asymptotically stable for all realistic values of the parameters if no time delays are considered. Moreover, if $$f_1,f_2$$ are given by (2.2) (such as in Remark 4.2), with $$c=x^\star$$, and the following inequality is satisfied   2α>η2−η+μ2−μ, the equilibrium point $$E$$ of system (2.1) is asymptotically stable for any choice of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$ and of the parameters $$k_1,k_2,k_3,w_1,w_2,w_3$$. In the special case $$\alpha=1$$ considered in (Savić & Jelić, 2005; Savić et al., 2006), it is easy to see that the above inequality is fulfilled for any $$\eta,\mu\in(0,1)$$, implying that the equilibrium point $$E$$ is locally asymptotically stable for any choice of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$ and of the parameters $$k_1,k_2,k_3,w_1,w_2,w_3$$. 5. Bifurcation analysis of system (2.1) The bifurcation analysis presented in this section takes into consideration the average time delays of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$. A first observation is that the time required by CRH to travel from the hypothalamus to the pituitary through the hypophyseal portal blood vessels is extremely short (Bairagi et al., 2008) and therefore, in most numerical simulations the average time delay $$\tau_1$$ is considered close to $$0$$. Moreover, $$\tau_{31}$$ and $$\tau_{32}$$ are comparable, as they represent the average time delays due to the negative feedback effect of the adrenal glands on the hypothalamus and pituitary, respectively, which are closely situated. For this reason, in this section, we will assume for simplicity that   H32(z)=H1(z)H31(z), and we denote   H(z)=H2(z)H32(z)=H1(z)H2(z)H31(z). In fact, $$H(z)$$ is the Laplace transform of the convolution of the delay kernels $$h_2$$ and $$h_{32}$$ defined as   h(t)=∫0th2(s)h32(t−s)ds, with the mean   τ=∫0∞sh(s)ds=τ2+τ32, (5.1) where $$\tau_2$$ and $$\tau_{32}$$ represent the average delays of the kernels $$h_2$$ and $$h_{32}$$ respectively. This results from the fact that the probability density function of the sum of two independent random variables is the convolution of their separate probability density functions. Therefore, the characteristic equation (4.2) becomes   (z+w1)(z+w2)(z+w3)+[a(z+w1)+b]H(z)=0, (5.2) or equivalently:   H(z)−1=Q(z), where   Q(z)=−aw1+b+az(z+w1)(z+w2)(z+w3). The function $$Q(z)$$ defined above will play an important role in the bifurcation analysis presented in this section. We summarize its most important properties in the following Lemma. Lemma 5.1 The function   ω↦|Q(iω)|=a2ω2+(aw1+b)2(ω2+w12)(ω2+w22)(ω2+w32) is strictly decreasing on $$[0,\infty)$$ and the equation   |Q(iω)|=1 has a unique positive real root $$\omega_0$$ if and only if inequality $$(\overline{I2})$$ is satisfied. Moreover, the following inequality holds:   ℑ(Q′(iω)Q(iω))>0∀ω>0. Proof We have   |Q(iω)|2=a2(ω2+w22)(ω2+w32)+2aw1+b2(ω2+w12)(ω2+w22)(ω2+w32) and it is easy to see that $$\omega\mapsto|Q(i\omega)|$$ is strictly decreasing on $$[0,\infty)$$, approaching $$0$$ as $$\omega\rightarrow\infty$$. Therefore, the equation $$|Q(i\omega)|=1$$ has a unique solution if and only if $$|Q(0)|>1$$, or $$w_1w_2w_3<aw_1+b$$, which is equivalent to inequality $$(\overline{I2})$$. Moreover, we have:   ddω|Q(iω)|2=ddω[Q(iω)Q(iω)¯]=2ℜ[Q(iω)¯ddωQ(iω)]=2ℜ[iQ(iω)¯Q′(iω)]=−2ℑ[Q(iω)¯Q′(iω)]=−2|Q(iω)|2ℑ(Q′(iω)Q(iω)). As $$\omega\mapsto |Q(i\omega)|^2$$ is strictly decreasing on $$(0,\infty)$$, its derivative is strictly negative, and hence, $$\Im\left(\displaystyle\frac{Q'(i\omega)}{Q(i\omega)}\right)>0$$, for any $$\omega>0$$. ■ Remark 5.1 A simple computation shows that   ℜ[Q(iω)]=aω4+[b(w1+w2+w3)+a(w12−w2w3)]ω2−w1w2w3(aw1+b)(ω2+w12)(ω2+w22)(ω2+w32). This formula will be useful in the framework of the bifurcation results that follow. Due to the high complexity of the problem, we are unable to perform the bifurcation analysis for general kernels $$h_i$$, $$i\in\{1,2,31,32\}$$. Thus, we focus our attention on the following cases: all delay kernels are Dirac kernels; all delay kernels are Gamma kernels; some delay kernels are Dirac kernels while others are Gamma kernels. 5.1 Dirac kernels Consider that all the delay kernels are Dirac kernels: $$h_1(t)=\delta(t-\tau_1)$$, $$h_2(t)=\delta(t-\tau_2)$$, $$h_{31}(t)=\delta(t-\tau_{31})$$, $$h_{32}(t)=\delta(t-\tau_{32})$$, where $$\tau_1,\tau_2,\tau_{31},\tau_{32}\geq 0$$ satisfy the property   τ2+τ32=τ1+τ2+τ31=τ>0. (5.3) In this case, the characteristic equation (5.2) becomes:   (z+w1)(z+w2)(z+w3)+[a(z+w1)+b]e−τz=0, (5.4) or equivalently:   eτz=Q(z). We choose $$\tau$$ as bifurcation parameter. Theorem 5.1 (Hopf bifurcations in the case of Dirac kernels) Assume that inequalities $$(I_1)$$ and $$(\overline{I_2})$$ are satisfied. For any $$p\in\mathbb{Z}^+$$, consider   τp=arccos⁡[ℜ(Q(iω0))]+2pπω0, (5.5) where $$\omega_0>0$$ is given by Lemma 5.1. The equilibrium point $$E$$ is asymptotically stable if any only if $$\tau\in[0,\tau_0)$$. For any $$p\in\mathbb{Z}^+$$, at $$\tau=\tau_p$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. Proof Equation (5.4) has a pair of complex conjugated solutions $$z=\pm i\omega$$ on the imaginary axis ($$\omega>0$$) if and only if   eiωτ=Q(iω). (5.6) Taking the absolute value in (5.6) we obtain that $$|Q(i\omega)|=1$$, and hence, we obtain that $$\omega=\omega_0$$, where $$\omega_0$$ is the unique positive real solution given by Lemma 5.1. From (5.6) we can also deduce   cos⁡(τω0)=ℜ[Q(iω0)], which leads us to the values $$\tau_p$$ given by (5.5). From Theorem 4.1, we know that the equilibrium point $$E$$ is asymptotically stable if $$\tau=0$$ (since inequality $$(I_1)$$ holds). The number of the roots of the characteristic equation from the left half-plane can change only if a root (or pair of complex conjugated roots) crosses the imaginary axis, or more precisely, whenever $$\tau=\tau_p$$, $$p\in\mathbb{Z}^+$$ (in which case, $$\pm i \omega_0$$ are roots of the characteristic equation). Therefore, for any $$\tau\in[0,\tau_0)$$ the equilibrium point $$E$$ is asymptotically stable. Let $$z_p(\tau)$$ denote the root of the characteristic equation (5.4) satisfying $$z_p(\tau_p)=i\omega_0$$. The function $$z_p(\tau)$$ satisfies   eτzp(τ)=Q(zp(τ)). Taking the derivative with respect to $$\tau$$, it follows that   (τzp′(τ)+zp(τ))eτzp(τ)=zp′(τ)Q′(zp(τ)). We obtain   zp′(τ)=zp(τ)Q′(zp(τ))e−τzp(τ)−τ=zp(τ)Q′(zp(τ))Q(zp(τ))−τ, and hence:   zp′(τp)=iω0Q′(iω0)Q(iω0)−τp. Taking the real part and using Lemma 5.1, we obtain:   dℜ(zp)dτ|τ=τp=ω0ℑ(Q′(iω0)Q(iω0))|Q′(iω0)Q(iω0)−τp|2>0. This nondegeneracy condition for the Hopf bifurcation shows that the equilibrium point $$E$$ can only be asymptotically stable if and only if $$\tau\in[0,\tau_0)$$ and that for any $$p\in\mathbb{Z}^+$$, at $$\tau=\tau_p$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. ■ 5.2 Gamma kernels We now consider that the delay kernels are Gamma distribution kernels: $$h_1(t)=\displaystyle\frac{t^{n_1-1}e^{-t/\beta}}{\beta^{n_1}(n_1-1)!}$$, $$h_2(t)=\displaystyle\frac{t^{n_2-1}e^{-t/\beta}}{\beta^{n_2}(n_2-1)!}$$, $$h_{31}(t)=\displaystyle\frac{t^{n_{31}-1}e^{- t/\beta}}{\beta^{n_{31}} (n_{31}-1)!}$$, $$h_{32}(t)=\displaystyle\frac{ t^{n_{32}-1}e^{-t/\beta}}{\beta^{n_{32}}(n_{32}-1)!}$$, where $$\beta>0$$ and $$n_1,n_2,n_{31},n_{32} \in\mathbb{Z}^+\setminus\{0\}$$ satisfy:   n2+n32=n1+n2+n31=n≥2. The characteristic equation (4.2) becomes:   (z+w1)(z+w2)(z+w3)+a(z+w1)+b(βz+1)n=0, (5.7) or equivalently   (βz+1)n=Q(z). Choosing $$\beta$$ as bifurcation parameter, we obtain the following result: Theorem 5.2 (Hopf bifurcations in the case of Gamma kernels) Assume that inequalities $$(I_1)$$ and $$(\overline{I_2})$$ are satisfied. Let $$\omega_n$$ denote the largest real root of the equation   Tn(1|Q(iω)|1/n)=ℜ(Q(iω))|Q(iω)| (5.8) from the interval $$(0,\omega_0)$$, where $$T_n$$ is the Chebyshev polynomial of the first kind of order $$n$$, and consider   βn=1ωn|Q(iωn)|2/n−1. (5.9) The equilibrium point $$E$$ is asymptotically stable if $$\beta\in(0,\beta_n)$$. At $$\beta=\beta_n$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. Proof Equation (5.7) has a pair of complex conjugated solutions $$z=\pm i\omega$$ on the imaginary axis ($$\omega>0$$) if and only if   (iβω+1)n=Q(iω). (5.10) Applying the modulus to both sides of equation, we obtain:   (β2ω2+1)n=|Q(iω)|2, which means that $$|Q(i\omega)|>1$$. Based on Lemma 5.1, we deduce that $$\omega\in(0,\omega_0)$$ and:   β=1ω|Q(iω)|2/n−1. (5.11) Replacing (5.11) in (5.10) we can also deduce   (1|Q(iω)|1/n+i1−1|Q(iω)|2/n)n=Q(iω)|Q(iω)|. Denoting $$\theta(\omega)=\arccos\left(\displaystyle\frac{1}{|Q(i\omega)|^{1/n}}\right)\in\left(0,\displaystyle\frac{\pi}{2}\right)$$, the above equation becomes:   cos⁡(nθ(ω))+isin⁡(nθ(ω))=Q(iω)|Q(iω)|. Taking the real part in the above equation, we obtain:   cos⁡(nθ(ω))=ℜ(Q(iω))|Q(iω)|, or equivalently   Tn(1|Q(iω)|1/n)=ℜ(Q(iω))|Q(iω)|, where $$T_n$$ denotes the Chebyshev polynomial of the first kind of order $$n$$. Therefore, we have obtained equation (5.8). As inequality $$(I_1)$$ holds, from the proof of the first part of Theorem 4.1, we know that if $$\beta=0$$, all the roots of the characteristic equation (5.7) are in the left half-plane. The number of the roots of the characteristic equation from the left half-plane can change only if a root (or pair of complex conjugated roots) crosses the imaginary axis. From (5.11) and Lemma 5.1 we can easily see that the critical values of $$\beta$$ decrease with respect to $$\omega$$, and hence, the smallest critical value of $$\beta$$ will be the one corresponding to the largest root $$\omega_n$$ of equation (5.8) from the interval $$(0,\omega_0)$$. Therefore, we obtain the smallest critical value $$\beta_n$$ of the bifurcation parameter, given by (5.9), and we deduce that for any $$\beta\in(0,\beta_n)$$ the equilibrium point $$E$$ is asymptotically stable. Let $$z(\beta)$$ denote the root of the characteristic equation (5.7) satisfying $$z(\beta_n)=i\omega_n$$. Based on the characteristic equation (5.7), we obtain:   (βz(β)+1)n=Q(z(β)). Taking the derivative with respect to $$\beta$$, it follows that   n(βz(β)+1)n−1(βz′(β)+z(β))=z′(β)Q′(z(β)). We obtain   z′(β)=nz(β)(βz(β)+1)Q′(z(β))Q(z(β))−nβ, and hence:   z′(βn)=inωn(iβnωn+1)Q′(iωn)Q(iωn)−nβn. Taking the real part, we obtain:   dℜ(z)dβ|β=βn=nωnℑ((iβnωn+1)Q′(iωn)Q(iωn))|(iβnωn+1)Q′(iωn)Q(iωn)−nβn|2=nωn[βnωnℜ(Q′(iωn)Q(iωn))+ℑ(Q′(iωn)Q(iωn))]|(iβnωn+1)Q′(iωn)Q(iωn)−nβn|2. A laborious computation shows that the term $$\beta_n\omega_n\Re\left(\frac{Q'(i\omega_n)}{Q(i\omega_n)}\right)+\Im\left(\frac{Q'(i\omega_n)}{Q(i\omega_n)}\right)$$ is positive, and hence, $$\left.\displaystyle\frac{d\Re( z)}{d\beta}\right|_{\beta=\beta_n}>0$$, implying that the equilibrium point $$E$$ is asymptotically stable if $$\beta\in(0,\beta_n)$$. System (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$ at $$\beta=\beta_n$$. ■ 5.3 A mix of Dirac and Gamma kernels For simplicity, we will assume that the delay kernels $$h_1,h_2,h_{31},h_{32}$$ are either Dirac or Gamma kernels, such that   H2(z)H32(z)=H1(z)H2(z)H31(z)=e−τz(βz+1)n where $$\tau\geq 0$$, $$\beta>0$$ and $$n\in\mathbb{Z}^+\setminus\{0\}$$. The characteristic equation (4.2) becomes:   (z+w1)(z+w2)(z+w3)+[a(z+w1)+b]e−τz(βz+1)n=0, (5.12) or equivalently   eτz=Q(z)(βz+1)n. Choosing $$\tau$$ as bifurcation parameter, we obtain the following result: Theorem 5.3 (Hopf bifurcations for a mix of delay kernels) Assume that inequalities $$(I_1)$$ and $$(\overline{I_2})$$ are satisfied and that $$\beta\in(0,\beta_n)$$, where $$\beta_n$$ is given by (5.9). Let $$\tilde{\omega}_{n}$$ denote the unique real positive root of the equation   |Q(iω)|2=(β2ω2+1)n (5.13) and consider   τ~np=1ω~n[arccos⁡(ℜ(Q(iω~n)(iβω~n+1)n))+2pπ]. (5.14) The equilibrium point $$E$$ is asymptotically stable if and only if $$\tau\in[0,\tilde{\tau}_{n0})$$. At $$\tau=\tilde{\tau}_{np}$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. Proof Equation (5.12) has a pair of complex conjugated solutions $$z=\pm i\omega$$ on the imaginary axis ($$\omega>0$$) if and only if   eiτω=Q(iω)(iβω+1)n. (5.15) Applying the modulus to both sides of equation, we obtain:   |Q(iω)|2=(β2ω2+1)n. Based on Lemma 5.1, the left hand-side of this equation is strictly decreasing, while the right hand-side is strictly increasing on $$(0,\infty)$$, implying that there is a unique solution $$\tilde{\omega}_n\in(0,\omega_0)$$. Taking the real parts of both sides of equation (5.15), we deduce the critical values of the bifurcation parameter $$\tau$$, given by (5.14). For $$\tau=0$$, taking into account that $$\beta\in(0,\beta_n)$$, from Theorem 5.2 we know that the equilibrium point $$E$$ is asymptotically stable. The number of the roots of the characteristic equation from the left half-plane can change only if a root (or pair of complex conjugated roots) crosses the imaginary axis, i.e. at the critical values $$\tilde{\tau}_{np}$$. Let $$z_{np}(\tau)$$ denote the root of the characteristic equation (5.12) satisfying $$z_{np}(\tilde{\tau}_{np})=i\tilde{\omega}_n$$. Based on the characteristic equation (5.12), we obtain:   (βznp(τ)+1)neτznp(τ)=Q(znp(τ)). Taking the derivative with respect to $$\tau$$, it follows that   n(βznp(τ)+1)n−1βznp′(τ)eτznp(τ)+(βznp(τ)+1)n(τznp′(τ)+znp(τ))eτznp(τ)=znp′(τ)Q′(znp(τ)). We obtain   znp(τ)=znp(τ)Q′(znp(τ))Q(znp(τ))−τ−nββznp(τ)+1, and hence:   znp′(τ~np)=iω~nQ′(iω~n)Q(iω~n)−τ−nβiβω~n+1. Taking the real part, we obtain:   dℜ(znp)dτ|τ=τ~np=ω~n[ℑ(Q′(iω~n)Q(iω~n))+nβ2ω~nβ2ω~n2+1]|Q′(iω~n)Q(iω~n)−τ−nβiβω~n+1|2>0. The positivity follows from Lemma 5.1. Thus, the equilibrium point $$E$$ can only be asymptotically stable if and only if $$\tau\in[0,\tilde{\tau}_{n0})$$ and for any $$p\in\mathbb{Z}^+$$, at $$\tau=\tilde{\tau}_{np}$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. ■ 6. A fractional-order model of the HPA axis In this section, the following fractional-order mathematical model of the HPA axis will be briefly discussed:   {cDqx1(t)=f1(x3(t−τ31))−w1x1(t)cDqx2(t)=f2(x3(t−τ32))x1(t−τ1)−w2x2(t)cDqx3(t)=k3x2(t−τ2)−w3x3(t) (6.1) where $$q\in(0,1)$$ and $$\tau_1,\tau_2,\tau_{31},\tau_{32}\geq 0$$ represent discrete time delays. The classical integer order derivative is generalized by the fractional-order Caputo derivative, defined below (Kilbas et al., 2006; Lakshmikantham et al., 2009; Podlubny, 1999). Definition 6.1 For a continuous function $$x(t)$$, with $$x'\in L^1_{loc}(\mathbb{R}^+)$$, the Caputo fractional-order derivative of order $$q\in(0,1)$$ of $$f$$ is defined by   cDqx(t)=1Γ(1−q)∫0t(t−s)−qx′(s)ds. It is obvious that the fractional-order system (6.1) and the integer-order system (2.1) have the same equilibrium state $$E$$. In order to study the stability of the equilibrium state $$E$$ in the framework of system (6.1) without delays, we rely on the linearization theorem recently proved in Li & Ma (2013). This linearization theorem is an analogue of the classical Hartman theorem for nonlinear integer-order dynamical systems. Moreover, for the corresponding linearized fractional-order system, the following stability result holds Matignon (1996): Theorem 6.1 The linear fractional-order autonomous system   cDqx=Axwhere  A∈Rn×n where $$q\in(0,1)$$ is asymptotically stable if and only if   |arg⁡(λ)|>qπ2∀λ∈σ(A), (6.2) where $$\sigma(A)$$ denotes the spectrum of the matrix $$A$$ (i.e. the set of all eigenvalues). The above theorem shows that in the case of linear fractional-order systems, the necessary and sufficient conditions for the asymptotic stability of the equilibrium state are weaker than the corresponding conditions from the classical integer-order case. Therefore, taking into account Theorem 4.1, we can easily obtain the following result: Proposition 6.1 In the non-delayed case ($$\tau_1=\tau_2=\tau_{31}=\tau_{32}=0$$), if inequality $$(I_1)$$ is satisfied, the equilibrium state $$E$$ of the fractional order system (6.1) is locally asymptotically stable. At this point, the general stability and bifurcation theory for nonlinear fractional-order systems with discrete delays is still an open problem, and an active area of research. Because of the lack of theoretical tools, we have to rely on numerical simulations to exemplify the existence of oscillatory solutions of system (6.1), which will be presented in the next section. 7. Numerical results and discussions 7.1 Parameter values For the numerical simulations, the values of the elimination constants $$w_i$$, $$i\in\{1,2,3\}$$ are computed according to the formula $$w_i=\displaystyle\frac{\ln(2)}{T_i}$$, where $$T_i$$ represent the plasma half-life of hormones, and are given by: $$T_1\approx 4$$ min, $$T_2\approx 19.9$$ min, $$T_3\approx 76.4$$ min (Carroll et al., 2007; Vinther et al., 2011). The equilibrium point $$E$$ of the system should be at the mean values of the hormones: $$\bar{x}_1=7.659$$ pg/ml (24-h mean value of CRH), $$\bar{x}_2=21$$ pg/ml (24-h mean value of ACTH) and $$\bar{x}_3=3.055$$ ng/ml (24-h mean value of free CORT) Carroll et al. (2007). Based on equation (3.3), this leads us to the following relationship:   (f1(x⋆)w1,w3x⋆k3,x⋆)=(x¯1,x¯2,x¯3) and hence:   x⋆=x¯3=3.055 ng/ml;k3=w3x⋆x¯2=ln⁡(2)76.4 min⋅3.055 ng/ml21pg/ml=1.31985 min−1;f1(x⋆)=w1x¯1=ln⁡(2)⋅7.659 pg/ml4 min=1.3272pgml⋅min. Moreover, due to the fact that $$x^\star$$ is the fixed point of the function $$\displaystyle\frac{k_3}{w_1w_2w_3}f_1(x)f_2(x)$$ (see (3.2)), it follows that   f2(x⋆)=w1w2w3k3⋅x⋆f1(x⋆)=0.0955 min−1. For the numerical simulations, the feedback functions $$f_1$$ and $$f_2$$ are considered as in equation (2.2). For fixed values of the parameters $$c$$ (given in ng/ml, as $$\bar{x}_3$$) $$\alpha\in[1,7]$$, $$\eta,\mu\in(0,1]$$ (dimensionless), the parameters $$k_1$$ and $$k_2$$ are uniquely determined, based on the numerical values of $$x^\star$$, $$f_1(x^\star)$$ and $$f_2(x^\star)$$ determined above. Hence:   k1=1.32721−η(3.055)αcα+(3.055)αpgml⋅min,k2=0.09551−μ(3.055)αcα+(3.055)α min−1. In the following, we will assume for simplicity that $$\eta=\mu$$. As for the average time delays, we first observe that the time required by CRH to travel from the hypothalamus to the pituitary through the hypophyseal portal blood vessels is extremely short (Bairagi et al., 2008) and therefore, for simplicity, we will assume a mean time delay $$\tau_1=0$$. The human inhibitory time course concerning the negative feedback of CORT on the production of ACTH shows great variability and has been described in the past as anything between 15 and 60 min (Boscaro et al., 1998; Posener et al., 1997). However, more recently, it has been shown that humans show fast HPA negative feedback (Russell et al., 2010), suggesting that both GR (glucocorticoid receptors) and MR (mineralocorticoid receptors) are involved in this mechanism, with GR effecting a rapid nongenomic feedback at the level of the anterior pituitary and MR sensing higher glucocorticoid levels while levels are still rising (Karst et al., 2005). Hence, we consider a mean delay $$\tau_{32}\in(0,60]$$. In Hermus et al. (1984), a 30-min delay has been reported in the positive-feedforward effect of ACTH on plasma CORT level, leading to the assumption that $$\tau_2\in(0,30]$$. 7.2 Dirac kernels Based on the previous observations and (5.3), we choose the following discrete time delays: average time delay accounting for the positive feedback of the hypothalamus on the pituitary: $$\tau_1=0$$; average time delay due to the positive feedback of the pituitary on the adrenal glands: $$\tau_2\leq 30$$ (min); average time delay due to the negative feedback effect of the adrenal glands on the hypothalamus and pituitary, respectively: $$\tau_{31}=\tau_{32}\leq 60$$ (min). Our aim is to observe periodic solutions for sufficiently large values of the bifurcation parameter $$\tau=\tau_2+\tau_{32}\leq 90$$ (min), depending on the choice of the parameters $$\alpha$$, $$\eta=\mu$$ and $$c$$. Inequality $$(\overline{I_2})$$ has to be fulfilled, because it is a necessary condition for the occurrence of bifurcations. Based on Theorem 5.1 and equation (5.5), we can numerically determine the critical value $$\tau_0$$ of the bifurcation parameter corresponding to the occurrence of a Hopf bifurcation and we are looking for critical values which are smaller that $$90$$ (min). In Fig. 2, for different values of $$\alpha\in[2,7]$$, we have represented the regions in the parameter plane $$(c,\eta=\mu)$$ for which inequality $$(\overline{I_2})$$ is satisfied, and the computed critical value $$\tau_0$$ is within $$90$$ min, $$60$$ min, $$30$$ min and $$15$$ min, respectively. Fig. 2. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_0\in(60,90]$$ (orange, the lowest region), $$\tau_0\in(30, 60]$$ (green, the second lowest region), $$\tau_0\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_0\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of Dirac kernels). Fig. 2. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_0\in(60,90]$$ (orange, the lowest region), $$\tau_0\in(30, 60]$$ (green, the second lowest region), $$\tau_0\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_0\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of Dirac kernels). We observe that for $$\alpha\geq 5$$, for suitable choices of the parameters $$c$$ and $$\eta=\mu$$, the critical value $$\tau_0$$ is smaller than $$15$$ min. This means that periodic solutions can be obtained in system (2.1) with discrete delays satisfying $$\tau_0<\tau_2+\tau_{31}=\tau_2+\tau_{32}<15$$ (min), which is in accordance with the fast feedback observed in humans (Russell et al., 2010). For example, when $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml), we compute $$k_1=18.18$$ (pg/(ml$$\cdot$$min)), $$k_2=1.3$$ (min$$^{-1}$$) and the critical value $$\tau_0=11.4732$$ (min). Therefore, periodic solutions can be observed for $$\tau_2+\tau_{31}=\tau_2+\tau_{32}>11.4372$$ (min) (see Fig. 3). Fig. 3. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=5$$ (min), $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=11.4732$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 3. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=5$$ (min), $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=11.4732$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. However, if $$\alpha<5$$, we note that for any combination of parameters $$c$$ and $$\eta=\mu$$ we obtain $$\tau_0>15$$ (min). For example, when $$\alpha=3$$, $$\mu=\eta=0.95$$ and $$c=2$$ (ng/ml), we compute $$k_1=5.14$$ (pg/(ml$$\cdot$$min)), $$k_2=0.36$$ (min$$^{-1}$$) and the critical value $$\tau_0=46.5028$$ (min). Hence, periodic solutions can only be observed for larger discrete delays (see Fig. 4). Fig. 4. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=27$$ (min), $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=46.5028$$ (min) in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Fig. 4. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=27$$ (min), $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=46.5028$$ (min) in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Numerical simulations show that when the bifurcation parameter $$\tau$$ passes through the critical value $$\tau_0$$, oscillations appear due to Hopf bifurcation phenomena. We notice that oscillations corresponding to the case of smaller critical values (fast feedback) have higher amplitudes and higher frequency (over a 24 hour range) than those corresponding to the case of larger critical values (slow feedback). In Vinther et al. (2011), where the minimal model of the HPA axis has been considered with discrete time delays, it has been reported that individual time delays need to exceed 19 min in order to observe oscillating solutions. Our numerical simulations show that, for a suitable choice of parameters, it is possible to obtain oscillations for time delays much smaller than 19 min, corresponding to the case of fast feedback noticed in Russell et al. (2010). 7.3 Gamma kernels For numerical simulations in the case of gamma kernels, assuming that there is no time delay accounting for the positive feedback of the hypothalamus on the pituitary, we consider $$h_1(t)=\delta(t)$$. Moreover, the other kernels are chosen to be strong gamma kernels (i.e. $$n_2=n_{31}=n_{32}=2$$):   h2(t)=h31(t)=h32(t)=te−t/ββ2. Strong kernels represent the case when the maximum influence on the growth rates of CRH, ACTH and CORT concentrations at any time $$t$$ is due to hormone concentrations at the previous time $$t-\tau$$, where $$\tau=2\beta$$ is the average time-delay. On the other hand, weak kernels (with $$n_2=n_{31}=n_{32}=1$$) would indicate that the maximum weighted response of the growth rates is affected by the current hormone concentration level, while past concentrations have exponentially decreasing influence, which is less plausible than the case of strong kernels, from biological point of view. In this case, $$n=n_2+n_{32}=n_1+n_2+n_{31}=4$$, and hence, the total average time-delay of the system is $$\tau=4\beta$$ (see equation (5.1)). Based on Theorem 5.2 and equation (5.9), we can numerically determine the critical value $$\beta_4$$ of the bifurcation parameter corresponding to the occurrence of a Hopf bifurcation. Therefore, we can compute the critical value of the total average time-delay $$\tau_g=4\beta_4$$. Our aim is to find critical values satisfying $$\tau_g\leq 90$$ (min). In Fig. 5, for different values of $$\alpha\in[2,7]$$, we have represented the regions in the parameter plane $$(c,\eta=\mu)$$ for which inequality $$(\overline{I_2})$$ is satisfied, and the computed critical value $$\tau_g$$ is within $$90$$ min, $$60$$ min, $$30$$ min and $$15$$ min, respectively. Fig. 5. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_g\in(60,90]$$ (orange, the lowest region), $$\tau_g\in(30, 60]$$ (green, the second lowest region), $$\tau_g\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_g\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of strong Gamma kernels). Fig. 5. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_g\in(60,90]$$ (orange, the lowest region), $$\tau_g\in(30, 60]$$ (green, the second lowest region), $$\tau_g\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_g\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of strong Gamma kernels). We observe that for $$\alpha\geq 6$$, for suitable choices of the parameters $$c$$ and $$\eta=\mu$$, the critical value $$\tau_g$$ is smaller than $$15$$ min. This means that periodic solutions can be obtained in system (2.1) with gamma delay kernels satisfying $$\tau=4\beta<15$$ (min), which is in accordance with the fast feedback observed in humans (Russell et al., 2010). When $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml), providing $$k_1=18.18$$ (pg/(ml$$\cdot$$min)) and $$k_2=1.3$$ (min$$^{-1}$$), we compute the critical value of the bifurcation parameter $$\beta_4=3.084$$ (min) and the critical value of the total average time-delay $$\tau_g=12.336$$ (min). In Fig. 6, periodic solutions can be observed for $$\beta=3.5$$ (min). Fig. 6. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=3.5$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=3.084$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 6. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=3.5$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=3.084$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. However, if $$\alpha<6$$, we note that for any combination of parameters $$c$$ and $$\eta=\mu$$ we obtain $$\tau_g>15$$ (min). When $$\alpha=3$$, $$\mu=\eta=0.95$$ and $$c=2$$ (ng/ml), providing $$k_1=5.14$$ (pg/(ml$$\cdot$$min)) and $$k_2=0.36$$ (min$$^{-1}$$), we determine the critical value $$\beta_4=16.9753$$ (min) and hence, $$\tau_g=67.9$$ (min). Periodic solutions can only be seen for $$\beta>16.9753$$ (min) (see Fig. 7). Fig. 7. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=17$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=34$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=16.9753$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Fig. 7. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=17$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=34$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=16.9753$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Similarly as in the case of Dirac kernels, we notice that oscillations corresponding to smaller critical values (fast feedback) have higher amplitudes and higher frequency (over a 24 hour range) than those corresponding to larger critical values (slow feedback). 7.4 Mixed kernels For numerical simulations in the case of mixed kernels, we choose: $$h_1(t)=\delta(t)$$ - no time-delay; $$h_2(t)=\delta(t-\tau_2)$$ - Dirac kernel; $$h_{31}(t)=h_{32}(t)=\displaystyle\frac{t e^{-t/\beta}}{\beta^{2}}$$ strong Gamma-kernels with $$n=n_{31}=n_{32}=2$$. From Theorem 5.3 and equation (5.14), we can numerically determine the critical value $$\tilde{\tau}_{20}$$ of the average time delay due the positive feedback of the pituitary on the adrenal glands $$\tau=\tau_2$$, representing the Hopf bifurcation parameter. When $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml) and $$\beta=3.5$$, with $$k_1=18.18$$ (pg/(ml$$\cdot$$min)) and $$k_2=1.3$$ (min$$^{-1}$$), the critical value is $$\tilde{\tau}_{20}=5.042$$ (min). In Fig. 8, periodic solutions are displayed for $$\tau_2=6$$ (min). Fig. 8. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=6$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=3.5$$ and mean delays $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=5.042$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 8. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=6$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=3.5$$ and mean delays $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=5.042$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. On the other hand, when $$\alpha=3$$, $$\mu=\eta=0.95$$ and $$c=2$$ (ng/ml), with $$k_1=5.14$$ (pg/(ml$$\cdot$$min)) and $$k_2=0.36$$ (min$$^{-1}$$), the critical value is $$\tilde{\tau}_{20}=22.13$$ (min). In Fig. 9, periodic solutions are shown for $$\tau_2=23$$ (min). Fig. 9. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=23$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=20$$ and mean delays $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=22.13$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Fig. 9. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=23$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=20$$ and mean delays $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=22.13$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. As in the previous two cases, oscillations corresponding to smaller critical values (fast feedback) have higher amplitudes and higher frequency (over a 24 h range) than those corresponding to larger critical values (slow feedback). 7.5 Fractional order model The numerical simulations for the fractional-order system (6.1) have been performed using an extension of the Adams-Bashforth-Moulton predictor-corrector method presented in Diethelm et al. (2002). The delay accounting for the positive feedback of the hypothalamus on the pituitary is $$\tau_1=0$$. The delays due to the positive feedback of the pituitary on the adrenal glands and to the negative feedback effect of the adrenal glands on the hypothalamus and pituitary, respectively, are chosen to be equal: $$\tau_2=\tau_{31}=\tau_{32}$$. When $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml), with $$k_1=18.18$$ (pg/(ml$$\cdot$$min)) and $$k_2=1.3$$ (min$$^{-1}$$), a stable limit cycle has been observed numerically for $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min) and the fractional order $$q=0.9$$ (see Fig. 10). As it can been expected, if a smaller fractional order is taken into account (e.g. $$q=0.8$$) for the same delays, the equilibrium point $$E$$ is asymptotically stable (see Fig. 11). Smaller fractional orders are associated with a more pronounced asymptotically stable behaviour of the system. Fig. 10. View largeDownload slide Stable periodic orbit of system (6.1) with fractional order $$q=0.9$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 10. View largeDownload slide Stable periodic orbit of system (6.1) with fractional order $$q=0.9$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 11. View largeDownload slide Trajectories of system (2.1) with fractional order $$q=0.8$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), converge to the asymptotically stable equilibrium point $$E$$, in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 11. View largeDownload slide Trajectories of system (2.1) with fractional order $$q=0.8$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), converge to the asymptotically stable equilibrium point $$E$$, in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. 8. Conclusions and future work In this article, we have generalized the existing minimal model of the HPA axis, firstly, by including distributed time delays and secondly, by considering fractional-order derivatives. This approach to the modelling of the biological processes is more realistic because it involves memory properties, taking into account the whole past history of the variables. These models are able to capture the vital mechanisms of the HPA system. The existence of a unique equilibrium point of the considered models has been proved. Considering general delay kernels in the model with distributed delays, delay-independent sufficient conditions for the local asymptotic stability of the unique equilibrium point have been obtained. These findings are useful if one is unable to estimate the time delays in the system. A thorough bifurcation analysis has been undertaken in three cases: Dirac kernels, Gamma kernels, and finally, a mixed choice of Dirac and Gamma kernels. Critical values of the appropriately chosen bifurcation parameters have been found which account for the occurrence of Hopf bifurcations. Studying the criticality of Hopf bifurcations is a laborious mathematical task, which will be addressed in a future paper. Extensive numerical simulations show that when the bifurcation parameters pass through the critical values, periodic solutions appear which reproduce the ultradian rhythm of the HPA axis. It has been observed that oscillations corresponding to smaller critical values of the bifurcation parameters (fast feedback) have higher amplitudes and higher frequency (over a 24 h range) than those generated by larger critical values (slow feedback). For the fractional-order mathematical model of the HPA axis, it has been shown that if no time delays are considered, the unique equilibrium point is asymptotically stable. When discrete time delays are introduced, we rely on numerical simulations to exemplify the existence of oscillatory solutions for sufficiently large subunitary values of the fractional order. Numerical simulations show that, in the presence of discrete time delays, smaller fractional orders are associated with a more pronounced asymptotically stable behaviour of the system, in a neighbourhood of the equilibrium point. Different approaches of the HPA axis including environmental and physiological perturbations (e.g. in the form of white or coloured noise, or time-varying input) which can be modelled by stochastic and/or impulsive terms, will be developed as future research, with the aim of reproducing both the circadian and ultradian rhythms underlying CORT secretion within the HPA system. Acknowledgements The authors are especially grateful to the editor and the referees for helpful comments and suggestions. Funding This work was supported by grants of the Romanian National Authority for Scientific Research and Innovation, CNCS-UEFISCDI, project no. PN-II-ID-PCE-2011-3-0198 and project no. PN-II-RU-TE-2014-4-0270. References Adimy M. & Crauste F. (2003) Global stability of a partial differential equation with distributed delay due to cellular replication. Nonlinear Anal. Theor. Methods Appl. , 54, 1469– 1491. Google Scholar CrossRef Search ADS   Adimy M. Crauste F. & Ruan S. (2005) Stability and Hopf bifurcation in a mathematical model of pluripotent stem cell dynamics. Nonlinear Anal. Real World Appl. , 6, 651– 670. Google Scholar CrossRef Search ADS   Adimy M. Crauste F. Halanay A. Neamţu M. & Opriş D. (2006). Stability of limit cycles in a pluripotent stem cell dynamics model. Chaos, Solitons Fract. , 27, 1091– 1107. Google Scholar CrossRef Search ADS   Andersen M. Vinther F. & Ottesen J. T. (2013) Mathematical modeling of the hypothalamic–pituitary–adrenal gland (HPA) axis, including hippocampal mechanisms. Math. Biosci. , 246, 122– 138. Google Scholar CrossRef Search ADS PubMed  Bairagi N. Chatterjee S. & Chattopadhyay J. (2008) Variability in the secretion of corticotropin-releasing hormone, adrenocorticotropic hormone and cortisol and understandability of the hypothalamic-pituitary-adrenal axis dynamics - a mathematical study based on clinical evidence. Math Med. Biol. , 25, 37– 63. Google Scholar CrossRef Search ADS PubMed  Bernard S. Bélair, J. & Mackey M. C. (2001) Sufficient conditions for stability of linear differential equations with distributed delay. Discrete Continuous Dyn. Syst. Ser. B. , 1, 233– 256. Google Scholar CrossRef Search ADS   Boscaro M. Paoletta A. Scarpa E. Barzon L. Fusaro P. Fallo F. & Sonino N. (1998) Age-related changes in glucocorticoid fast feedback inhibition of adrenocorticotropin in man 1. J. Clin. Endocrinol. Metab. , 83, 1380– 1383. Google Scholar PubMed  Campbell S. A. & Jessop R. (2009) Approximating the stability region for a differential equation with a distributed delay. Math. Model. Nat. Phenom. , 4, 1– 27. Google Scholar CrossRef Search ADS   Carroll B. J., Cassidy F. Naftolowitz D. Tatham N. E., Wilson W. H., Iranmanesh A. Liu P. Y. & Veldhuis J. D. (2007) Pathophysiology of hypercortisolism in depression. Acta Psychiat. Scand. , 115, 90– 103. Google Scholar CrossRef Search ADS   Conrad M. Hubold C. Fischer B. & Peters A. (2009) Modeling the hypothalamus–pituitary–adrenal system: homeostasis by interacting positive and negative feedback. J. Biol. Phys. , 35, 149– 162. Google Scholar CrossRef Search ADS PubMed  Cushing J. M. (2013) Integrodifferential Equations and Delay Models in Population Dynamics , Vol. 20. Heidelberg, Berlin: Springer Science & Business Media. Diekmann O. & Gyllenberg M. (2012) Equations with infinite delay: blending the abstract and the concrete. J. Differ. Equations , 252, 819– 851. Google Scholar CrossRef Search ADS   Diethelm K. Ford N. J. & Freed A. D. (2002) A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. , 29, 3– 22. Google Scholar CrossRef Search ADS   Faria T. & Oliveira J. J. (2008) Local and global stability for Lotka–Volterra systems with distributed delays and instantaneous negative feedbacks. J. Differ. Equations , 244, 1049– 1079. Google Scholar CrossRef Search ADS   Gudmand-Hoeyer J. Timmermann S. & Ottesen J. T. (2014) Patient-specific modeling of the neuroendocrine HPA-axis and its relation to depression: Ultradian and circadian oscillations. Math. Biosci. , 257, 23– 32. Google Scholar CrossRef Search ADS PubMed  Hermus A. R. M. M. Pieters G. F. F. M. Smals A. G. H. Benraad Th. J. & Kloppenborg P. W. C. (1984) Plasma adrenocorticotropin, cortisol, and aldosterone responses to corticotropin-releasing factor: modulatory effect of basal cortisol levels. J. Clin. Endocrinol. Metab. , 58, 187– 191. Google Scholar CrossRef Search ADS PubMed  Jelić S. Ćupić Ž. & Kolar-Anić L. (2005) Mathematical modeling of the hypothalamic–pituitary–adrenal system activity. Math. Biosci. , 197, 173– 187. Google Scholar CrossRef Search ADS PubMed  Jessop R. & Campbell S. A. (2010) Approximating the stability region of a neural network with a general distribution of delays. Neural Networks , 23, 1187– 1201. Google Scholar CrossRef Search ADS PubMed  Karst H. Berger S. Turiault M. Tronche F. Schütz G. & Joëls M. (2005) Mineralocorticoid receptors are indispensable for nongenomic modulation of hippocampal glutamate transmission by corticosterone. Proc. Natl. Acad. Sci. USA. , 102, 19204– 19207. Google Scholar CrossRef Search ADS   Kilbas A. A., Srivastava H. M. & Trujillo J. J. (2006) Theory and Applications of Fractional Differential Equations . Amsterdam: Elsevier. Kyrylov V. Severyanov L. & Vieira A. (2005) Modeling robust oscillatory behavior of the hypothalamic-pituitary-adrenal axis. IEEE Trans. Biomed. Eng.,   52, 1977– 1983. Google Scholar CrossRef Search ADS   Lakshmikantham V. Leela S. & Devi J. Vasundhara. (2009) Theory of Fractional Dynamic Systems . Cambridge, UK: Cambridge Scientific Publishers. Landsberg L. Young J. B., Wilson J. D. & Foster D. W. (1992) Williams Textbook of Endocrinology . New Jersey: Prentice Hall International. Lenbury Y. & Pornsawad P. (2005) A delay-differential equation model of the feedback-controlled hypothalamus–pituitary–adrenal axis in humans. Math. Med. Biol. , 22, 15– 33. Google Scholar CrossRef Search ADS PubMed  Li C. & Ma Y. (2013) Fractional dynamical system and its linearization theorem. Nonlinear Dyn. , 71, 621– 633. Google Scholar CrossRef Search ADS   Markovic V. M., Cupic Z. Vukojevic V. & Kolar-Anic L. (2011) Predictive modeling of the hypothalamic-pituitary-adrenal (HPA) axis response to acute and chronic stress. Endocrine J. , 58, 889– 904. Google Scholar CrossRef Search ADS   Matignon D. (1996) Stability results for fractional differential equations with applications to control processing. Comput. Eng. Syst. Appl. , 2, 963– 968. Murray J. D. (2002) Mathematical Biology I: An Introduction . Interdisciplinary Applied Mathematics, vol. 17. New York, USA: Springer. Özbay H. Bonnet C. & Clairambault J. (2008) Stability analysis of systems with distributed delays and application to hematopoietic cell maturation dynamics. 47th IEEE Conference on Decision and Control , Cancun, Mexico: IEEE, pp. 2050– 2055. Podlubny I. (1999) Fractional Differential Equations . San Diego, California, USA and London, UK: Academic Press. Pornsawad P. (2013) The feedforward-feedback system of the hypothalamus-pituitary-adrenal axis. Advances in Computing, Communications and Informatics (ICACCI), 2013 International Conference on . Mysore, India: IEEE, pp. 1374– 1379. Posener J. A., Schildkraut J. J., Wilfams G. H., & Schatzberg A. F. (1997) Cortisol feedback effects on plasma corticotropin levels in healthy subjects. Psychoneuroendocrinology , 22, 169– 176. Google Scholar CrossRef Search ADS PubMed  Ruan S. & Wolkowicz G. S. K. (1996) Bifurcation analysis of a chemostat model with a distributed delay. J. Math. Anal. Appl. , 204, 786– 812. Google Scholar CrossRef Search ADS   Russell G. M., Henley D. E., Leendertz J. Douthwaite J. A., Wood S. A., Stevens A. Woltersdorf W. W., Peeters B. W. M. M. Ruigt G. S. F. White A. et al.   (2010) Rapid glucocorticoid receptor-mediated inhibition of hypothalamic–pituitary–adrenal ultradian activity in healthy males. J. Neurosci. , 30, 6106– 6115. Google Scholar CrossRef Search ADS PubMed  Savić D. & Jelić S. (2005) A mathematical model of the hypothalamo-pituitary-adrenocortical system and its stability analysis. Chaos, Solitons Fract. , 26, 427– 436. Google Scholar CrossRef Search ADS   Savić D. Jelić S. & Burić N. (2006) Stability of a general delay differential model of the hypothalamo-pituitary-adrenocortical system. Int. J. Bifurc. Chaos , 16, 3079– 3085. Google Scholar CrossRef Search ADS   Swanson L. W. (2000) Cerebral hemisphere regulation of motivated behavior. Brain Res. , 886, 113– 164. Google Scholar CrossRef Search ADS PubMed  Veldhuis J. D., Keenan D. M., & Pincus S. M. (2008) Motivations and methods for analyzing pulsatile hormone secretion. Endocrine Rev. , 29, 823– 864. Google Scholar CrossRef Search ADS   Vinther F. Andersen M. & Ottesen J. T. (2011) The minimal model of the hypothalamic–pituitary–adrenal axis. J. Math. Biol. , 63, 663– 690. Google Scholar CrossRef Search ADS PubMed  Yuan Y. & Bélair J. (2011) Stability and Hopf bifurcation analysis for functional differential equation with distributed delay. SIAM J. Appl. Dyn. Syst. , 10, 551– 581. Google Scholar CrossRef Search ADS   © The authors 2016. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical Medicine And Biology: A Journal Of The Ima Oxford University Press

# Stability and Hopf bifurcation analysis for the hypothalamic-pituitary-adrenal axis model with memory

, Volume 35 (1) – Mar 1, 2018
30 pages

/lp/ou_press/stability-and-hopf-bifurcation-analysis-for-the-hypothalamic-pituitary-cCaQn08CyR
Publisher
Institute of Mathematics and its Applications
ISSN
1477-8599
eISSN
1477-8602
D.O.I.
10.1093/imammb/dqw020
Publisher site
See Article on Publisher Site

### Abstract

Abstract This article generalizes the existing minimal model of the hypothalamic-pituitary-adrenal (HPA) axis in a realistic way, by including memory terms: distributed time delays, on one hand and fractional-order derivatives, on the other hand. The existence of a unique equilibrium point of the mathematical models is proved and a local stability analysis is undertaken for the system with general distributed delays. A thorough bifurcation analysis for the distributed delay model with several types of delay kernels is provided. Numerical simulations are carried out for the distributed delays models and for the fractional-order model with discrete delays, which substantiate the theoretical findings. It is shown that these models are able to capture the vital mechanisms of the HPA system. 1. Introduction The hypothalamuspituitaryadrenal (HPA) axis is a self-regulated dynamic feedback neuroendocrine system i.e. engaged in the rapid response to stressful stimuli and is responsible for the return to homeostasis through complex feedback mechanisms (Conrad et al., 2009; Swanson, 2000). By regulating the plasma levels of corticosteroids secreted from adrenal glands, it also controls many bodily processes, including mood and emotions, digestion, sexuality, the immune system, energy storage and expenditure. The HPA axis is organized into three distinct regions: the hypothalamus, pituitary gland and adrenal gland, with a complex set of direct influences and feedback interactions among the three endocrine glands. These glands work together by producing and secreting, or responding to common hormones including corticotropin-releasing hormone (CRH), corticotropin (ACTH), and cortisol (CORT) Kyrylov et al. (2005). Both physical stressors (e.g. infection, thermal exposure, dehydration) and psychological stressors (e.g. fear, anticipation) activate the hypothalamus to release CRH, which induces the ACTH production in the pituitary. Then, ACTH is transported by the blood to the adrenal cortex of the adrenal gland, where it stimulates the production of CORT, which in turn suppresses the production of both CRH and ACTH (see Fig. 1). Fig. 1. View largeDownload slide A simple schematic representation of the HPA axis with the negative feedback. Fig. 1. View largeDownload slide A simple schematic representation of the HPA axis with the negative feedback. It is important to keep CORT concentration within a certain physiological range. CORT overproduction, which is often due to a pituitary tumour causing high levels of circulating ACTH, leads to Cushings disease. CORT underproduction, which generates Addisons disease, is most frequently the result of adrenal destruction. In the past few decades, mathematical modelling has started to play an increasingly important role in the study of metabolic and endocrine processes, both in physiology and in clinical medicine (Bairagi et al., 2008). The theory of nonlinear dynamical systems has become a promising research tool for studying rhythmicity in humans. Moreover, it is often mandatory to introduce time delays in the mathematical models describing real life phenomena. Analysing a mathematical model of the HPA axis is helpful in understanding the simultaneous feedback mechanisms and revealing the different ways in which a malfunction can occur (Andersen et al., 2013). It is important to emphasize that time delays of up to 60 min (according to Boscaro et al. (1998); Posener et al. (1997)) unavoidably exist in the HPA axis, due to the transportation of the hormones among the three glands. Several mathematical models of the HPA axis have been developed and analysed by Jelić et al. (2005); Kyrylov et al. (2005); Lenbury & Pornsawad (2005); Savić & Jelić (2005); Savić et al. (2006); Bairagi et al. (2008); Conrad et al. (2009); Markovic et al. (2011); Vinther et al. (2011); Andersen et al. (2013); Pornsawad (2013); Gudmand-Hoeyer et al. (2014). Experimental data show the circadian as well as ultradian rhythm of hormone levels (Carroll et al., 2007; Veldhuis et al., 2008), which should be reflected by the mathematical models of the HPA axis, through the existence of oscillatory solutions. The ultradian rhythm is usually considered an inherent behaviour of the HPA axis, whereas the circadian rhythm is regarded as an external input to the axis (Andersen et al., 2013). A frequently studied model of the HPA axis, called ‘minimal model’ (Vinther et al., 2011), consists of a system of three coupled, non-linear differential equations, with the hormones CRH, ACTH and free CORT as variables. While the mathematical model investigated by Kyrylov et al. (2005) included two more differential equations, accounting for corticosteroid-binding globulin (CBG) bound CORT, and albumin bound CORT, (Vinther et al., 2011) pointed out that only free CORT is capable of interacting with the rest of the HPA-axis, which is the reason for considering only three differential equations in the formulation of the minimal model. If time delays are not incorporated in the mathematical model, no oscillatory behaviour has been detected (Savić & Jelić, 2005; Vinther et al., 2011; Andersen et al., 2013). For systems of ordinary differential equations, sufficiently large time delays are often linked with generating oscillatory solutions. A particular case of the minimal model with discrete time-delays and exponential negative-feedback mechanism has been studied in Lenbury & Pornsawad (2005) and later, in Pornsawad (2013), revealing the occurrence of Hopf bifurcations resulting in the emergence of periodic orbits. Other variants of the minimal model with Hill-type feedback functions and discrete time-delays have been considered in Savić & Jelić (2005); Savić et al. (2006), but no oscillations have been reported. More recently, Vinther et al. (2011) have included discrete time-delays in the minimal model and have observed periodic solutions by computer simulations. Their investigations show that time delays of at least 18 min in the feedback mechanisms are needed for generating oscillations. A slightly modified version of the minimal model with discrete time-delays has been analysed in Bairagi et al. (2008), successfully obtaining the desired Hopf bifurcation and thereby, oscillating solutions for sufficiently large time-delays. In this article, we generalize the existing minimal model of the HPA axis in two realistic ways: firstly, by including distributed time delays and secondly, by considering fractional-order derivatives. On one hand, distributed time delays represent the situation where the delays occur in certain ranges of values with some associated probability distributions, taking into account the whole past history of the variables. In many real world applications, distributed time delays are more realistic and more accurate than discrete time delays (Cushing, 2013). Distributed delay models appear in a wide range of applications such as, population biology (Ruan & Wolkowicz, 1996; Faria & Oliveira, 2008), haematopoiesis (Adimy & Crauste, 2003; Adimy et al., 2005, 2006; Özbay et al., 2008), neural networks (Jessop & Campbell, 2010). On the other hand, the main benefit of fractional-order models in comparison with classical integer-order models is that fractional derivatives provide a good tool for the description of memory and hereditary properties of various processes (Podlubny, 1999; Kilbas et al., 2006; Lakshmikantham et al., 2009). Due to the fact that the whole past history of the variables is accounted for in the formulation of both distributed time-delays as well as fractional-order derivatives, these generalizations of the minimal model are able to reflect biological variability in a better way than other approaches. The article is structured as follows. Section 2 provides the mathematical model of the HPA axis, where instead of considering the transportation of different hormones as an instantaneous process, we introduce distributed time delays to account for the time needed by the hormones to travel from source to destination. In Section 3, the existence of a unique equilibrium point of the system is shown. Local stability analysis of the system with general distributed delays is analysed in Section 4. In Section 5, we undertake a bifurcation analysis for the distributed delay model in the case of several types of delay kernels. The fractional order mathematical model of the HPA axis is presented and shortly analysed in Section 6. Numerical simulations are carried out and discussed in Section 7, followed by concluding remarks in Section 8. 2. Mathematical model of HPA with distributed delays In formulating the mathematical model which describes the variation in time of the concentrations of the three hormones CRH, ACTH and CORT, the following sequence of typical events is considered, according to the schema presented in Fig. 1. CRH is secreted from the hypothalamus and released into the portal blood vessel of the hypophyseal stalk, and then transported to the anterior pituitary where it stimulates the secretion of ACTH, with an average time delay $$\tau_1$$. Then, in the cortex of the adrenal glands, ACTH stimulates the secretion of the stress hormone CORT with the average time delay $$\tau_2$$. CORT has a negative feedback effect on the hypothalamus and the pituitary, expressed by two feedback functions $$f_1$$ and $$f_2$$, affecting the synthesis and release of CRH and ACTH, respectively. On one hand, CORT inhibits the secretion of CRH through glucocorticoid receptors (GRs) situated in the hypothalamus Landsberg et al. (1992), with an average time delay $$\tau_{31}$$. On the other hand, CORT also performs a negative feedback on the secretion of ACTH through GRs situated in the pituitary, with an average time delay $$\tau_{32}$$. The hormone concentrations of CRH, ACTH and CORT are depleted through the rate constants $$w_1$$, $$w_2$$ and $$w_3$$, respectively. The mathematical model of the HPA axis studied in this article is based on the minimal model introduced in Vinther et al. (2011). The main improvement is that we consider distributed time delays to account for the transport of the hormones among the glands. Denoting the hormone concentrations, for simplicity, by $$CRH(t)=x_1(t)$$, $$ACTH(t)=x_2(t)$$, $$CORT(t)=x_3(t)$$, the following system of differential equations with distributed delays is considered:   {x˙1(t)=f1(∫−∞tx3(s)h31(t−s)ds)−w1x1(t),x˙2(t)=f2(∫−∞tx3(s)h32(t−s)ds)∫−∞tx1(s)h1(t−s)ds−w2x2(t),x˙3(t)=k3∫−∞tx2(s)h2(t−s)ds−w3x3(t), (2.1) where all the first terms on the right hand side represent production and all the second terms represent depletion of hormones. The constant $$k_3$$ as well as the elimination constants $$w_1,w_2,w_3$$ are positive. The functions $$f_1,f_2:[0,\infty)\to (0,\infty)$$, which represent the negative feedback from CORT on CRH and ACTH, respectively, are assumed to be strictly decreasing, smooth and bounded on $$[0,\infty)$$. In particular, the results presented in this article are also applicable when Hill functions are being used in the expression of the feedback functions (Vinther et al., 2011; Andersen et al., 2013):   f1(u)=k1(1−ηuαcα+uα),f2(u)=k2(1−μuαcα+uα), (2.2) with $$\alpha\geq 1$$, $$k_1,k_2>0$$, $$\eta,\mu\in(0,1)$$, $$c>0$$. It is easy to verify that functions (2.2) satisfy all the properties mentioned above. However, it may be possible to model the negative feedback using different types of functions $$f_1$$ and $$f_2$$. In this article, our aim is to obtain general results which will also be applicable to other choices of negative feedback functions, besides functions (2.2), often used in the literature. In system (2.1), the delay kernels $$h_1,h_2,h_{31},h_{32}:[0,\infty)\to[0,\infty)$$ are probability density functions representing the probability that a particular time delay occurs. They are assumed to be bounded, piecewise continuous and satisfy   ∫0∞h(s)ds=1. (2.3) The average delay of a delay kernel $$h(t)$$ is given by   τ=∫0∞sh(s)ds<∞. Two important classes of delay kernels, which are often used in the literature, are worth mentioning: Dirac kernels: $$h(s)=\delta(s-\tau)$$, where $$\tau\geq 0$$. In this particular case, the distributed delay is reduced to a discrete time delay:   ∫−∞tx(s)h(t−s)ds=∫0∞x(t−s)δ(s−τ)ds=x(t−τ). Gamma kernels: $$h(s)=\displaystyle\frac{ s^{p-1}e^{-s/\beta}}{\beta^p\Gamma(p)}$$, where $$p,\beta>0$$. The average delay of a Gamma kernel is $$\tau=p\beta$$. The analysis of the mathematical models including particular classes of delay kernels (such as weak Gamma kernels with $$p=1$$ or strong Gamma kernels with $$p=2$$) may shed a light on how distributed delays affect the dynamics differently from discrete delays. However, in the modelling of real world phenomena, one usually does not have access to the exact distribution, and approaches using general kernels may be more useful (Bernard et al., 2001; Campbell & Jessop, 2009; Yuan & Bélair, 2011; Diekmann & Gyllenberg, 2012). Initial conditions associated with system (2.1) are as follows:   xi(s)=φi(s),∀s∈(−∞,0],i=1,2,3, where $$\varphi_i$$ are bounded continuous functions defined on $$(-\infty,0]$$, with values in $$[0,\infty)$$. 3. Existence of a unique equilibrium point An equilibrium point of system (2.1) is a solution of the following algebraic system:   {f1(x3)=w1x1,f2(x3)x1=w2x2,k3x2=w3x3, (3.1) which is equivalent to   {x1=f1(x3)w1,x2=w3k3x3,x3=k3w1w2w3f1(x3)f2(x3). (3.2) Due to the properties of $$f_1$$ and $$f_2$$, the function $$\displaystyle\frac{k_3}{w_1w_2w_3}f_1(x)f_2(x)$$ which appears in the right hand side of the last equation of system (3.2), is strictly positive and strictly decreasing on $$[0,\infty)$$, and therefore, it has a unique fixed point $$x^\star>0$$. It follows that system (2.1) has a unique equilibrium point   E=(f1(x⋆)w1,w3x⋆k3,x⋆). (3.3) In the following, necessary and sufficient conditions will be explored for the local asymptotic stability of the equilibrium point $$E$$ and the occurrence of limit cycles in a neighbourhood of $$E$$ (due to Hopf bifurcations) that can explain the ultradian rhythm. 4. Local stability analysis of system (2.1) In this section, considering general delay kernels, we seek to obtain delay independent sufficient conditions for the local asymptotic stability of the equilibrium point $$E$$. Such results prove to be useful if one is unable to accurately estimate the time delays in system (2.1). Using the transformation $$y_1(t)=x_1(t)-\displaystyle\frac{f_1(x^\star)}{w_1}$$, $$y_2(t)=x_2(t)-\displaystyle\frac{w_3x^\star}{k_3}$$ and $$y_3(t)=x_3(t)-x^\star$$, the linearized system of (2.1) at the equilibrium point $$E$$ is:   {y˙1(t)=f1′(x⋆)∫−∞ty3(s)h31(t−s)ds−w1y1(t),y˙2(t)=f2(x⋆)∫−∞ty1(s)h1(t−s)ds+1w1f1(x⋆)f2′(x⋆)∫−∞ty3(s)h32(t−s)ds−w2y2(t),y˙3(t)=k3∫−∞ty2(s)h2(t−s)ds−w3y3(t). (4.1) The associated characteristic equation of the linearized system (4.1) is:   (z+w1)(z+w2)(z+w3)+a(z+w1)H2(z)H32(z)+bH1(z)H2(z)H31(z)=0, (4.2) where $$H_i( z)=\int_{0}^\infty e^{- z s}h_i(s)ds$$ represent the Laplace transforms of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$ and   a=−k3w1f1(x⋆)f2′(x⋆)=−w2w3x⋆f2′(x⋆)f2(x⋆)>0, (4.3)  b=−k3f1′(x⋆)f2(x⋆)=−w1w2w3x⋆f1′(x⋆)f1(x⋆)>0. (4.4) The following inequalities will be useful for the theoretical analysis:   (I1)8f1(x⋆)+x⋆f1′(x⋆)≥0;(I2)1+x⋆f1′(x⋆)f1(x⋆)+x⋆f2′(x⋆)f2(x⋆)>0;(I2¯)1+x⋆f1′(x⋆)f1(x⋆)+x⋆f2′(x⋆)f2(x⋆)≤0. Remark 4.1 Consider the feedback function   f1(u)=k1(1−ηuαcα+uα) with $$k_1>0$$, $$\eta\in(0,1)$$, $$c>0$$, $$\alpha\geq 1$$. A straightforward computation yields:   8f1(u)+uf1′(u)=k18(1−η)u2α+(16−8η−αη)uαcα+8c2α(cα+uα)2. It can be easily seen that if $$\alpha\leq 8$$, we have   8f1(u)+uf1′(u)≥k18(1−η)u2α+16(1−η)uαcα+8c2α(cα+uα)2≥0, and hence, the inequality $$8f_1(u)+uf_1'(u)\geq 0$$ is satisfied for any $$u\geq 0$$ (regardless of the choice of parameters $$k_1$$, $$\eta$$ or $$c$$). Therefore, if $$\alpha\leq 8$$, inequality $$(I_1)$$ holds as well. It is worth mentioning that according to Murray (2002), $$\alpha>7$$ is considered unphysiological. Remark 4.2 Consider the particular case of feedback functions given by (2.2), with $$k_1,k_2>0$$, $$\eta,\mu\in(0,1)$$, $$\alpha\geq 1$$ and with the constant $$c$$ given by   c=x⋆. This assumption comes from the fact that the constant $$c$$ is chosen to be equal to the mean value of free CORT (Vinther et al., 2011), which, in turn, is equal to the last component of the equilibrium point $$E$$. In this case, the term that appears in the left hand-side of inequalities $$(I_2)$$ and $$(\overline{I_2})$$ becomes   1+x⋆f1′(x⋆)f1(x⋆)+x⋆f2′(x⋆)f2(x⋆)=1−α2[η2−η+μ2−μ]. It is important to note that in this case, inequalities $$(I_2)$$ and $$(\overline{I_2})$$ only depend on the parameters $$\alpha,\eta,\mu$$, and do not depend on the choice of the parameters $$k_1,k_2$$. For instance, if $$\eta=\mu=0.5$$, inequality $$(I_2)$$ is equivalent to $$\alpha<3$$. Theorem 4.1 (Local asymptotic stability) (1) In the non-delayed case, if inequality $$(I_1)$$ is satisfied, then the equilibrium point $$E$$ of system (2.1) is locally asymptotically stable. (2) For any delay kernels $$h_i(t)$$, $$i\in\{1,2,31,32\}$$, if inequality $$(I_2)$$ holds, then the equilibrium point $$E$$ of system (2.1) is locally asymptotically stable. Proof 1. In the non-delayed case, the characteristic equation (4.2) becomes:   λ3+c1λ2+c2λ+c3=0, (4.5) where   c1=w1+w2+w3>0,c2=w1w2+w2w3+w1w3+a>0,c3=w1w2w3+aw1+b>0. We can easily compute   c1c2−c3=a(w2+w3)+γ(w1,w2,w3)+2w1w2w3−b, where   γ(w1,w2,w3)=w12w2+w1w22+w12w3+w1w32+w2w32+w22w3. By the inequality of arithmetic and geometric means, we deduce   γ(w1,w2,w3)≥6w1w2w3,for any w1,w2,w3>0, and hence, based on (4.4) and $$(I_1)$$, we obtain:   c1c2−c3≥a(w2+w3)+8w1w2w3−b==a(w2+w3)+w1w2w3(8+x⋆f1′(x⋆)f1(x⋆))>0. By the Routh–Hurwitz stability test, the equilibrium point $$E$$ is asymptotically stable. 2. From (4.3) and (4.4) it can be easily seen that inequality $$(I_2)$$ is equivalent to   a+bw1<w2w3. The characteristic equation (4.2) can be written as   φ(z)=ψ(z), where the functions $$\varphi$$ and $$\psi$$ are given by   φ(z)=−(z+w1)(z+w2)(z+w3),ψ(z)=a(z+w1)H2(z)H32(z)+bH1(z)H2(z)H31(z). These functions are holomorphic in the right half-plane. Let $$z\in \mathbb{C}$$ with $$\Re(z)\geq 0$$. For any $$i\in\{1,2,31,32\}$$, from (2.3) we obtain:   |Hi(z)|=|∫0∞e−zshi(s)ds|≤∫0∞|e−zs|hi(s)ds=∫0∞e−ℜ(z)shi(s)ds≤∫0∞hi(s)ds=1, and hence, we have:   |ψ(z)|≤a|z+w1||H2(z)||H32(z)|+b|H1(z)||H2(z)||H31(z)|≤a|z+w1|+b=|z+w1|(a+b|z+w1|)=|z+w1|(a+b|z|2+w12+2ℜ(z)w1)≤|z+w1|(a+bw1)<|z+w1|w2w3≤|z+w1|[(|z|2+w22+2ℜ(z)w2)(|z|2+w32+2ℜ(z)w3)]1/2=|z+w1||z+w2||z+w3|=|φ(z)|. Hence, the inequality $$|\psi(z)|<|\varphi(z)|$$ holds for any $$z\in \mathbb{C}$$, $$\Re(z)\geq 0$$. Therefore, the characteristic equation $$\varphi(z)=\psi(z)$$ does not have any root in the right half-plane (or the imaginary axis). This means that all the roots of the characteristic equation (4.2) have strictly negative real part, and the equilibrium $$E$$ is asymptotically stable. ■ Corollary 4.1 For any delay kernels $$h_i(t)$$, $$i\in\{1,2,31,32\}$$, if the equilibrium point $$E$$ of system (2.1) is unstable, then inequality $$(\overline{I_2})$$ holds. In other words, inequality $$(\overline{I_2})$$ is a necessary condition for the occurrence of bifurcations in system (2.1). Remark 4.3 If $$f_1$$ is given by (2.2) with $$\alpha\leq 8$$, it follows from Remark 4.1 that inequality $$(I_1)$$ holds. Based on Theorem 4.1, the equilibrium point $$E$$ is asymptotically stable in the non-delayed case. This improves the sufficient condition $$\alpha\leq 7.46$$ presented in Vinther et al. (2011). As $$\alpha>7$$ is considered unphysiological (Murray, 2002), the equilibrium point $$E$$ is locally asymptotically stable for all realistic values of the parameters if no time delays are considered. Moreover, if $$f_1,f_2$$ are given by (2.2) (such as in Remark 4.2), with $$c=x^\star$$, and the following inequality is satisfied   2α>η2−η+μ2−μ, the equilibrium point $$E$$ of system (2.1) is asymptotically stable for any choice of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$ and of the parameters $$k_1,k_2,k_3,w_1,w_2,w_3$$. In the special case $$\alpha=1$$ considered in (Savić & Jelić, 2005; Savić et al., 2006), it is easy to see that the above inequality is fulfilled for any $$\eta,\mu\in(0,1)$$, implying that the equilibrium point $$E$$ is locally asymptotically stable for any choice of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$ and of the parameters $$k_1,k_2,k_3,w_1,w_2,w_3$$. 5. Bifurcation analysis of system (2.1) The bifurcation analysis presented in this section takes into consideration the average time delays of the delay kernels $$h_i$$, $$i\in\{1,2,31,32\}$$. A first observation is that the time required by CRH to travel from the hypothalamus to the pituitary through the hypophyseal portal blood vessels is extremely short (Bairagi et al., 2008) and therefore, in most numerical simulations the average time delay $$\tau_1$$ is considered close to $$0$$. Moreover, $$\tau_{31}$$ and $$\tau_{32}$$ are comparable, as they represent the average time delays due to the negative feedback effect of the adrenal glands on the hypothalamus and pituitary, respectively, which are closely situated. For this reason, in this section, we will assume for simplicity that   H32(z)=H1(z)H31(z), and we denote   H(z)=H2(z)H32(z)=H1(z)H2(z)H31(z). In fact, $$H(z)$$ is the Laplace transform of the convolution of the delay kernels $$h_2$$ and $$h_{32}$$ defined as   h(t)=∫0th2(s)h32(t−s)ds, with the mean   τ=∫0∞sh(s)ds=τ2+τ32, (5.1) where $$\tau_2$$ and $$\tau_{32}$$ represent the average delays of the kernels $$h_2$$ and $$h_{32}$$ respectively. This results from the fact that the probability density function of the sum of two independent random variables is the convolution of their separate probability density functions. Therefore, the characteristic equation (4.2) becomes   (z+w1)(z+w2)(z+w3)+[a(z+w1)+b]H(z)=0, (5.2) or equivalently:   H(z)−1=Q(z), where   Q(z)=−aw1+b+az(z+w1)(z+w2)(z+w3). The function $$Q(z)$$ defined above will play an important role in the bifurcation analysis presented in this section. We summarize its most important properties in the following Lemma. Lemma 5.1 The function   ω↦|Q(iω)|=a2ω2+(aw1+b)2(ω2+w12)(ω2+w22)(ω2+w32) is strictly decreasing on $$[0,\infty)$$ and the equation   |Q(iω)|=1 has a unique positive real root $$\omega_0$$ if and only if inequality $$(\overline{I2})$$ is satisfied. Moreover, the following inequality holds:   ℑ(Q′(iω)Q(iω))>0∀ω>0. Proof We have   |Q(iω)|2=a2(ω2+w22)(ω2+w32)+2aw1+b2(ω2+w12)(ω2+w22)(ω2+w32) and it is easy to see that $$\omega\mapsto|Q(i\omega)|$$ is strictly decreasing on $$[0,\infty)$$, approaching $$0$$ as $$\omega\rightarrow\infty$$. Therefore, the equation $$|Q(i\omega)|=1$$ has a unique solution if and only if $$|Q(0)|>1$$, or $$w_1w_2w_3<aw_1+b$$, which is equivalent to inequality $$(\overline{I2})$$. Moreover, we have:   ddω|Q(iω)|2=ddω[Q(iω)Q(iω)¯]=2ℜ[Q(iω)¯ddωQ(iω)]=2ℜ[iQ(iω)¯Q′(iω)]=−2ℑ[Q(iω)¯Q′(iω)]=−2|Q(iω)|2ℑ(Q′(iω)Q(iω)). As $$\omega\mapsto |Q(i\omega)|^2$$ is strictly decreasing on $$(0,\infty)$$, its derivative is strictly negative, and hence, $$\Im\left(\displaystyle\frac{Q'(i\omega)}{Q(i\omega)}\right)>0$$, for any $$\omega>0$$. ■ Remark 5.1 A simple computation shows that   ℜ[Q(iω)]=aω4+[b(w1+w2+w3)+a(w12−w2w3)]ω2−w1w2w3(aw1+b)(ω2+w12)(ω2+w22)(ω2+w32). This formula will be useful in the framework of the bifurcation results that follow. Due to the high complexity of the problem, we are unable to perform the bifurcation analysis for general kernels $$h_i$$, $$i\in\{1,2,31,32\}$$. Thus, we focus our attention on the following cases: all delay kernels are Dirac kernels; all delay kernels are Gamma kernels; some delay kernels are Dirac kernels while others are Gamma kernels. 5.1 Dirac kernels Consider that all the delay kernels are Dirac kernels: $$h_1(t)=\delta(t-\tau_1)$$, $$h_2(t)=\delta(t-\tau_2)$$, $$h_{31}(t)=\delta(t-\tau_{31})$$, $$h_{32}(t)=\delta(t-\tau_{32})$$, where $$\tau_1,\tau_2,\tau_{31},\tau_{32}\geq 0$$ satisfy the property   τ2+τ32=τ1+τ2+τ31=τ>0. (5.3) In this case, the characteristic equation (5.2) becomes:   (z+w1)(z+w2)(z+w3)+[a(z+w1)+b]e−τz=0, (5.4) or equivalently:   eτz=Q(z). We choose $$\tau$$ as bifurcation parameter. Theorem 5.1 (Hopf bifurcations in the case of Dirac kernels) Assume that inequalities $$(I_1)$$ and $$(\overline{I_2})$$ are satisfied. For any $$p\in\mathbb{Z}^+$$, consider   τp=arccos⁡[ℜ(Q(iω0))]+2pπω0, (5.5) where $$\omega_0>0$$ is given by Lemma 5.1. The equilibrium point $$E$$ is asymptotically stable if any only if $$\tau\in[0,\tau_0)$$. For any $$p\in\mathbb{Z}^+$$, at $$\tau=\tau_p$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. Proof Equation (5.4) has a pair of complex conjugated solutions $$z=\pm i\omega$$ on the imaginary axis ($$\omega>0$$) if and only if   eiωτ=Q(iω). (5.6) Taking the absolute value in (5.6) we obtain that $$|Q(i\omega)|=1$$, and hence, we obtain that $$\omega=\omega_0$$, where $$\omega_0$$ is the unique positive real solution given by Lemma 5.1. From (5.6) we can also deduce   cos⁡(τω0)=ℜ[Q(iω0)], which leads us to the values $$\tau_p$$ given by (5.5). From Theorem 4.1, we know that the equilibrium point $$E$$ is asymptotically stable if $$\tau=0$$ (since inequality $$(I_1)$$ holds). The number of the roots of the characteristic equation from the left half-plane can change only if a root (or pair of complex conjugated roots) crosses the imaginary axis, or more precisely, whenever $$\tau=\tau_p$$, $$p\in\mathbb{Z}^+$$ (in which case, $$\pm i \omega_0$$ are roots of the characteristic equation). Therefore, for any $$\tau\in[0,\tau_0)$$ the equilibrium point $$E$$ is asymptotically stable. Let $$z_p(\tau)$$ denote the root of the characteristic equation (5.4) satisfying $$z_p(\tau_p)=i\omega_0$$. The function $$z_p(\tau)$$ satisfies   eτzp(τ)=Q(zp(τ)). Taking the derivative with respect to $$\tau$$, it follows that   (τzp′(τ)+zp(τ))eτzp(τ)=zp′(τ)Q′(zp(τ)). We obtain   zp′(τ)=zp(τ)Q′(zp(τ))e−τzp(τ)−τ=zp(τ)Q′(zp(τ))Q(zp(τ))−τ, and hence:   zp′(τp)=iω0Q′(iω0)Q(iω0)−τp. Taking the real part and using Lemma 5.1, we obtain:   dℜ(zp)dτ|τ=τp=ω0ℑ(Q′(iω0)Q(iω0))|Q′(iω0)Q(iω0)−τp|2>0. This nondegeneracy condition for the Hopf bifurcation shows that the equilibrium point $$E$$ can only be asymptotically stable if and only if $$\tau\in[0,\tau_0)$$ and that for any $$p\in\mathbb{Z}^+$$, at $$\tau=\tau_p$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. ■ 5.2 Gamma kernels We now consider that the delay kernels are Gamma distribution kernels: $$h_1(t)=\displaystyle\frac{t^{n_1-1}e^{-t/\beta}}{\beta^{n_1}(n_1-1)!}$$, $$h_2(t)=\displaystyle\frac{t^{n_2-1}e^{-t/\beta}}{\beta^{n_2}(n_2-1)!}$$, $$h_{31}(t)=\displaystyle\frac{t^{n_{31}-1}e^{- t/\beta}}{\beta^{n_{31}} (n_{31}-1)!}$$, $$h_{32}(t)=\displaystyle\frac{ t^{n_{32}-1}e^{-t/\beta}}{\beta^{n_{32}}(n_{32}-1)!}$$, where $$\beta>0$$ and $$n_1,n_2,n_{31},n_{32} \in\mathbb{Z}^+\setminus\{0\}$$ satisfy:   n2+n32=n1+n2+n31=n≥2. The characteristic equation (4.2) becomes:   (z+w1)(z+w2)(z+w3)+a(z+w1)+b(βz+1)n=0, (5.7) or equivalently   (βz+1)n=Q(z). Choosing $$\beta$$ as bifurcation parameter, we obtain the following result: Theorem 5.2 (Hopf bifurcations in the case of Gamma kernels) Assume that inequalities $$(I_1)$$ and $$(\overline{I_2})$$ are satisfied. Let $$\omega_n$$ denote the largest real root of the equation   Tn(1|Q(iω)|1/n)=ℜ(Q(iω))|Q(iω)| (5.8) from the interval $$(0,\omega_0)$$, where $$T_n$$ is the Chebyshev polynomial of the first kind of order $$n$$, and consider   βn=1ωn|Q(iωn)|2/n−1. (5.9) The equilibrium point $$E$$ is asymptotically stable if $$\beta\in(0,\beta_n)$$. At $$\beta=\beta_n$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. Proof Equation (5.7) has a pair of complex conjugated solutions $$z=\pm i\omega$$ on the imaginary axis ($$\omega>0$$) if and only if   (iβω+1)n=Q(iω). (5.10) Applying the modulus to both sides of equation, we obtain:   (β2ω2+1)n=|Q(iω)|2, which means that $$|Q(i\omega)|>1$$. Based on Lemma 5.1, we deduce that $$\omega\in(0,\omega_0)$$ and:   β=1ω|Q(iω)|2/n−1. (5.11) Replacing (5.11) in (5.10) we can also deduce   (1|Q(iω)|1/n+i1−1|Q(iω)|2/n)n=Q(iω)|Q(iω)|. Denoting $$\theta(\omega)=\arccos\left(\displaystyle\frac{1}{|Q(i\omega)|^{1/n}}\right)\in\left(0,\displaystyle\frac{\pi}{2}\right)$$, the above equation becomes:   cos⁡(nθ(ω))+isin⁡(nθ(ω))=Q(iω)|Q(iω)|. Taking the real part in the above equation, we obtain:   cos⁡(nθ(ω))=ℜ(Q(iω))|Q(iω)|, or equivalently   Tn(1|Q(iω)|1/n)=ℜ(Q(iω))|Q(iω)|, where $$T_n$$ denotes the Chebyshev polynomial of the first kind of order $$n$$. Therefore, we have obtained equation (5.8). As inequality $$(I_1)$$ holds, from the proof of the first part of Theorem 4.1, we know that if $$\beta=0$$, all the roots of the characteristic equation (5.7) are in the left half-plane. The number of the roots of the characteristic equation from the left half-plane can change only if a root (or pair of complex conjugated roots) crosses the imaginary axis. From (5.11) and Lemma 5.1 we can easily see that the critical values of $$\beta$$ decrease with respect to $$\omega$$, and hence, the smallest critical value of $$\beta$$ will be the one corresponding to the largest root $$\omega_n$$ of equation (5.8) from the interval $$(0,\omega_0)$$. Therefore, we obtain the smallest critical value $$\beta_n$$ of the bifurcation parameter, given by (5.9), and we deduce that for any $$\beta\in(0,\beta_n)$$ the equilibrium point $$E$$ is asymptotically stable. Let $$z(\beta)$$ denote the root of the characteristic equation (5.7) satisfying $$z(\beta_n)=i\omega_n$$. Based on the characteristic equation (5.7), we obtain:   (βz(β)+1)n=Q(z(β)). Taking the derivative with respect to $$\beta$$, it follows that   n(βz(β)+1)n−1(βz′(β)+z(β))=z′(β)Q′(z(β)). We obtain   z′(β)=nz(β)(βz(β)+1)Q′(z(β))Q(z(β))−nβ, and hence:   z′(βn)=inωn(iβnωn+1)Q′(iωn)Q(iωn)−nβn. Taking the real part, we obtain:   dℜ(z)dβ|β=βn=nωnℑ((iβnωn+1)Q′(iωn)Q(iωn))|(iβnωn+1)Q′(iωn)Q(iωn)−nβn|2=nωn[βnωnℜ(Q′(iωn)Q(iωn))+ℑ(Q′(iωn)Q(iωn))]|(iβnωn+1)Q′(iωn)Q(iωn)−nβn|2. A laborious computation shows that the term $$\beta_n\omega_n\Re\left(\frac{Q'(i\omega_n)}{Q(i\omega_n)}\right)+\Im\left(\frac{Q'(i\omega_n)}{Q(i\omega_n)}\right)$$ is positive, and hence, $$\left.\displaystyle\frac{d\Re( z)}{d\beta}\right|_{\beta=\beta_n}>0$$, implying that the equilibrium point $$E$$ is asymptotically stable if $$\beta\in(0,\beta_n)$$. System (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$ at $$\beta=\beta_n$$. ■ 5.3 A mix of Dirac and Gamma kernels For simplicity, we will assume that the delay kernels $$h_1,h_2,h_{31},h_{32}$$ are either Dirac or Gamma kernels, such that   H2(z)H32(z)=H1(z)H2(z)H31(z)=e−τz(βz+1)n where $$\tau\geq 0$$, $$\beta>0$$ and $$n\in\mathbb{Z}^+\setminus\{0\}$$. The characteristic equation (4.2) becomes:   (z+w1)(z+w2)(z+w3)+[a(z+w1)+b]e−τz(βz+1)n=0, (5.12) or equivalently   eτz=Q(z)(βz+1)n. Choosing $$\tau$$ as bifurcation parameter, we obtain the following result: Theorem 5.3 (Hopf bifurcations for a mix of delay kernels) Assume that inequalities $$(I_1)$$ and $$(\overline{I_2})$$ are satisfied and that $$\beta\in(0,\beta_n)$$, where $$\beta_n$$ is given by (5.9). Let $$\tilde{\omega}_{n}$$ denote the unique real positive root of the equation   |Q(iω)|2=(β2ω2+1)n (5.13) and consider   τ~np=1ω~n[arccos⁡(ℜ(Q(iω~n)(iβω~n+1)n))+2pπ]. (5.14) The equilibrium point $$E$$ is asymptotically stable if and only if $$\tau\in[0,\tilde{\tau}_{n0})$$. At $$\tau=\tilde{\tau}_{np}$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. Proof Equation (5.12) has a pair of complex conjugated solutions $$z=\pm i\omega$$ on the imaginary axis ($$\omega>0$$) if and only if   eiτω=Q(iω)(iβω+1)n. (5.15) Applying the modulus to both sides of equation, we obtain:   |Q(iω)|2=(β2ω2+1)n. Based on Lemma 5.1, the left hand-side of this equation is strictly decreasing, while the right hand-side is strictly increasing on $$(0,\infty)$$, implying that there is a unique solution $$\tilde{\omega}_n\in(0,\omega_0)$$. Taking the real parts of both sides of equation (5.15), we deduce the critical values of the bifurcation parameter $$\tau$$, given by (5.14). For $$\tau=0$$, taking into account that $$\beta\in(0,\beta_n)$$, from Theorem 5.2 we know that the equilibrium point $$E$$ is asymptotically stable. The number of the roots of the characteristic equation from the left half-plane can change only if a root (or pair of complex conjugated roots) crosses the imaginary axis, i.e. at the critical values $$\tilde{\tau}_{np}$$. Let $$z_{np}(\tau)$$ denote the root of the characteristic equation (5.12) satisfying $$z_{np}(\tilde{\tau}_{np})=i\tilde{\omega}_n$$. Based on the characteristic equation (5.12), we obtain:   (βznp(τ)+1)neτznp(τ)=Q(znp(τ)). Taking the derivative with respect to $$\tau$$, it follows that   n(βznp(τ)+1)n−1βznp′(τ)eτznp(τ)+(βznp(τ)+1)n(τznp′(τ)+znp(τ))eτznp(τ)=znp′(τ)Q′(znp(τ)). We obtain   znp(τ)=znp(τ)Q′(znp(τ))Q(znp(τ))−τ−nββznp(τ)+1, and hence:   znp′(τ~np)=iω~nQ′(iω~n)Q(iω~n)−τ−nβiβω~n+1. Taking the real part, we obtain:   dℜ(znp)dτ|τ=τ~np=ω~n[ℑ(Q′(iω~n)Q(iω~n))+nβ2ω~nβ2ω~n2+1]|Q′(iω~n)Q(iω~n)−τ−nβiβω~n+1|2>0. The positivity follows from Lemma 5.1. Thus, the equilibrium point $$E$$ can only be asymptotically stable if and only if $$\tau\in[0,\tilde{\tau}_{n0})$$ and for any $$p\in\mathbb{Z}^+$$, at $$\tau=\tilde{\tau}_{np}$$, system (2.1) undergoes a Hopf bifurcation at the equilibrium point $$E$$. ■ 6. A fractional-order model of the HPA axis In this section, the following fractional-order mathematical model of the HPA axis will be briefly discussed:   {cDqx1(t)=f1(x3(t−τ31))−w1x1(t)cDqx2(t)=f2(x3(t−τ32))x1(t−τ1)−w2x2(t)cDqx3(t)=k3x2(t−τ2)−w3x3(t) (6.1) where $$q\in(0,1)$$ and $$\tau_1,\tau_2,\tau_{31},\tau_{32}\geq 0$$ represent discrete time delays. The classical integer order derivative is generalized by the fractional-order Caputo derivative, defined below (Kilbas et al., 2006; Lakshmikantham et al., 2009; Podlubny, 1999). Definition 6.1 For a continuous function $$x(t)$$, with $$x'\in L^1_{loc}(\mathbb{R}^+)$$, the Caputo fractional-order derivative of order $$q\in(0,1)$$ of $$f$$ is defined by   cDqx(t)=1Γ(1−q)∫0t(t−s)−qx′(s)ds. It is obvious that the fractional-order system (6.1) and the integer-order system (2.1) have the same equilibrium state $$E$$. In order to study the stability of the equilibrium state $$E$$ in the framework of system (6.1) without delays, we rely on the linearization theorem recently proved in Li & Ma (2013). This linearization theorem is an analogue of the classical Hartman theorem for nonlinear integer-order dynamical systems. Moreover, for the corresponding linearized fractional-order system, the following stability result holds Matignon (1996): Theorem 6.1 The linear fractional-order autonomous system   cDqx=Axwhere  A∈Rn×n where $$q\in(0,1)$$ is asymptotically stable if and only if   |arg⁡(λ)|>qπ2∀λ∈σ(A), (6.2) where $$\sigma(A)$$ denotes the spectrum of the matrix $$A$$ (i.e. the set of all eigenvalues). The above theorem shows that in the case of linear fractional-order systems, the necessary and sufficient conditions for the asymptotic stability of the equilibrium state are weaker than the corresponding conditions from the classical integer-order case. Therefore, taking into account Theorem 4.1, we can easily obtain the following result: Proposition 6.1 In the non-delayed case ($$\tau_1=\tau_2=\tau_{31}=\tau_{32}=0$$), if inequality $$(I_1)$$ is satisfied, the equilibrium state $$E$$ of the fractional order system (6.1) is locally asymptotically stable. At this point, the general stability and bifurcation theory for nonlinear fractional-order systems with discrete delays is still an open problem, and an active area of research. Because of the lack of theoretical tools, we have to rely on numerical simulations to exemplify the existence of oscillatory solutions of system (6.1), which will be presented in the next section. 7. Numerical results and discussions 7.1 Parameter values For the numerical simulations, the values of the elimination constants $$w_i$$, $$i\in\{1,2,3\}$$ are computed according to the formula $$w_i=\displaystyle\frac{\ln(2)}{T_i}$$, where $$T_i$$ represent the plasma half-life of hormones, and are given by: $$T_1\approx 4$$ min, $$T_2\approx 19.9$$ min, $$T_3\approx 76.4$$ min (Carroll et al., 2007; Vinther et al., 2011). The equilibrium point $$E$$ of the system should be at the mean values of the hormones: $$\bar{x}_1=7.659$$ pg/ml (24-h mean value of CRH), $$\bar{x}_2=21$$ pg/ml (24-h mean value of ACTH) and $$\bar{x}_3=3.055$$ ng/ml (24-h mean value of free CORT) Carroll et al. (2007). Based on equation (3.3), this leads us to the following relationship:   (f1(x⋆)w1,w3x⋆k3,x⋆)=(x¯1,x¯2,x¯3) and hence:   x⋆=x¯3=3.055 ng/ml;k3=w3x⋆x¯2=ln⁡(2)76.4 min⋅3.055 ng/ml21pg/ml=1.31985 min−1;f1(x⋆)=w1x¯1=ln⁡(2)⋅7.659 pg/ml4 min=1.3272pgml⋅min. Moreover, due to the fact that $$x^\star$$ is the fixed point of the function $$\displaystyle\frac{k_3}{w_1w_2w_3}f_1(x)f_2(x)$$ (see (3.2)), it follows that   f2(x⋆)=w1w2w3k3⋅x⋆f1(x⋆)=0.0955 min−1. For the numerical simulations, the feedback functions $$f_1$$ and $$f_2$$ are considered as in equation (2.2). For fixed values of the parameters $$c$$ (given in ng/ml, as $$\bar{x}_3$$) $$\alpha\in[1,7]$$, $$\eta,\mu\in(0,1]$$ (dimensionless), the parameters $$k_1$$ and $$k_2$$ are uniquely determined, based on the numerical values of $$x^\star$$, $$f_1(x^\star)$$ and $$f_2(x^\star)$$ determined above. Hence:   k1=1.32721−η(3.055)αcα+(3.055)αpgml⋅min,k2=0.09551−μ(3.055)αcα+(3.055)α min−1. In the following, we will assume for simplicity that $$\eta=\mu$$. As for the average time delays, we first observe that the time required by CRH to travel from the hypothalamus to the pituitary through the hypophyseal portal blood vessels is extremely short (Bairagi et al., 2008) and therefore, for simplicity, we will assume a mean time delay $$\tau_1=0$$. The human inhibitory time course concerning the negative feedback of CORT on the production of ACTH shows great variability and has been described in the past as anything between 15 and 60 min (Boscaro et al., 1998; Posener et al., 1997). However, more recently, it has been shown that humans show fast HPA negative feedback (Russell et al., 2010), suggesting that both GR (glucocorticoid receptors) and MR (mineralocorticoid receptors) are involved in this mechanism, with GR effecting a rapid nongenomic feedback at the level of the anterior pituitary and MR sensing higher glucocorticoid levels while levels are still rising (Karst et al., 2005). Hence, we consider a mean delay $$\tau_{32}\in(0,60]$$. In Hermus et al. (1984), a 30-min delay has been reported in the positive-feedforward effect of ACTH on plasma CORT level, leading to the assumption that $$\tau_2\in(0,30]$$. 7.2 Dirac kernels Based on the previous observations and (5.3), we choose the following discrete time delays: average time delay accounting for the positive feedback of the hypothalamus on the pituitary: $$\tau_1=0$$; average time delay due to the positive feedback of the pituitary on the adrenal glands: $$\tau_2\leq 30$$ (min); average time delay due to the negative feedback effect of the adrenal glands on the hypothalamus and pituitary, respectively: $$\tau_{31}=\tau_{32}\leq 60$$ (min). Our aim is to observe periodic solutions for sufficiently large values of the bifurcation parameter $$\tau=\tau_2+\tau_{32}\leq 90$$ (min), depending on the choice of the parameters $$\alpha$$, $$\eta=\mu$$ and $$c$$. Inequality $$(\overline{I_2})$$ has to be fulfilled, because it is a necessary condition for the occurrence of bifurcations. Based on Theorem 5.1 and equation (5.5), we can numerically determine the critical value $$\tau_0$$ of the bifurcation parameter corresponding to the occurrence of a Hopf bifurcation and we are looking for critical values which are smaller that $$90$$ (min). In Fig. 2, for different values of $$\alpha\in[2,7]$$, we have represented the regions in the parameter plane $$(c,\eta=\mu)$$ for which inequality $$(\overline{I_2})$$ is satisfied, and the computed critical value $$\tau_0$$ is within $$90$$ min, $$60$$ min, $$30$$ min and $$15$$ min, respectively. Fig. 2. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_0\in(60,90]$$ (orange, the lowest region), $$\tau_0\in(30, 60]$$ (green, the second lowest region), $$\tau_0\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_0\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of Dirac kernels). Fig. 2. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_0\in(60,90]$$ (orange, the lowest region), $$\tau_0\in(30, 60]$$ (green, the second lowest region), $$\tau_0\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_0\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of Dirac kernels). We observe that for $$\alpha\geq 5$$, for suitable choices of the parameters $$c$$ and $$\eta=\mu$$, the critical value $$\tau_0$$ is smaller than $$15$$ min. This means that periodic solutions can be obtained in system (2.1) with discrete delays satisfying $$\tau_0<\tau_2+\tau_{31}=\tau_2+\tau_{32}<15$$ (min), which is in accordance with the fast feedback observed in humans (Russell et al., 2010). For example, when $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml), we compute $$k_1=18.18$$ (pg/(ml$$\cdot$$min)), $$k_2=1.3$$ (min$$^{-1}$$) and the critical value $$\tau_0=11.4732$$ (min). Therefore, periodic solutions can be observed for $$\tau_2+\tau_{31}=\tau_2+\tau_{32}>11.4372$$ (min) (see Fig. 3). Fig. 3. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=5$$ (min), $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=11.4732$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 3. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=5$$ (min), $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=11.4732$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. However, if $$\alpha<5$$, we note that for any combination of parameters $$c$$ and $$\eta=\mu$$ we obtain $$\tau_0>15$$ (min). For example, when $$\alpha=3$$, $$\mu=\eta=0.95$$ and $$c=2$$ (ng/ml), we compute $$k_1=5.14$$ (pg/(ml$$\cdot$$min)), $$k_2=0.36$$ (min$$^{-1}$$) and the critical value $$\tau_0=46.5028$$ (min). Hence, periodic solutions can only be observed for larger discrete delays (see Fig. 4). Fig. 4. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=27$$ (min), $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=46.5028$$ (min) in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Fig. 4. View largeDownload slide Stable periodic orbit of system (2.1) with Dirac kernels (discrete delays $$\tau_1=0$$, $$\tau_2=27$$ (min), $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau$$ exceeds the critical value $$\tau_0=46.5028$$ (min) in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Numerical simulations show that when the bifurcation parameter $$\tau$$ passes through the critical value $$\tau_0$$, oscillations appear due to Hopf bifurcation phenomena. We notice that oscillations corresponding to the case of smaller critical values (fast feedback) have higher amplitudes and higher frequency (over a 24 hour range) than those corresponding to the case of larger critical values (slow feedback). In Vinther et al. (2011), where the minimal model of the HPA axis has been considered with discrete time delays, it has been reported that individual time delays need to exceed 19 min in order to observe oscillating solutions. Our numerical simulations show that, for a suitable choice of parameters, it is possible to obtain oscillations for time delays much smaller than 19 min, corresponding to the case of fast feedback noticed in Russell et al. (2010). 7.3 Gamma kernels For numerical simulations in the case of gamma kernels, assuming that there is no time delay accounting for the positive feedback of the hypothalamus on the pituitary, we consider $$h_1(t)=\delta(t)$$. Moreover, the other kernels are chosen to be strong gamma kernels (i.e. $$n_2=n_{31}=n_{32}=2$$):   h2(t)=h31(t)=h32(t)=te−t/ββ2. Strong kernels represent the case when the maximum influence on the growth rates of CRH, ACTH and CORT concentrations at any time $$t$$ is due to hormone concentrations at the previous time $$t-\tau$$, where $$\tau=2\beta$$ is the average time-delay. On the other hand, weak kernels (with $$n_2=n_{31}=n_{32}=1$$) would indicate that the maximum weighted response of the growth rates is affected by the current hormone concentration level, while past concentrations have exponentially decreasing influence, which is less plausible than the case of strong kernels, from biological point of view. In this case, $$n=n_2+n_{32}=n_1+n_2+n_{31}=4$$, and hence, the total average time-delay of the system is $$\tau=4\beta$$ (see equation (5.1)). Based on Theorem 5.2 and equation (5.9), we can numerically determine the critical value $$\beta_4$$ of the bifurcation parameter corresponding to the occurrence of a Hopf bifurcation. Therefore, we can compute the critical value of the total average time-delay $$\tau_g=4\beta_4$$. Our aim is to find critical values satisfying $$\tau_g\leq 90$$ (min). In Fig. 5, for different values of $$\alpha\in[2,7]$$, we have represented the regions in the parameter plane $$(c,\eta=\mu)$$ for which inequality $$(\overline{I_2})$$ is satisfied, and the computed critical value $$\tau_g$$ is within $$90$$ min, $$60$$ min, $$30$$ min and $$15$$ min, respectively. Fig. 5. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_g\in(60,90]$$ (orange, the lowest region), $$\tau_g\in(30, 60]$$ (green, the second lowest region), $$\tau_g\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_g\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of strong Gamma kernels). Fig. 5. View largeDownload slide Regions in the parameter plane $$(c,\eta=\mu)$$ resulting in critical values $$\tau_g\in(60,90]$$ (orange, the lowest region), $$\tau_g\in(30, 60]$$ (green, the second lowest region), $$\tau_g\in(15,30]$$ (blue, the third lowest region, if available) and $$\tau_g\in(0,15]$$ (pink, the fourth lowest region, if available) respectively, while inequality $$(\overline{I_2})$$ is fulfilled, for different values of $$\alpha$$. (case of strong Gamma kernels). We observe that for $$\alpha\geq 6$$, for suitable choices of the parameters $$c$$ and $$\eta=\mu$$, the critical value $$\tau_g$$ is smaller than $$15$$ min. This means that periodic solutions can be obtained in system (2.1) with gamma delay kernels satisfying $$\tau=4\beta<15$$ (min), which is in accordance with the fast feedback observed in humans (Russell et al., 2010). When $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml), providing $$k_1=18.18$$ (pg/(ml$$\cdot$$min)) and $$k_2=1.3$$ (min$$^{-1}$$), we compute the critical value of the bifurcation parameter $$\beta_4=3.084$$ (min) and the critical value of the total average time-delay $$\tau_g=12.336$$ (min). In Fig. 6, periodic solutions can be observed for $$\beta=3.5$$ (min). Fig. 6. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=3.5$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=3.084$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 6. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=3.5$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=3.084$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. However, if $$\alpha<6$$, we note that for any combination of parameters $$c$$ and $$\eta=\mu$$ we obtain $$\tau_g>15$$ (min). When $$\alpha=3$$, $$\mu=\eta=0.95$$ and $$c=2$$ (ng/ml), providing $$k_1=5.14$$ (pg/(ml$$\cdot$$min)) and $$k_2=0.36$$ (min$$^{-1}$$), we determine the critical value $$\beta_4=16.9753$$ (min) and hence, $$\tau_g=67.9$$ (min). Periodic solutions can only be seen for $$\beta>16.9753$$ (min) (see Fig. 7). Fig. 7. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=17$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=34$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=16.9753$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Fig. 7. View largeDownload slide Stable periodic orbit of system (2.1) with strong Gamma kernels ($$n_1=0$$, $$n_2=n_{31}=n_{32}=2$$, $$\beta=17$$ (min), mean delays: $$\tau_2=\tau_{31}=\tau_{32}=34$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\beta$$ exceeds the critical value $$\beta_4=16.9753$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Similarly as in the case of Dirac kernels, we notice that oscillations corresponding to smaller critical values (fast feedback) have higher amplitudes and higher frequency (over a 24 hour range) than those corresponding to larger critical values (slow feedback). 7.4 Mixed kernels For numerical simulations in the case of mixed kernels, we choose: $$h_1(t)=\delta(t)$$ - no time-delay; $$h_2(t)=\delta(t-\tau_2)$$ - Dirac kernel; $$h_{31}(t)=h_{32}(t)=\displaystyle\frac{t e^{-t/\beta}}{\beta^{2}}$$ strong Gamma-kernels with $$n=n_{31}=n_{32}=2$$. From Theorem 5.3 and equation (5.14), we can numerically determine the critical value $$\tilde{\tau}_{20}$$ of the average time delay due the positive feedback of the pituitary on the adrenal glands $$\tau=\tau_2$$, representing the Hopf bifurcation parameter. When $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml) and $$\beta=3.5$$, with $$k_1=18.18$$ (pg/(ml$$\cdot$$min)) and $$k_2=1.3$$ (min$$^{-1}$$), the critical value is $$\tilde{\tau}_{20}=5.042$$ (min). In Fig. 8, periodic solutions are displayed for $$\tau_2=6$$ (min). Fig. 8. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=6$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=3.5$$ and mean delays $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=5.042$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 8. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=6$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=3.5$$ and mean delays $$\tau_{31}=\tau_{32}=7$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=5.042$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. On the other hand, when $$\alpha=3$$, $$\mu=\eta=0.95$$ and $$c=2$$ (ng/ml), with $$k_1=5.14$$ (pg/(ml$$\cdot$$min)) and $$k_2=0.36$$ (min$$^{-1}$$), the critical value is $$\tilde{\tau}_{20}=22.13$$ (min). In Fig. 9, periodic solutions are shown for $$\tau_2=23$$ (min). Fig. 9. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=23$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=20$$ and mean delays $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=22.13$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. Fig. 9. View largeDownload slide Stable periodic orbit of system (2.1) with mixed kernels (Dirac kernel $$h_2(t)=\delta(t-\tau_2)$$, with $$\tau_2=23$$ (min), and strong Gamma kernels $$h_{31}=h_{32}$$ with $$n_{31}=n_{32}=2$$, $$\beta=20$$ and mean delays $$\tau_{31}=\tau_{32}=40$$ (min)) due to the Hopf bifurcation taking place when the bifurcation parameter $$\tau=\tau_2$$ exceeds the critical value $$\tilde{\tau}_{20}=22.13$$ (min), in the case of feedback functions $$f_1(u)=5.14\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$, $$f_2(u)=0.36\left(1-\frac{u^{3}}{2000^{3}+u^{3}}\right)$$. As in the previous two cases, oscillations corresponding to smaller critical values (fast feedback) have higher amplitudes and higher frequency (over a 24 h range) than those corresponding to larger critical values (slow feedback). 7.5 Fractional order model The numerical simulations for the fractional-order system (6.1) have been performed using an extension of the Adams-Bashforth-Moulton predictor-corrector method presented in Diethelm et al. (2002). The delay accounting for the positive feedback of the hypothalamus on the pituitary is $$\tau_1=0$$. The delays due to the positive feedback of the pituitary on the adrenal glands and to the negative feedback effect of the adrenal glands on the hypothalamus and pituitary, respectively, are chosen to be equal: $$\tau_2=\tau_{31}=\tau_{32}$$. When $$\alpha=6$$, $$\mu=\eta=1$$ and $$c=2$$ (ng/ml), with $$k_1=18.18$$ (pg/(ml$$\cdot$$min)) and $$k_2=1.3$$ (min$$^{-1}$$), a stable limit cycle has been observed numerically for $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min) and the fractional order $$q=0.9$$ (see Fig. 10). As it can been expected, if a smaller fractional order is taken into account (e.g. $$q=0.8$$) for the same delays, the equilibrium point $$E$$ is asymptotically stable (see Fig. 11). Smaller fractional orders are associated with a more pronounced asymptotically stable behaviour of the system. Fig. 10. View largeDownload slide Stable periodic orbit of system (6.1) with fractional order $$q=0.9$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 10. View largeDownload slide Stable periodic orbit of system (6.1) with fractional order $$q=0.9$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 11. View largeDownload slide Trajectories of system (2.1) with fractional order $$q=0.8$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), converge to the asymptotically stable equilibrium point $$E$$, in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. Fig. 11. View largeDownload slide Trajectories of system (2.1) with fractional order $$q=0.8$$ and discrete delays $$\tau_1=0$$ and $$\tau_2=\tau_{31}=\tau_{32}=14$$ (min), converge to the asymptotically stable equilibrium point $$E$$, in the case of feedback functions $$f_1(u)=18.18\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$, $$f_2(u)=1.3\left(1-\frac{u^{6}}{2000^{6}+u^{6}}\right)$$. 8. Conclusions and future work In this article, we have generalized the existing minimal model of the HPA axis, firstly, by including distributed time delays and secondly, by considering fractional-order derivatives. This approach to the modelling of the biological processes is more realistic because it involves memory properties, taking into account the whole past history of the variables. These models are able to capture the vital mechanisms of the HPA system. The existence of a unique equilibrium point of the considered models has been proved. Considering general delay kernels in the model with distributed delays, delay-independent sufficient conditions for the local asymptotic stability of the unique equilibrium point have been obtained. These findings are useful if one is unable to estimate the time delays in the system. A thorough bifurcation analysis has been undertaken in three cases: Dirac kernels, Gamma kernels, and finally, a mixed choice of Dirac and Gamma kernels. Critical values of the appropriately chosen bifurcation parameters have been found which account for the occurrence of Hopf bifurcations. Studying the criticality of Hopf bifurcations is a laborious mathematical task, which will be addressed in a future paper. Extensive numerical simulations show that when the bifurcation parameters pass through the critical values, periodic solutions appear which reproduce the ultradian rhythm of the HPA axis. It has been observed that oscillations corresponding to smaller critical values of the bifurcation parameters (fast feedback) have higher amplitudes and higher frequency (over a 24 h range) than those generated by larger critical values (slow feedback). For the fractional-order mathematical model of the HPA axis, it has been shown that if no time delays are considered, the unique equilibrium point is asymptotically stable. When discrete time delays are introduced, we rely on numerical simulations to exemplify the existence of oscillatory solutions for sufficiently large subunitary values of the fractional order. Numerical simulations show that, in the presence of discrete time delays, smaller fractional orders are associated with a more pronounced asymptotically stable behaviour of the system, in a neighbourhood of the equilibrium point. Different approaches of the HPA axis including environmental and physiological perturbations (e.g. in the form of white or coloured noise, or time-varying input) which can be modelled by stochastic and/or impulsive terms, will be developed as future research, with the aim of reproducing both the circadian and ultradian rhythms underlying CORT secretion within the HPA system. Acknowledgements The authors are especially grateful to the editor and the referees for helpful comments and suggestions. Funding This work was supported by grants of the Romanian National Authority for Scientific Research and Innovation, CNCS-UEFISCDI, project no. PN-II-ID-PCE-2011-3-0198 and project no. PN-II-RU-TE-2014-4-0270. References Adimy M. & Crauste F. (2003) Global stability of a partial differential equation with distributed delay due to cellular replication. Nonlinear Anal. Theor. Methods Appl. , 54, 1469– 1491. Google Scholar CrossRef Search ADS   Adimy M. Crauste F. & Ruan S. (2005) Stability and Hopf bifurcation in a mathematical model of pluripotent stem cell dynamics. Nonlinear Anal. Real World Appl. , 6, 651– 670. Google Scholar CrossRef Search ADS   Adimy M. Crauste F. Halanay A. Neamţu M. & Opriş D. (2006). Stability of limit cycles in a pluripotent stem cell dynamics model. Chaos, Solitons Fract. , 27, 1091– 1107. Google Scholar CrossRef Search ADS   Andersen M. Vinther F. & Ottesen J. T. (2013) Mathematical modeling of the hypothalamic–pituitary–adrenal gland (HPA) axis, including hippocampal mechanisms. Math. Biosci. , 246, 122– 138. Google Scholar CrossRef Search ADS PubMed  Bairagi N. Chatterjee S. & Chattopadhyay J. (2008) Variability in the secretion of corticotropin-releasing hormone, adrenocorticotropic hormone and cortisol and understandability of the hypothalamic-pituitary-adrenal axis dynamics - a mathematical study based on clinical evidence. Math Med. Biol. , 25, 37– 63. Google Scholar CrossRef Search ADS PubMed  Bernard S. Bélair, J. & Mackey M. C. (2001) Sufficient conditions for stability of linear differential equations with distributed delay. Discrete Continuous Dyn. Syst. Ser. B. , 1, 233– 256. Google Scholar CrossRef Search ADS   Boscaro M. Paoletta A. Scarpa E. Barzon L. Fusaro P. Fallo F. & Sonino N. (1998) Age-related changes in glucocorticoid fast feedback inhibition of adrenocorticotropin in man 1. J. Clin. Endocrinol. Metab. , 83, 1380– 1383. Google Scholar PubMed  Campbell S. A. & Jessop R. (2009) Approximating the stability region for a differential equation with a distributed delay. Math. Model. Nat. Phenom. , 4, 1– 27. Google Scholar CrossRef Search ADS   Carroll B. J., Cassidy F. Naftolowitz D. Tatham N. E., Wilson W. H., Iranmanesh A. Liu P. Y. & Veldhuis J. D. (2007) Pathophysiology of hypercortisolism in depression. Acta Psychiat. Scand. , 115, 90– 103. Google Scholar CrossRef Search ADS   Conrad M. Hubold C. Fischer B. & Peters A. (2009) Modeling the hypothalamus–pituitary–adrenal system: homeostasis by interacting positive and negative feedback. J. Biol. Phys. , 35, 149– 162. Google Scholar CrossRef Search ADS PubMed  Cushing J. M. (2013) Integrodifferential Equations and Delay Models in Population Dynamics , Vol. 20. Heidelberg, Berlin: Springer Science & Business Media. Diekmann O. & Gyllenberg M. (2012) Equations with infinite delay: blending the abstract and the concrete. J. Differ. Equations , 252, 819– 851. Google Scholar CrossRef Search ADS   Diethelm K. Ford N. J. & Freed A. D. (2002) A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. , 29, 3– 22. Google Scholar CrossRef Search ADS   Faria T. & Oliveira J. J. (2008) Local and global stability for Lotka–Volterra systems with distributed delays and instantaneous negative feedbacks. J. Differ. Equations , 244, 1049– 1079. Google Scholar CrossRef Search ADS   Gudmand-Hoeyer J. Timmermann S. & Ottesen J. T. (2014) Patient-specific modeling of the neuroendocrine HPA-axis and its relation to depression: Ultradian and circadian oscillations. Math. Biosci. , 257, 23– 32. Google Scholar CrossRef Search ADS PubMed  Hermus A. R. M. M. Pieters G. F. F. M. Smals A. G. H. Benraad Th. J. & Kloppenborg P. W. C. (1984) Plasma adrenocorticotropin, cortisol, and aldosterone responses to corticotropin-releasing factor: modulatory effect of basal cortisol levels. J. Clin. Endocrinol. Metab. , 58, 187– 191. Google Scholar CrossRef Search ADS PubMed  Jelić S. Ćupić Ž. & Kolar-Anić L. (2005) Mathematical modeling of the hypothalamic–pituitary–adrenal system activity. Math. Biosci. , 197, 173– 187. Google Scholar CrossRef Search ADS PubMed  Jessop R. & Campbell S. A. (2010) Approximating the stability region of a neural network with a general distribution of delays. Neural Networks , 23, 1187– 1201. Google Scholar CrossRef Search ADS PubMed  Karst H. Berger S. Turiault M. Tronche F. Schütz G. & Joëls M. (2005) Mineralocorticoid receptors are indispensable for nongenomic modulation of hippocampal glutamate transmission by corticosterone. Proc. Natl. Acad. Sci. USA. , 102, 19204– 19207. Google Scholar CrossRef Search ADS   Kilbas A. A., Srivastava H. M. & Trujillo J. J. (2006) Theory and Applications of Fractional Differential Equations . Amsterdam: Elsevier. Kyrylov V. Severyanov L. & Vieira A. (2005) Modeling robust oscillatory behavior of the hypothalamic-pituitary-adrenal axis. IEEE Trans. Biomed. Eng.,   52, 1977– 1983. Google Scholar CrossRef Search ADS   Lakshmikantham V. Leela S. & Devi J. Vasundhara. (2009) Theory of Fractional Dynamic Systems . Cambridge, UK: Cambridge Scientific Publishers. Landsberg L. Young J. B., Wilson J. D. & Foster D. W. (1992) Williams Textbook of Endocrinology . New Jersey: Prentice Hall International. Lenbury Y. & Pornsawad P. (2005) A delay-differential equation model of the feedback-controlled hypothalamus–pituitary–adrenal axis in humans. Math. Med. Biol. , 22, 15– 33. Google Scholar CrossRef Search ADS PubMed  Li C. & Ma Y. (2013) Fractional dynamical system and its linearization theorem. Nonlinear Dyn. , 71, 621– 633. Google Scholar CrossRef Search ADS   Markovic V. M., Cupic Z. Vukojevic V. & Kolar-Anic L. (2011) Predictive modeling of the hypothalamic-pituitary-adrenal (HPA) axis response to acute and chronic stress. Endocrine J. , 58, 889– 904. Google Scholar CrossRef Search ADS   Matignon D. (1996) Stability results for fractional differential equations with applications to control processing. Comput. Eng. Syst. Appl. , 2, 963– 968. Murray J. D. (2002) Mathematical Biology I: An Introduction . Interdisciplinary Applied Mathematics, vol. 17. New York, USA: Springer. Özbay H. Bonnet C. & Clairambault J. (2008) Stability analysis of systems with distributed delays and application to hematopoietic cell maturation dynamics. 47th IEEE Conference on Decision and Control , Cancun, Mexico: IEEE, pp. 2050– 2055. Podlubny I. (1999) Fractional Differential Equations . San Diego, California, USA and London, UK: Academic Press. Pornsawad P. (2013) The feedforward-feedback system of the hypothalamus-pituitary-adrenal axis. Advances in Computing, Communications and Informatics (ICACCI), 2013 International Conference on . Mysore, India: IEEE, pp. 1374– 1379. Posener J. A., Schildkraut J. J., Wilfams G. H., & Schatzberg A. F. (1997) Cortisol feedback effects on plasma corticotropin levels in healthy subjects. Psychoneuroendocrinology , 22, 169– 176. Google Scholar CrossRef Search ADS PubMed  Ruan S. & Wolkowicz G. S. K. (1996) Bifurcation analysis of a chemostat model with a distributed delay. J. Math. Anal. Appl. , 204, 786– 812. Google Scholar CrossRef Search ADS   Russell G. M., Henley D. E., Leendertz J. Douthwaite J. A., Wood S. A., Stevens A. Woltersdorf W. W., Peeters B. W. M. M. Ruigt G. S. F. White A. et al.   (2010) Rapid glucocorticoid receptor-mediated inhibition of hypothalamic–pituitary–adrenal ultradian activity in healthy males. J. Neurosci. , 30, 6106– 6115. Google Scholar CrossRef Search ADS PubMed  Savić D. & Jelić S. (2005) A mathematical model of the hypothalamo-pituitary-adrenocortical system and its stability analysis. Chaos, Solitons Fract. , 26, 427– 436. Google Scholar CrossRef Search ADS   Savić D. Jelić S. & Burić N. (2006) Stability of a general delay differential model of the hypothalamo-pituitary-adrenocortical system. Int. J. Bifurc. Chaos , 16, 3079– 3085. Google Scholar CrossRef Search ADS   Swanson L. W. (2000) Cerebral hemisphere regulation of motivated behavior. Brain Res. , 886, 113– 164. Google Scholar CrossRef Search ADS PubMed  Veldhuis J. D., Keenan D. M., & Pincus S. M. (2008) Motivations and methods for analyzing pulsatile hormone secretion. Endocrine Rev. , 29, 823– 864. Google Scholar CrossRef Search ADS   Vinther F. Andersen M. & Ottesen J. T. (2011) The minimal model of the hypothalamic–pituitary–adrenal axis. J. Math. Biol. , 63, 663– 690. Google Scholar CrossRef Search ADS PubMed  Yuan Y. & Bélair J. (2011) Stability and Hopf bifurcation analysis for functional differential equation with distributed delay. SIAM J. Appl. Dyn. Syst. , 10, 551– 581. Google Scholar CrossRef Search ADS   © The authors 2016. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

### Journal

Mathematical Medicine And Biology: A Journal Of The ImaOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations

Abstract access only

18 million full-text articles

Print

20 pages / month

PDF Discount

20% off