## A Novel Use of Pseudospectra in Mathematical Biology: Understanding HPA Axis Sensitivity

Catherine Drysdale 1 and Matthew Colbrook 2

1 SMQB, University of Birmingham, Edgbaston, Birmingham, B15 2TT c.n.d.drysdale@bham.ac.uk ,

- 2 DAMTP, University of Cambridge, Wilberforce Rd, Cambridge, CB3 0WA

Abstract. The Hypothalamic-Pituitary-Adrenal (HPA) axis is a major neuroendocrine system, and its dysregulation is implicated in various diseases. This system also presents interesting mathematical challenges for modeling. We consider a nonlinear delay differential equation model and calculate pseudospectra of three different linearizations: a timedependent Jacobian, linearization around the limit cycle, and dynamic mode decomposition (DMD) analysis of Koopman operators (global linearization). The time-dependent Jacobian provided insight into experimental phenomena, explaining why rats respond differently to perturbations during corticosterone secretion's upward versus downward slopes. We developed new mathematical techniques for the other two linearizations to calculate pseudospectra on Banach spaces and apply DMD to delay differential equations, respectively. These methods helped establish local and global limit cycle stability and study transients. Additionally, we discuss using pseudospectra to substantiate the model in experimental contexts and establish bio-variability via data-driven methods. This work is the first to utilize pseudospectra to explore the HPA axis.

Keywords: HPA axis, pseudospectra, nonlinear delay differential equations, dynamic mode decomposition (DMD)

## 1 Introduction

The Hypothalamic-Pituitary-Adrenal (HPA) axis, a major neuroendocrine system, comprises the hypothalamus, pituitary gland, and adrenal glands. The hypothalamus releases corticotropin-releasing hormone (CRH), which stimulates the pituitary gland to secrete adrenocorticotropic hormone (ACTH). ACTH prompts the adrenal glands to produce cortisol, completing a negative feedback loop that inhibits further secretion of ACTH and CRH. These hormones follow circadian and ultradian rhythms. Three interacting behaviors of the HPA axis are noted: the circadian rhythm initiated by the hypothalamic suprachiasmatic nucleus (SCN) through CRH production, the ultradian pulsatility of ACTH and cortisol - arguably the axis's intrinsic dynamics as they persist without CRH, and activation by external stressors, which increase hormone amplitudes but not frequency. Fig. 1 illustrates this negative feedback loop. Dysregulation of

Fig. 1. A schematic of the HPA axis. The sun represents that CRH is entrained by the day and night cycle, which triggers an enveloping of cortisol and ACTH on a circadian level. On an ultradian level, ACTH triggers cortisol, which then inhibits ACTH. Both operate on regular pulses, the amplitude of which varies according to the other and the wider system. The lightning bolt represents external stressors that cause ACTH and, thereby, cortisol to rise. The dashed line represents the restriction we make to the negative feedback loop between ACTH and cortisol that we consider in this paper.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-28647f84-887e-4dda-810d-a8446e76c8fd.jpg)

the HPA axis, such as chronic over-activation resulting in sustained high cortisol levels or insufficient cortisol production, is implicated in various diseases, including anxiety disorders [1], depression [2], PTSD [3], Addison's disease [4], Cushing's syndrome [5], and chronic fatigue syndrome [6].

Mathematical modeling of the HPA axis has varied over the years, often targeting specific diseases or aligning with experimental observations. These models are typically nonlinear ordinary differential equations (ODEs) or delay differential equations (DDEs). Multiple mechanisms, such as time delays, negative feedback loops, nonlinearities, and varying timescales, can induce oscillations, as overviewed in [7]. The 'minimal model' [8] incorporates CRH, ACTH, and cortisol and assumes a differentiable negative feedback effect on ACTH and CRH. Under minimal assumptions, this model demonstrates oscillations through a Hopf bifurcation, albeit with unphysical parameters. Another model [9] features similar components but attributes oscillations to Hill-function type nonlinearities, noting the absence of oscillations with physiologically reasonable parameters. Interestingly, perturbations in this model led to two additional fixed points, indicative of 'hypercortisolemic' and 'hypocortisolemic' depressive groups, suggesting that mechanisms such as delays, rather than nonlinearity, might drive oscillations. This hypothesis is supported by [10], where a model including ACTH, cortisol, and glucocorticoid receptors demonstrated self-sustained oscillations with sufficient delay, despite CRH being modeled as a constant input.

While a normal form adeptly captures the bifurcation structure, it often fails to align well with empirical data. This issue may persist across the models above, regardless of the chosen model's correctness. The tension arises from the inter-

action between non-normality and nonlinearity, where existing tools focus either on one aspect or the other but seldom both. For example, studies in [11, 12] have shown that a non-normal linear operator, despite having stable eigenvalues, can give rise to a wave packet that interacts with the nonlinearity. Such interactions render first-order approximations, such as normal forms that depict bifurcation phenomena, not representative of the entire solution. Moreover, relying exclusively on linear tools, such as spectra and pseudospectra, can fail to capture the transient behaviors of the complete nonlinear system.

We seek a balanced approach by exploring different linearizations that help us comprehensively understand the complete nonlinear system. We explore transient behaviors for nonlinear systems by computing different types of pseudospectra for various linearizations, gaining different insights from each linearization.

We consider the following nonlinear DDE model of the HPA axis:

$$\frac { d A } { d t } = - e _ { a } A + \frac { h c ^ { m _ { 1 } } } { c ^ { m _ { 1 } } + [ C ( t - \tau _ { 1 } ) ] ^ { m _ { 1 } } }, \ \frac { d C } { d t } = - e _ { c } C + \frac { \beta [ A ( t - \tau _ { 2 } ) ] ^ { m _ { 2 } } } { a ^ { m _ { 2 } } + [ A ( t - \tau _ { 2 } ) ] ^ { m _ { 2 } } }, \ ( 1 )$$

where A denotes ACTH and C denotes cortisol. The initial conditions and constants are described in Section 2. This model, used in [13, 14], where analysis demonstrates the existence of a Hopf bifurcation, was chosen for its ability to exhibit oscillatory phenomena at physiologically reasonable parameter values. It enables focused analysis of the intrinsic dynamics of the HPA axis, particularly the ultradian rhythms of ACTH and cortisol. We seek to answer the question, 'What should we linearize around?' and explore three linearizations:

- · A time-dependent Jacobian with time-dependent spectra;
- · Linearization around the limit cycle, investigating whether there is inherent non-normality that could suggest transient effects;
- · A dynamic mode decomposition (DMD) approach, which transforms a finitedimensional nonlinear system (that approximates our original DDE system) into an infinite-dimensional linear system (with no delays). We apply the data-driven method ResDMD [15] to compute pseudospectra.

To our knowledge, this work is the first to use pseudospectra to explore the HPA axis. Different linearizations allow us to investigate different biological phenomena. For instance, by considering the time-dependent Jacobian, we see that there is instability corresponding to the upward slopes of cortisol, as in rat models when they are seen to become more aggressive if presented with a stressor on the upward slope of corticosterone 3 [16, 17]. On the downward slope, the non-normality index is low, suggesting the system is not sensitive to perturbation during this phase. This is the first mathematical modeling to match the experimental findings.

3 Corticosterone is the equivalent of cortisol in rats. We use mathematical models fitted to human data. Hence, cortisol best describes the variable C in Eq. (1). However, relevant experiments were done on rats, so we often focus on the similarities between the behavior of cortisol in our models and corticosterone in rats.

The paper is organized as follows. Section 2 provides a detailed description of the model, including the motivations for the parameters and the process of dedimensionalizing the model. Section 3, titled 'pseudospectra of the model three ways', presents the pseudospectra calculated from a time-dependent Jacobian, linearization around the limit cycle, and DMD. In the conclusion, we discuss our results and the contexts in which each can be useful for explaining various experiments and observations. Additionally, we explore how pseudospectra might be used for model validation, i.e., determining whether the perturbed model captures experimental results. We also consider the potential of using data-driven methods to capture the spectrum and pseudospectrum. Code producing all of the figures of the paper can be found here: https://github.com/miinadietrich/ hpa axis pseudospectra.

## 2 Model

We use the following constants from [13, 14] for the system in Eq. (1):

- -e a = 0 04, the rate of elimination of ACTH; .
- -e c = 0 01, the rate of elimination of cortisol; .
- -m 1 = 4, the number of cortisol modules that bind to the free receptors in the adrenal gland;
- -m 2 = 4, the number of ACTH modules that bind to the receptors in the pituitary gland;
- -a = 21, the Hill half-maximum constant for ACTH;
- -c = 6 11, the Hill half-maximum constant for cortisol; .
- -h = 7 66, the positive stimulation of CRH on ACTH secretion; .
- -β = 1, the rate of production of cortisol from the release of ACTH;
- -τ 1 = 15, the delay for secreted cortisol to have negative feedback on ACTH;
- -τ 2 = 15, the delay for cortisol to stimulate ACTH production in the adrenal glands.

We non-dimensionalize by letting A = ax and C = cy and t = τ/e c . We also scale the delays accordingly with t 1 = e τ c 1 and t 2 = e τ c 2 and set c 1 = e a /e c = 4, c 2 = 1 ( / ae c ) ≈ 4 76, . c 3 = β/ ce ( c ) ≈ 16 37 . . The result is our system of study:

$$\frac { d x } { d \tau } = - c _ { 1 } x + \frac { h c _ { 2 } } { 1 + [ y ( \tau - t _ { 1 } ) ] ^ { m _ { 1 } } }, \quad \frac { d y } { d \tau } = - y + c _ { 3 } \frac { [ x ( \tau - t _ { 2 } ) ] ^ { m _ { 2 } } } { 1 + [ x ( \tau - t _ { 2 } ) ] ^ { m _ { 2 } } } \quad ( 2 )$$

with initial conditions x t ( ) = x 0 for t ∈ -[ t 1 , 0] and y t ( ) = y 0 for t ∈ -[ t 2 , 0]. A default initial condition is ( x , y 0 0 ) = (0 8858 1 7461), arising . , . from solving Eq. (1) for a fixed point by linearizing and setting the left-hand-side equal to zero [13]. In the same vein as the bifurcation analysis, we also consider a range of physical h values, i.e., for different CRH inputs that would occur throughout the day as CRH operates on a circadian cycle. In [14], the default value of h was 7 66, but the full range was from 3 068 to 23. For the considered range of values . . of h , the full system has a limit cycle behavior that we plot in Fig. 2. This limit cycle varies continuously with h .

Fig. 2. Left: Trajectories with 100 randomly selected constant initial conditions (red dots) for h = 7 66 and limit cycle (green). Right: The limit cycle as . h is varied.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-3fa13a0c-0c86-45e7-9bc5-bbe4b789013d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-61395fb4-0139-458d-b618-71dd4f2d3a06.jpg)

## 3 Pseudospectra of the model three ways

We now consider three different families of pseudospectra associated with the above model. Readers from mathematical biology may be more familiar with spectra or eigenvalues, which are used to study a system's long-term stability (or instability). In other words, spectra capture the post-transient behavior of linear systems. In contrast, pseudospectra enable us to quantify sensitivity to perturbations and capture the transient dynamics of linear systems. An excellent account can be found in [18, Chapters 14-19].

For instance, consider a finite matrix A whose eigenvalues all have magnitude less than 1. The linear discrete-time dynamical system x k +1 = Ax k is stable since lim k →∞ ∥ A k ∥ → 0. However, the pseudospectral contours may protrude beyond the unit circle, indicating that there could be some transient growth before decay when the system is perturbed. In the case of pseudospectra calculated on a limit cycle, this could imply that the limit cycle takes longer to stabilize after a perturbation. We will discuss this point below and compute so-called Kreiss constants that capture this behavior.

Before looking at the different types of pseudospectra of linearizations of Eq. (2), it is instructive to consider the following three equivalent definitions of pseudospectra of finite matrices.

Definition 1 (Pseudospectra of Matrices). Let A ∈ C N × N and ϵ &gt; 0 . The ϵ -pseudospectrum Sp ( ϵ A ) is the set of z ∈ C defined equivalently by any of:

- · ∥ ( A -z ) -1 ∥ -1 ≤ ϵ
- · z ∈ Sp( A + E ) for some matrix E ∈ C N × N with ∥ E ∥ ≤ ϵ
- · ∥ ( A -z u ) ∥ ≤ ϵ for some vector u ∈ C N with ∥ u ∥ = 1 .

↦

We can see from the first point in the definition that eigenvalues are contained within the ϵ -contours of the function z →∥ ( A -z ) -1 ∥ -1 . It is often an excellent numerical check to ensure that spectral computations are accurate by verifying that the spectrum is contained within the pseudospectrum. The second condition describes how far the eigenvalues move in response to a size ϵ perturbation.

The third condition associates spatial structures (i.e., the vector u ) with the perturbed eigenvalues, which often have physical meaning.

Often, for nonlinear systems, there is a choice regarding how to obtain the linearization or the matrix A in the above definitions. This is more nuanced in the context of delay-differential equations. Furthermore, the ϵ -pseudospectra form distinct sets; therefore, comparing the pseudospectra of one operator to another is also nuanced. In each section, we provide a short preamble explaining what we hope to understand from this linearization and the method used with numerical details. Then, we present the results and their biological ramifications.

## 3.1 Time-dependent Jacobian: Pointwise linearization

We linearize Eq. (2) with a time-dependent Jacobian around an arbitrary time point rather than an equilibrium point. This approach is often used in feedback control [19, Chapter 12]. It can be viewed as representing a nonlinear system by a linear system over an instantaneous time window.

Method: Taking a base trajectory ( x , y 0 0 ), we write x τ ( ) = x 0 ( τ ) + ˜( x τ ) and y τ ( ) = y 0 ( τ ) + ˜( y τ ) . To first order, this yields the following linear DDE system:

$$\left ( \frac { \frac { d \tilde { x } } { d \tau } } { \frac { d \tilde { y } } { d \tau } } \right ) = A \left ( \tilde { x } \right ) + B \left ( \tilde { x } ( \tau - t _ { 1 } ) \right ) + C \left ( \tilde { x } ( \tau - t _ { 2 } ) \right ),$$

where A,B and C are the following matrices:

$$A = \begin{pmatrix} - c _ { 1 } & 0 \\ 0 & - 1 \end{pmatrix}, \ B = \begin{pmatrix} 0 - h c _ { 2 } \frac { m _ { 1 } [ y _ { 0 } ( \tau - t _ { 1 } ) ] ^ { m _ { 1 } - 1 } } { ( 1 + [ y _ { 0 } ( \tau - t _ { 1 } ) ] ^ { m _ { 1 } } ) ^ { 2 } } \\ 0 & 0 \end{pmatrix}, \ C = \begin{pmatrix} 0 & 0 \\ c _ { 3 } \frac { m _ { 2 } [ x _ { 0 } ( \tau - t _ { 2 } ) ] ^ { m _ { 2 } - 1 } } { ( 1 + [ x _ { 0 } ( \tau - t _ { 2 } ) ] ^ { m _ { 2 } } ) ^ { 2 } } & 0 \end{pmatrix}.$$

Taking Laplace transforms of (3) leads to ∆ τ ( λ ) = λI -A -Be -λt 1 -Ce -λt 2 . In particular, we obtain a time-dependent spectrum

$$\text{Sp} ( \Delta _ { \tau } ) = \{ z \in \mathbb { C } \, \colon \det ( \Delta _ { \tau } ( z ) ) = 0 \}$$

corresponding to a nonlinear eigenvalue problem [20]. The spectrum is generally infinite (in contrast to linear finite matrix eigenvalue problems), but for any c ∈ R , there are only finitely many eigenvalues λ with Re( λ ) ≥ c [21, Theorem 1.5]. In what follows, we use the MATLAB package TDS-CONTROL [22] to compute the (finitely many) eigenvalues in a right-half plane.

We are interested in pseudospectra, defined as

$$\text{Sp} _ { \epsilon } ( \Delta _ { \tau } ) = \{ z \in \mathbb { C } \, \colon \exists u \in \mathbb { C } ^ { 2 } \text{ s.t. } \| \Delta _ { \tau } ( z ) u \| \leq \epsilon, \| u \| = 1 \}.$$

Note how this generalizes the matrix case in Definition 1. One can show that the other two equivalent definitions in Definition 1 also hold. 4 To compute Sp ( ϵ ∆ τ ), we compute the smallest singular values of ∆ τ ( z ) over a grid of z points.

4 For the second point in Definition 1, we consider pencils λ → E λ ( ) bounded by ϵ .

We are interested in how the eigenvalues and pseudospectra vary with time τ (for the selection of base points ( x , y 0 0 )) and the parameter h to see if they match the experimental findings of more transient growth regarding upward peaks versus downward peaks of cortisol. We consider the spectral abscissa, defined as α = α τ ( ) = max λ ∈ Sp( ∆ τ ) Re( λ , ) and the distance to instability, defined as d τ = inf { ϵ &gt; 0 : ∃ s ∈ R , ∥ ∆ τ ( is ) -1 ∥ -1 = ϵ } whenever α τ ( ) &lt; 0. The system in Eq. (3) is (asymptotically) stable if and only if α τ ( ) &lt; 0 [21, Chapter 1]. In that case, d τ tells us the minimum size of perturbation needed to cause the system to become unstable. We also consider the ratio -α/d τ , whenever α &lt; 0. This index tells us the ratio between the distance of the spectrum to the imaginary axis and the minimum value of ∥ ∆ τ ( z ) -1 ∥ -1 on the imaginary axis. With an abuse of terminology, we shall refer to this as a non-normality index. 5

Results: Fig. 3 (panels (a) - (d)) plots α d , τ , and -α/d τ for h = 7 66. We also . show pseudospectra plots at different time points (panels (e) - (h)).

The instability regions correspond to the upward rise of cortisol and ACTH (see panel (a)). They also coincide with points where x (ACTH) is maximized, and y (cortisol) is minimized before it increases again. These observations strongly correlate with the findings in [16, 17], where rats became more aggressive under a stressor during the upward slope of corticosterone rather than during the downward peak. Additionally, the downward slope of cortisol is entirely within the stable region of the graph.

The index -α /d τ τ has peaks before and after the unstable region (panel (d)). These correspond to the dips in the distance to instability (panel (c)) and the spectral abscissa (panel (b)). Therefore, a small perturbation can easily bring us into the unstable half-plane, and this in itself should be considered as a 'nonnormal feature'. In [23], it was said that the HPA axis needs to be both 'sensitive to environmental perturbations, and able to respond differently to both small and large stimuli', as well as be 'robust with preservation of dynamic behavior during these perturbations'. We can associate the two different-sized peaks with different responses to small and large stimuli at different times. The fact that there is a plateau between these two peaks possibly corresponds to robustness and the system's desire to return to equilibrium following a perturbation but being sufficiently 'non-normal' at other points to respond to perturbation. These ideas have to be experimentally validated.

The leading eigenvalue splits into two eigenvalues (panel (e) to panel (f)) on the stable side of the plane. These eigenvalues coalesce again (panel (g)) before splitting and going into the right-half plane (panel (h)). The coalescence corresponds to the two unstable peaks we saw before, whereas the separation corresponds to the instability. Arguably, we have an 'intermittent' Hopf bifurcation. Such bifurcations have been used for the relation between CRH drive and delays [24], but not so far in the study of the ultradian rhythms.

5 Since ∆ τ ( λ ) is nonlinear in the spectral parameter λ , the phrase non-normality here is used in the sense of large pseudospectral regions enclosing eigenvalues.

Fig. 3. (a) - (d): Solution trajectories, spectral abscissa, distance to instability, and ratio -α/d τ . (e) - (h): Pseudospectra of the linearized system for different τ . We display pseudospectra by plotting contours of ϵ (colorbars). The eigenvalues are shown in red. The stability line (imaginary axis) is the dotted magenta line.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-2d4eaf06-a154-4736-976d-ce31b2b3fcbb.jpg)

This behavior is consistent across different values of h , i.e., independent of CRH input. Fig. 4 plots the maximum index -α/d τ and α as h varies.

## 3.2 Floquet pseudospectrum: Linearization about limit cycle

We now consider a periodic linear DDE system, corresponding to linearization about the limit cycle. Unlike linear periodic ODE systems, a complete Floquet theory does not generally exist. Nevertheless, we can still characterize stability

Fig. 4. Maximum α and index -α/dτ (over the limit cycle) as h varies.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-e728525c-3e58-4b11-8baa-1ff00bb4a727.jpg)

and transients through the spectral properties of an infinite-dimensional operator. 6

Method: If the base trajectory ( x , y 0 0 ) is periodic, then Eq. (3) is a linear periodic DDE system, where the matrices A B , and C are periodic in time. The stability of such systems is characterized by characteristic multipliers [25, Chapter 8], which are similar to those found in Floquet theory. 7 Let ω denote the period of the limit cycle. Consider the operator U , which maps a continuous function ϕ defined over the interval [ -t 1 , 0] to the solution of Eq. (3) over the interval [ -t 1 + ω, ω ], starting from the initial condition ϕ . The operator U acts on the Banach space C ([ -t 1 , 0]) of continuous functions on [ -t 1 , 0] equipped with the L ∞ norm. The stability and transient behavior of the system are determined by the spectra and pseudospectra of U . Pseudospectra of U are given by 8

$$\text{Sp} _ { \epsilon } ( U ) = \{ z \in \mathbb { C } \, \colon \exists \phi \in C ( [ - t _ { 1 }, 0 ] ) \text{ s.t.} \, \| ( U - z I ) \phi \| _ { L ^ { \infty } } \leq \epsilon, \| \phi \| _ { L ^ { \infty } } = 1 \}.$$

Computing pseudospectra in Banach spaces generally poses a significant challenge but is feasible in our case. Consider a discrete grid of equally spaced points s , s 1 2 , . . . , s N with s 1 = -t 1 &lt; s 2 &lt; · · · &lt; s N = 0. We utilize piecewise affine hat functions, ϕ j , such that ϕ j ( s i ) = δ ij . These functions form a basis of the piecewise affine functions on the interval [ -t 1 , 0] with knots at the points s j . Using this basis for the variables x and y , we compute U ϕ [ j 1 ; ϕ j 2 ] by solving the periodic DDE system and then evaluate the solution at the same grid points to form T ∈ C 2 N × 2 N which approximates U . For a given z , we compute the inverse

6 In our case, this operator is completely continuous. This implies that its spectrum is countable, zero is the only accumulation point of the spectrum, and all non-zero parts of the spectrum are eigenvalues of finite multiplicity.

7 It should be noted that a complete Floquet theory does not generally exist.

8 A subtlety arises here. The spectrum of U is independent of the initial time for which we study the system in Eq. (3), and the eigenspaces for different initial times are diffeomorphic. However, the pseudospectra can change. For consistency, as we vary h , we take a point at which ACTH peaks as the initial time.

## 10

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-31a222a1-1ac2-4038-9ce8-3cef4c7fa2be.jpg)

Fig. 5. Left: The pseudospectrum of T for the default parameter h = 7 66. The eigen-. values are shown in red, and the stability circle (unit circle) is the dotted magenta circle. Right: The Kreiss constant for different values of h (the dashed line is h = 7 66). .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-0121fafb-90e2-4827-b011-7ba1d4b19bc7.jpg)

-1 ∞

( T -zI ) and subsequently calculate its matrix l norm. This provides an approximation of ∥ ( U -zI ) -1 ∥ -1 , which characterizes the level sets of Sp ( ϵ U ). We also consider the Kreiss constant of U , K c ( U ) = sup ϵ&gt; 0 [sup {| z | : z ∈ Sp ( ϵ U ) } -c /ϵ ] = sup | z &gt;c | ( | z | -c ) ∥ ( U -z ) -1 ∥ defined for c &gt; sup {| z | : z ∈ Sp ( ϵ U ) } . These constants give a lower bound on transient behavior about the limit cycle through the relation [18, Chapter 16]: sup k ≥ 0 c -k ∥ U k ∥ ≥ K c ( U . )

Results: We take N = 50, which corresponds to a total of 100 basis functions. Fig. 5 (left) displays the pseudospectra for the default h = 7 66. There is a . dominant eigenvalue λ ≈ 0 9946 (which has converged in . N ) very close to one, while all other eigenvalues cluster at zero; however, large regions of pseudospectra are present. We consider c = 1 and found K 1 ( T ) = 7 4014. Fig. 5 (right) plots the . Kreiss constant for different values of h . After an initial peak, we see a decrease in transient effects with increasing h .

The limit cycle is stable, as expected in a hormone system. However, it may exhibit significant transient growth under perturbation, as exhibited by the large regions of pseudospectra. In particular, following a perturbation, the dynamics may take longer to stabilize compared to a system described by a normal matrix (whose pseudospectra are well-behaved). Given that the dominant eigenvalue is very close to one, these transient behaviors will likely be critically important. The variation in the Kreiss constant in Fig. 5 suggests that the limit cycle may be more susceptible to perturbations when CRH concentrations vary at different times of the day.

## 3.3 DMD analysis: Global linearization

In our final example, we consider a global linearization through Koopman operators. These operators act on a lifted infinite-dimensional space of functions of the state ( x, y ). We use data-driven algorithms to compute a Galerkin approximation of this operator and its pseudospectra from trajectory data.

Method: We first augment the state ( x, y ) to form the infinite state vector x ( τ ) = ( x τ ( ) , y ( τ ) , x ( τ -t 1 ) , y ( τ -t 1 ) , x ( τ -2 t 1 ) , y ( τ -2 t 1 ) , . . . ). Formally, this leads to an infinite system of nonlinear ODEs for x . As an approximation, we truncate x to d -1 time lags so that x ∈ R 2 d ( d = 10 in what follows), leading to an approximate system d x dτ = f ( x ) , for some nonlinear function f : R 2 d → R 2 d . Letting ω denote the time period of the limit cycle (which depends on h ), we sample this system with ∆τ = ω/ 10, yielding a discrete-time dynamical system

$$\mathbf x _ { n + 1 } = \mathbf F ( \mathbf x _ { n } ), \ \ n = 0, 1, 2, \dots.$$

Given a suitable function (called an observable) g : R 2 d → C , the Koopman operator K acts on g by composition with F :

$$[ K g ] ( { \mathbf x } ) = g ( { \mathbf F } ( { \mathbf x } ) ).$$

The key point is that K is linear. We have traded the finite-dimensional nonlinear system in Eq. (4) for a linear operator K that acts on an infinite-dimensional function space of observables. In general, the choice of function space matters, see [26, 27], but it is typical to take an L 2 space with respect to some measure. Since g ( x n ) = [ K g n ]( x 0 ), we see that K characterizes the forward-time dynamics of the system and its pseudospectra capture transient behavior. In particular, observables g associated with Sp ( ϵ K ) so that ∥ ( K -zI g ) ∥ ≤ ϵ ∥ g ∥ are physically significant and correspond to approximate coherence since ∥ ( K n -z n I ) g ∥ = ∥ g ( x n ) -g ( x 0 ) ∥ = O ( nϵ g ∥ ∥ ) .

Koopman operators have received considerable attention over the last decade; see [28, 27]. As well as dealing with nonlinearity, the considerable advantage of this approach is that we do not have to know the map F . One can also use experimental data. In what follows, we will build approximations of K using trajectory data. We will employ the Residual DMD (ResDMD) algorithm, which computes spectra and pseudospectra of K with rigorous convergence guarantees. A complete discussion of this algorithm is beyond the scope or space of this paper, but we refer the reader to [15, 29].

We collect data from 10 4 initial points selected uniformly at random in [ -3 5] , × -[ 1 8]. , We run the system in Eq. (2) forward by ( d + 1) × ∆τ to generate an initial x 0 as well as x 1 , and x 2 . In the language of DMD, the data consists of M = 2 × 10 4 snapshots formed by length-three trajectories for each initial condition. As our basis of observable, we use k -means on the data to select N = 400 centers in R 2 d , which we use with Gaussian radial basis functions.

Results: Fig. 6 (left) shows the pseudospectrum computed by ResDMD, as well as the DMD eigenvalues (red dots). On the circle, a lattice of eigenvalues corresponds to harmonics on the limit cycle generated by an eigenfunction g 0 with eigenvalue λ 0 (green dot). The lattice structure arises since the Koopman operator is multiplicative: K g g [ 1 2 ] = [ Kg 1 ][ Kg 2 ]. These harmonics are shown in Fig. 7, where we plot the eigenfunctions corresponding to λ 0 , λ 2 0 , and λ 3 0 across the entire state space R 20 . Here, the interpretation is that each two-dimensional

Fig. 6. Left: The pseudospectrum of K for the default parameter h = 7 66. The eigen-. values are shown in red, and the stability circle (unit circle) is the dotted magenta circle. Right: The Kreiss constant for different values of h (the dashed line is h = 7 66). .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-7864960e-b4d9-4244-9204-06e3571b9f0a.jpg)

block of x ∈ R 20 corresponds to marching the eigenfunction backward in time, capturing dynamics off the limit cycle.

Fig. 6 (left) also shows large regions of pseudospectra corresponding to transient (off-limit cycle) behavior. Fig. 6 (right) shows the Kreiss constant of the Koopman operator as h varies.

The physical interpretation of these results aligns with those in Section 3.2 in that it is physically reasonable for the system to be stable while exhibiting large regions of pseudospectra. However, the Kreiss constants vary significantly between the methods and show a different type of dependency on h . We should not be alarmed by the differences in the Kreiss constants, as the norms they are based on are fundamentally different and incomparable. Furthermore, the Koopman operator provides a global linearization of a dynamical system, in contrast to the local linearization discussed in Section 3.2.

## 4 Conclusion

Two fundamental questions in using pseudospectra for mathematical biology are: 'How should the matrix be obtained?' and 'How are clinically relevant perturbations modeled?'. We addressed the first question in the context of the HPAaxis by calculating pseudospectra through three methods: a time-dependent Jacobian, linearization around the limit cycle, and global linearization via Koopman operators. In addressing 'How should the matrix be obtained?', our initial model [13, 14] represents one of the many possible approaches. Future work could consider models incorporating CRH as a variable (as opposed to a parameter) or including glucocorticoid receptors. However, we want to highlight the role of spectra and pseudospectra in what we term 'model substantiation,' assessing how well a model aligns with experimental evidence. Spectra and pseudospectra are useful for identifying the limits of or substantiating models based on their perturbation response. Other tests for model substantiation include assessing

Fig. 7. Eigenfunctions of K corresponding to harmonics on the limit cycle (left). Moving from left to right, we march backward in time, capturing the transient dynamics off the limit cycle.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00845v2-figure-17e13bf3-f245-4625-91f2-5ebb12b1356a.jpg)

the stability of a limit cycle and the existence of oscillations under physically reasonable parameters.

Our first method, the time-dependent Jacobian, provided physically relevant results consistent with early rat experiments [16, 17], where aggressive behavior was observed when rats were perturbed on the upward slope of corticosterone secretion. Additionally, fluctuations in the non-normality index indicated that the rhythms of ACTH and cortisol could produce a range of responses within a single period. Linearization around the limit cycle showed that, although the cycle was stable, the response to perturbations could be disproportionate, as evidenced by the pseudospectra. The DMD method similarly demonstrated stability across all eigenvalues and limit cycle-like structures in the eigenvectors yet highlighted sensitivity to perturbations in the pseudospectra. For computing pseudospectra in the latter two cases, we developed novel mathematical techniques for Banach spaces and for DMD in delay differential equations. Interestingly, increasing h resulted in diminishing transient effects in local linearization while enhancing them in global linearization, as captured by the Kreiss constants.

Future work will focus on modeling clinically relevant perturbations. In [17], a male rat faced a stressor (a competitor rat) during both the rising and falling phases of cortisol. In [16], adrenalectomized rats received either a constant or a pulsatile dose of corticosterone; while the ACTH response was dampened with a constant dose, oscillatory corticosterone plasma levels and corresponding aggressive behavior were still observed. Our current models and pseudospectra methods do not account for extended perturbation durations or differences between injected and plasma corticosterone. To refine our understanding of aggres- sive behavior during unstable periods, as observed in Fig. 3, we propose using structured pseudospectra - specifying the E in the pseudospectra definition to better replicate experimental outcomes. Additionally, linearization around the limit cycle and DMD could elucidate the effects of hormonal misalignment, exacerbated by external stressors like sleep deprivation and meal timing [30]. Converting these stressors into structured ACTH perturbations will be a key focus.

Finally, data-driven approaches like DMD help bypass challenges in model substantiation but introduce difficulties in distinguishing perturbations from underlying rhythms or deriving perturbations from data. With established prior rhythms and modelable perturbations, we can investigate how structured pseudospectra on matrices align models with reality. This opens avenues for studying individual variability in perturbation responses and advancing personalized medicine. We advocate for strong collaboration between mathematicians and experimentalists to ensure data suits these sophisticated techniques.

## Bibliography

- [1] Carlo Faravelli et al. The role of life events and HPA axis in anxiety disorders: A review. Curr. Pharm. Des. , 18(35):5663-5674, 2012.
- [2] Carmine M. Pariante and Stafford L. Lightman. The HPA axis in major depression: Classical theories and new developments. Trends Neurosci. , 31 (9):464-468, 2008.
- [3] C.S. de Kloet et al. Assessment of HPA-axis function in posttraumatic stress disorder: pharmacological and non-pharmacological challenge tests, a review. J. Psychiatr. Res. , 40(6):550-567, 2006.
- [4] Andreas Barthel et al. An update on Addison's disease. Exp. Clin. Endocrinol. Diabetes. , 127(02/03):165-175, 2019.
- [5] Roberta Giordano et al. Ghrelin, hypothalamus-pituitary-adrenal (HPA) axis and Cushing's syndrome. Pituitary , 7(4):243-248, 2004.
- [6] Anthony J. Cleare. The HPA axis and the genesis of chronic fatigue syndrome. Trends Endocrinol. Metab. , 15(2):55-59, 2004.
- [7] B´ ela Nov´k and John J. Tyson. Design principles of biochemical oscillators. a Nat. Rev. Mol. Cell Biol. , 9(12):981-991, 2008.
- [8] Frank Vinther, Morten Andersen, and Johnny T. Ottesen. The minimal model of the hypothalamic-pituitary-adrenal axis. J. Math. Biol. , 63(4): 663-690, 2011.
- [9] Morten Andersen, Frank Vinther, and Johnny T. Ottesen. Mathematical modeling of the hypothalamic-pituitary-adrenal gland (HPA) axis, including hippocampal mechanisms. Math. Biosci. , 246(1):122-138, 2013.

- [10] Jamie J. Walker, John R. Terry, and Stafford L. Lightman. Origin of ultradian pulsatility in the hypothalamic-pituitary-adrenal axis. Proc. R. Soc. B , 277(1688):1627-1633, 2010.
- [11] Peter J. Schmid and Dan S. Henningson. Stability and Transition in Shear Flows . Number 142 in Applied mathematical sciences. Springer, New York, 2001.
- [12] Jean-Marc Chomaz. Global instabilities in spatially developing flows: Nonnormality and nonlinearity. Annu. Rev. Fluid Mech. , 37(1):357-392, 2005.
- [13] Hamed Malek, Mohammad Mehdi Ebadzadeh, Reza Safabakhsh, and Alireza Razavi. Mathematical analysis of the role of pituitary-adrenal interactions in ultradian rhythms of the HPA axis. Comput. Biol. Med. , 135: 104580, August 2021.
- [14] Hamed Malek et al. Dynamics of the HPA axis and inflammatory cytokines: Insights from mathematical modeling. Comput. Biol. Med. , 67:1-12, 2015.
- [15] Matthew J. Colbrook and Alex Townsend. Rigorous data-driven computation of spectral properties of Koopman operators for dynamical systems. Commun. Pure Appl. Math. , 77(1):221-283, 2023.
- [16] Ratna A. Sarabdjitsingh et al. Stress responsiveness varies over the ultradian glucocorticoid cycle in a brain-region-specific manner. Endocrinology , 151(11):5369-5379, 2010.
- [17] J. Haller et al. Ultradian corticosterone rhythm and the propensity to behave aggressively in male rats. J. Neuroendocrinol. , 12(10):937-940, 2000.
- [18] Lloyd N. Trefethen and Mark Embree. Spectra and pseudospectra: The behavior of nonnormal matrices and operators . Princeton University Press, 2005.
- [19] Hassan K. Khalil. Nonlinear systems . Prentice Hall, Upper Saddle River, N.J., 3 edition, 2002.
- [20] Stefan G¨ttel and Fran¸ coise u Tisseur. The nonlinear eigenvalue problem. Acta Numer. , 26:1-94, 2017.
- [21] Keqin Gu, Vladimir L. Kharitonov, and Jie Chen. Stability of time-delay systems . Birkh¨user Boston, 2003. a
- [22] Pieter Appeltans, Haik Silm, and Wim Michiels. TDS-CONTROL: A MATLAB package for the analysis and controller-design of time-delay systems. IFAC-PapersOnLine , 55(16):272-277, 2022.
- [23] Stafford L Lightman, Matthew T Birnie, and Becky L Conway-Campbell. Dynamics of ACTH and cortisol secretion and implications for disease. Endocr. Rev. , 41(3):bnaa002, 2020.
- [24] Eder Zavala et al. Mathematical modelling of endocrine systems. Trends Endocrinol. Metab. , 30(4):244-257, 2019.
- [25] Jack K. Hale. Theory of Functional Differential Equations . Springer New York, 1977.
- [26] Igor Mezi´ c. Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry. J. Nonlinear Sci. , 30(5):20912145, 2020.
- [27] Matthew J Colbrook. The multiverse of dynamic mode decomposition algorithms. Handbook of Numerical Analysis , 2024.

## 16 Drysdale and Colbrook

- [28] Steven L. Brunton, Marko Budiˇ si´, c Eurika Kaiser, and J. Nathan Kutz. Modern Koopman theory for dynamical systems. SIAM Rev. , 64(2):229340, 2022.
- [29] Matthew J Colbrook, Igor Mezi´, and Alexei Stepanenko. Limits and powers c of Koopman learning. arXiv preprint arXiv:2407.06312 , 2024.
- [30] Eder Zavala. Misaligned hormonal rhythmicity: Mechanisms of origin and their clinical significance. J. Neuroendocrinol. , 34(6):e13144, 2022.
Catherine Drysdale Matthew J. Colbrook
A Novel Use of Pseudospectra in Mathematical Biology: Understanding HPA Axis Sensitivity
The Hypothalamic-Pituitary-Adrenal (HPA) axis is a major neuroendocrine
system, and its dysregulation is implicated in various diseases. This system
also presents interesting mathematical challenges for modeling. We consider a
nonlinear delay differential equation model and calculate pseudospectra of
three different linearizations: a time-dependent Jacobian, linearization around
the limit cycle, and dynamic mode decomposition (DMD) analysis of Koopman
operators (global linearization). The time-dependent Jacobian provided insight
into experimental phenomena, explaining why rats respond differently to
perturbations during corticosterone secretion's upward versus downward slopes.
We developed new mathematical techniques for the other two linearizations to
calculate pseudospectra on Banach spaces and apply DMD to delay differential
equations, respectively. These methods helped establish local and global limit
cycle stability and study transients. Additionally, we discuss using
pseudospectra to substantiate the model in experimental contexts and establish
bio-variability via data-driven methods. This work is the first to utilize
pseudospectra to explore the HPA axis.