A SIMPLIFICATION AND ANALYSIS OF THE HPA AXIS MODEL

The human stress response is controlled largely by the hypothalamic pituitary adrenal (HPA) axis. Models predicting the levels of the hormones involved are very often not analytically solvable because of nonlinear complexity. A simplified HPA axis model with the minimum nonlinear terms including the glucocorticoid receptor (GR) is proposed. Quantitative analysis and numerical simulation related oscillating solutions are both investigated in this model.


Introduction
The activity of the hypothalamic pituitary adrenal (HPA) axis involves the consecutive release of corticotrophin releasing hormone (CRH) in the hypothalamus, adrenocorticotropic hormone (ACTH) in the pituitary, and cortisol (O) in the adrenal glands.These are released in sequence in response to stressors, and their levels shortly after a stressor are of interest for medical purposes, specifically in treating depression and PTSD.Additionally, their baseline levels oscillate with a period of approximately 24 hours, another factor of interest for the same reasons.The traditional model of this axis exhibits feedback inhibition in that O down regulates the production of both CRH and ACTH [1].
Gupta et al. [2] proposed a model for the HPA axis that involved an additional variable, the glucocorticoid receptor (R or GR), whose levels are pertinent to those of the other hormones.
In particular, GR binds to cortisol, which causes GR to dimerize.This newly-bound glucocorticoid receptor (OR) then down regulates the production of ACTH and up regulates that of GR.Thus, there is an added step that controls the extent to which cortisol inhibits ACTH production in the pituitary, while cortisol still directly inhibits CRH production in the hypothalamus.The inclusion of this GR creates a bistability in the axis, with both the normal low GR concentration state and the alternate high GR concentration state.The latter is associated with a downstream decrease in the levels of O, a state known as hypocortisolism.Hypocortisolism is often associated with post-traumatic stress disorder (PTSD) [3], as well as a number of inflammatory and autoimmune diseases, including rheumatoid arthritis and certain allergic conditions [4].
A number of newer models have also incorporated a state of hypercortisolism, which is often associated with major depression [5].This is not always the case, however, because individuals with depression tend to have much greater variance in their cortisol levels than do individuals who are healthy or those who have PTSD [6] [7].Xiao [8] presented a theoretical model for HPA axis considering GR as the only negative feedback term for CRH and ACTH.The motivation behind this simplified model is to create one which retains the negative feedback feature of the HPA axis model using the newly-bound OR while providing more analytic insight into equilibrium stability and oscillation solutions related to circadian rhythm of hormone concentrations.
The organization of this paper is as follows.Section 2 introduces the formulation of the model.Section 3 describes the quantitative analysis of this model including existence and stability of equilibriums.Parameter estimations and numerical simulations are given in Section 4. The paper ends with a discussion in Section 5.

Builing the model
In this section, we build the model based on plausible assumptions first, then we apply nondimensionalzation to further simplify the system.

Model Formulation
In developing the simplified model of HPA axis, two assumptions were made: (i) this model takes into account a linear production of GR and the rate of its degradation.As for the other hormones involved, this model makes certain simplifications, with the rate of change of each enzyme essentially depending on either the level of stress (in the case of CRH) or the level of the enzyme that came before it (in the cases of ACTH and CORT); and (ii) the negative feedback to CRH uses the linear inhibition with bound glucocorticoid receptor (OR).The simplified model is formulated as follows: where CRH is the concentration of corticotrophin releasing hormone, O is the concentration of cortisol, R is the concentration of glucocorticoid receptor, K i1 is the negative feedback strength of OR on CRH, K c is the production rate constant of CRH, and K cd is the degradation constant of CRH.Based on the second assumption, we propose that CRH production is regulated by the bounded glucocorticoid receptor OR.ACTH is the concentration of adrenocorticotropic hormone, and it is positively regulated by the concentration of CRH, which stimulates its production; its production is down regulated by its own degradation as well as the presence of OR.Thus, we arrive at second equation in Equation 1. K i2 is the negative feedback strength of OR on ACTH, K a is the production rate constant of ACTH, and K ad is the degradation constant of ACTH.
The concentration of O is generally considered to be determined by the concentration of ACTH (which increases the rate of its production) as well as its own concentration because it expresses negative feedback with itself.Both of these concepts are present in the Sriram model [5], and they lead to the third equation.In this equation, K o is the production rate of O and K od is its degradation rate.According to the aforementioned Gupta model [2], when O binds to GR to create the bounded complex OR, the newly-bounded OR then up-regulates the production of GR.We therefore simplify this mechanism to say that O directly up regulates the production of GR.When we additionally include the degradation of GR, we arrive at the forth equation in this system.In this equation, K r is the linear production term of GR, and K rd is the linear degradation term thereof.

Nondimensionalization
Define new system parameters in order to simplify (1): And thus (1) becomes:

Model Analysis
In this section, we offer the qualitative analysis of this system.

Existence of Equilibrium
The equilibrium of system (2) satisfies the following system of algebra equation where all parameters are positive .
After substitutions in (3), we obtain where the first equation in (4) can be written as so if there exist positive real roots for the above equation, then the system will have positive equilibrium under condition, 1 − o 2 > 0. Set: then we have 2) has a unique positive equilibrium.
Proof.From second equation of ( 4), we can have only if 0 < o < 1 will make the equilibrium positive which stands for the concentration of hormone CRH.
Based on Descartes' Rule of Signs, there are two sign changes in positive -root case: and there are two sign changes in negative-root case: So we have the following conclusion: the corresponding polynomial will have 2 negative real roots, 2 complex roots, or 2 positive real roots, 2 complex roots, or 2 positive real roots and 2 negative real roots, or 4 complex roots.
Since f (0) = A > 0 and f (1) = B < 0, based on Intermediate Value Theorem, we can get that there exists at least one positive real root in the interval (0, 1).Then for the polynomial, there will be two positive real roots for sure,and there exists only one positive real root on the interval (0, 1).Hence we prove that system (2) has a unique positive equilibrium.

Determination of Oscillation Solution
Denote the equilibrium of ( 2) as (c 0 , a 0 , o 0 , r 0 ).The Jacobian matrix of ( 2) is as follows: The associated characteristic equation is given by where It's easy to see that a 1 , a 2 , a 3 , a 4 are all positive since all parameters are positive, so we can conclude there is no non-negative real root for the characteristic equation based on Descartes' Rule of Signs.And the corresponding equation will have 0 real negative roots, 4 complex roots, or 2 real negative roots, 2 complex roots, or 4 negative real roots.From the Routh-Hurwitz criterion, if two conditions, a 1 a 2 − a 3 > 0 and a 1 a 2 a 3 > (a 3 ) 2 + (a 1 ) 2 a 4 , are satisfied, then the equilibrium (c 0 , a 0 , o 0 , r 0 ) in ( 2) is local asymptotically stable.
However, since circadian rhythm is a key feature in HPA axis, we would like to see if periodic solution can appear in this system or not.Liu (1994) [9] presented a criterion for supercritical Hopf Bifurcation related to Routh-Hurwitz criterion and the appearance of stable limit cycles can be interpreted as the circadian rhythm in this system.The Hopf bifurcation occurs when , which for our model is when: and we have a 1 a 2 a 3 = (a 3 ) 2 + (a 1 ) 2 a 4 . .

Numerical Analysis
The derived condition above is very difficult to simplify, requiring numerical analysis under certain parameter values.In this section, we will use XPPAUT, Matlab, and R to implement the numerical analysis.

Bifurcation Diagram
Bifurcation analyses through software XPPAUT [10] is carried out by varying k c1 , the negative feedback term for CRH, as the bifurcation parameter and O as the dynamical variable.
Bifurcation diagram is shown in Figure 1, and Hopf point occurs when k c1 increase.Also, as we can see from the diagram, the periodic behavior of o is only valid for a small regime, since we require the solutions to be positive.

Discussion
We have derived a simplified HPA axis model with glucocorticoid receptor considered.After nondimensionalization for the system, when 0 < o < 1, there exists a unique positive equilibrium which is locally stable under a certain condition.Also, under another condition there exists a Hopf bifurcation point on bifurcation parameter k c1 , where the appearance of limit cycles can be interpreted as the circadian rhythm of the secretion of hormone.However, we can apply the parameter estimation based on the data, it shows that our system only capture the circadian rhythm in one period.After short oscillation the system will go out of bound.Two reasons can cause this situation to occur: 1) The structure of this model is over simplified.Most of the HPA axis models contain highly nonlinear terms such as Michaelis-Menten kinetics or Hill equations to describe the negative feedback and the degradation of those hormones, which lead the oscillation more stable; 2) The mechanism of Hopf bifurcation need more depth exploration.We

FIGURE 1 .
FIGURE 1. Bifurcation analysis was carried out with k c1 as the bifurcation parameter.All the other four parameters were kept constant throughout the simulation.All five parameters are set to 1 to get a stable steady state.As k c1 increases, Hopf bifurcation was found where k c1 is approximately 16.23.

FIGURE 2 .
FIGURE 2. Parameter Estimation for first normal subjects with scaled raw Cortisol data.All the data is scaled between 0 to 1 to fit the restriction of o shown above.The circadian rhythm can be captured during the first 24 hours and then the simulation goes out of bound quickly.

FIGURE 3 .
FIGURE 3. Parameter Estimation for second normal subjects with scaled raw Cortisol data.All the data is scaled between 0 to 1 to fit the restriction of o shown above.The circadian rhythm can be captured during the first 24 hours and then the simulation goes out of bound in a slow pace.

TABLE 1 .
Parameters of Simplified HPA axis model