  A four-dimensional mathematical model of the hypothalamus-pituitary-adrenal
(HPA) axis is investigated, incorporating the influence of the GR concentration
and general feedback functions. The inclusion of distributed time delays
provides a more realistic modeling approach, since the whole past history of
the variables is taken into account. The positivity of the solutions and the
existence of a positively invariant bounded region are proved. It is shown that
the considered four-dimensional system has at least one equilibrium state and a
detailed local stability and Hopf bifurcation analysis is given. Numerical
results reveal the fact that an appropriate choice of the system's parameters
leads to the coexistence of two asymptotically stable equilibria in the
non-delayed case. When the total average time delay of the system is large
enough, the coexistence of two stable limit cycles is revealed, which
successfully model the ultradian rhythm of the HPA axis both in a normal
disease-free situation and in a diseased hypocortisolim state, respectively.
Numerical simulations reflect the importance of the theoretical results.
