{'fraction_non_alphanumeric': 0.034992287082893564, 'fraction_numerical': 0.015425834212876513, 'mean_word_length': 4.793979303857008, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability analysis of a hypothalamic-pituitary-adrenal axis model with inclusion of glucocorticoid receptor and memory

This paper analyzes a four-dimensional model of the hypothalamic-pituitary-adrenal (HPA) axis that includes the influence of the glucocorticoid receptor in the pituitary. Due to the spatial separation between the hypothalamus, pituitary and adrenal glands, distributed time delays are introduced in the mathematical model. The existence of the positive equilibrium point is proved and a local stability and bifurcation analysis is provided, considering several types of delay kernels. The fractional-order model with discrete time delays is also taken into account. Numerical simulations are provided to illustrate the effectiveness of the theoretical findings.


INTRODUCTION
One of the most important systems in stress response is hypothalamus-pituitary-adrenal (HPA) axis, because the released hormones lead to energy directed to the organisms. The HPA axis is organized into three distinct regions: the hypothalamus, pituitary gland and adrenal gland. It is 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), adenocorticotropin (ACTH), and cortisol (CORT) [1]. The dynamics of the HPA axis is very important, because cortisol overproduction leads to Cushings disease; cortisol underproduction generates Addisons disease.
During the past decade, the mathematical modelling of the HPA axis has been intensively studied [1,2,3,4,5,6,7,8,9,10,11,12]. In 2005, Savic and Jelic [10] presented one of the mathematical models of HPA axis described by a system with three nonlinear differential equations with one delay. Three variables have been considered: the concentration of CRH, the concentration of ACT H and the concentration of free cortisol CORT . According to the results obtained by stability analysis, the system does not oscillate at all in the neighborhood of the equilibrium point. In 2006, Savic et al. [11] considered the same model, but four time delays have been introduced. Once more, their results showed that the system describing the HPA axis does not oscillate in a neighborhood of the equilibrium points.
In 2011, Vinther et al. [12] introduced feedback functions as generalized Hill functions, based on the receptor dynamics. They investigated the possibility of time delays as being capable of producing oscillations in accordance with the ultradian rhythm. By numercial simulations, oscillating solutions in a neighborhood of the equilibrium point have been observed.
In 2015, Kaslik and Neamtu [13] proved the occurrence of oscillating solutions in a neighborhood of the unique equilibrium point, by performing a Hopf bifurcation analysis, considering several distributed time delays and fractional derivatives in the minimal three-dimensional mathematical model previously analyzed in [12].
On the other hand, Gupta et al. [14] and Sriram et al. [15] incorporated expression of the glucocorticoid receptor (GR) in the mathematical model of the HPA axis and demonstrated that repeated stress and the level of GR concentration leads to a disorder; however, it is important to note that they did not consider time delays.
With the aim of improving the mathematical model of the HPA axis described in [14], we include distributed time delays and 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 [16]. Distributed delay models appear in a wide range of applications such as, population biology [17,18], hematopoiesis [19,20,21,22], neural networks [23]. 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 [24,25,26].
The paper is structured as follows. Section 2 provides the mathematical model of the HPA axis, with the influence of the GR concentration, where we introduce distributed time delays to account for the time needed by the hormones to travel from source to destination, as well as fractional-order derivatives to account for the memory properties of the system. In Section 3, a local stability analysis of the system with distributed delays is provided. Numerical simulations are carried out and discussed in Section 4, followed by concluding remarks in Section 5.

THE MATHEMATICAL MODEL
The mechanism of the HPA axis can be shortly described as follows. Paraventricular nuclei (PVN) of hypothalamus generate corticotropin releasing hormone (CRH) which induces the adenocorticotropin (ACTH) production in the pituitary. Then, ACTH stimulates the adrenal cortex to produce glucocorticoids, which in turn suppress the production of both CRH and ACTH.
In formulating the mathematical model which describes the variation in time of the concentrations of the three hormones CRH, ACTH and CORT, as well as the glucocorticoid receptors GR, the following sequence of typical events is considered, according to the schema presented in Figure 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. In our mathematical model, the evolution of the CRH hormone concentration in the hypothalamus is described by: where a 1 models a circadian production and stress term, a 2 is the inhibition constant, a 3 is the degradation rate of cortisol.
In the pituitary, cortisol enters the cells and binds the glucocorticoid receptor in the cytoplasm, causing the receptor to dimerize. This leads to the complex GR × CORT that inhibits the production of ACT H, modeled by: where b 1 stands for a production term, b 2 is the inhibition constant, b 3 is the degradation rate of ACT H, anḋ where c 1 is the dimerization constant, c 2 the binding affinity constant, c 3 the production constant and c 4 the degradation term for pituitary GR production.
As for the evolution of cortisol concentration produced by the adrenal, the following equation is considered: where d 1 is a production coefficient and d 2 is the degradation rate of CORT . Therefore, the four-dimensional mathematical model with distributed time delay is: with the initial conditions 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. ( The average delay of a delay kernel h(t) is given by Two important classes of delay kernels often used in the literature are worth mentioning: • Dirac kernels: h(s) = δ(s − τ), where τ ≥ 0. In this particular case, the distributed delay is reduced to a discrete time delay: On the other hand, we also consider the four-dimensional mathematical model with fractional-order derivatives of Caputo type and discrete time delays τ 1 , τ 2 > 0: with the initial conditions

STABILITY OF THE EQUILIBRIUM POINT
In what follows we use the notations: The coordinates of the equilibrium point of systems (1) and (3) are the positive solution of the following algebraic system: In the case of the distributed time-delay system (1), considering the transformation: the linearized system at an equilibrium point E 0 becomes: where The characteristic equation of (4) is given by: If there are no delays, the characteristic equation (5) is given by: Using the Hurwitz criteria for equation (6) we obtain:
In the presence of different types of time delays, we analyze three cases, as follows. Case 1. Dirac kernels: h i (s) = δ(s − τ i ), i = 1, 2.

Proposition 3
Suppose that z 0 = ω 2 0 , and The sign of Re dλ(τ) dτ | τ=τ j is the same as that of dh(z) dz . Therefore, the following transversality condition holds: which corresponds to the occurrence of a Hopf bifurcation in a neighborhood of the equilibrium point E 0 .

Proposition 4
Suppose that z 0 = ω 2 0 , and The sign of Re dλ(τ) dτ | τ=τ 1 j is the same as that of dh(z) dz . Therefore, the following transversality condition holds: which corresponds to the occurrence of a Hopf bifurcation in a neighborhood of the equilibrium point E 0 .

NUMERICAL SIMULATIONS
For the numerical simulations, we consider the scaled parameter values used in [14]: For these values of the system parameters, three positive equilibrium states coexist because of the homodimerization of the GR with cortisol: one equilibrium state with a low GR concentration which is asymptotically stable for any choice of time delays, one equilibrium state with a medium GR concentration which is unstable for any choice of time delays and one equilibrium state with a high GR concentration whose stability depends on the choice of time delays. According to [14], the low GR concentration represents the normal state, and high GR concentration reflects a dysregulated HPA axis resulting in persistently low cortisol levels.
In the following, we investigate the stability of the equilibrium point with high GR concentration: In the case of discrete time-delays: h 1 (s) = δ(s − τ 1 ), h 2 (s) = δ(s − τ 2 ), with τ 1 = 25 min, τ 2 = τ − τ 1 min, we compute the critical value τ = 32.8043 min, according to Proposition 3. Figure 2 displays the trajectories of system (1) for τ = τ 1 + τ 2 = 31 min and τ = τ 1 + τ 2 = 33 min, respectively. When τ = 31 min (below the Hopf bifurcation value), we observe that the trajectories of the system converge to the asymptotically stable equilibrium state E 0 . On the other hand, when τ = 33 min (above the Hopf bifurcation value), the trajectories of the system converge to the asymptotically stable limit cycle which appears in a neighborhood of E 0 due to the supercritical Hopf bifurcation. In the case of weak gamma kernels: h 1 (s) = h 2 (s) = ae −as , with a = τ −1 , the equilibrium state E 0 is asymptotically stable, for any a > 0. Figure 3 shows the trajectories of system (1) for τ = 50 min, where the trajectories of the system converge to the asymptotically stable equilibrium state E 0 .
For the fractional-order model (3) with q = 0.8, we consider the discrete time-delays: h 1 (s) = δ(s − τ 1 ), h 2 (s) = δ(s − τ 2 ), with τ 1 = 25 min, τ 2 = τ − τ 1 min. Figure 4 displays the trajectories of system (3) for τ = τ 1 + τ 2 = 39 min and τ = τ 1 + τ 2 = 40 min, respectively. When τ = 39 min (below the Hopf bifurcation value), we observe that the trajectories of the system converge to the asymptotically stable equilibrium state E 0 . On the other hand, when τ = 40 min (above the Hopf bifurcation value), the trajectories of the system converge to the asymptotically stable limit cycle which appears in a neighborhood of E 0 due to the supercritical Hopf bifurcation. Numerical simulations show that, compared to the integer-order model with discrete time-delays, as the fractional order decreases, the critical value for the Hopf bifurcation parameter τ = τ 1 + τ 2 increases.   In this paper, we have considered a four-dimensional mathematical model that describes the hypothalamus-pituitaryadrenal axis with the influence of the GR concentration. Due to the fact that the involved processes are not instantaneous, we have incorporated distributed delays, and we have also taken into account fractional-order derivatives for the modelling of memory properties. This approach to the modelling of the biological processes is more realistic because it takes into account the whole past history of the variables, capturing the vital mechanisms of the HPA system. Sufficient conditions for the local asymptotic stability of the equilibrium points have been obtained and a Hopf bifurcation analysis has been undertaken. Numerical simulations reflect the importance of the theoretical results. As a direction for future work, this model can be generalized by considering general feedback functions to account for the interactions within the HPA axis.