Hormonal dysregulation after prolonged HPA axis activation can be explained by 2 changes of adrenal and corticotroph masses

11 Stress activates a complex network of hormones known as the Hypothalamic-Pituitary- 12 Adrenal (HPA) axis. The HPA axis regulates multiple metabolic, immune, and behavioral end- 13 points. Dysregulation of the HPA axis is a hallmark of several psychiatric conditions, including 14 depression and substance-abuse disorders. However, little is known about the origin of this 15 dysregulation, and we currently lack mathematical models that can explain its dynamics on the 16 timescale of weeks. Specifically, accumulated experimental evidence indicates that after 17 prolonged activation of the HPA axis, ACTH responses become blunted, and this blunting 18 persists for weeks even after normalization of cortisol levels and responses. Here, we use 19 mathematical modelling to show that this dysregulation can be explained by changes in the 20 functional masses of the glands (total mass of the cells) that secrete ACTH and cortisol. These 21 mass changes occur because the hormones CRH and ACTH regulate the growth of corticotroph 22 and adrenal cortex cells, respectively. Furthermore, we show that impaired glucocorticoid 23 receptor (GR) feedback exacerbates this dysregulation, providing a rationale for the role of GR 24 in depression. We propose that this dysregulation is a side-effect of a circuit with a 25 physiological function, in which gland-mass changes provide dynamical compensation to 26 physiological variation. These findings suggest that gland-mass dynamics may play a role in the pathophysiology of stress-related disorders.


Introduction 30
Physical and psychological challenges activate a network of endocrine interactions, 31 known as the hypothalamic-pituitary-adrenal (HPA) axis ( Figure 1A, (1-4)). The HPA axis 32 generates an adaptive stress response, and its dysregulation is implicated in a range of diseases 33 with elevated cortisol levels due to CRH secretion by the placenta. 3-6 weeks after delivery, cortisol levels and dynamics return to normal, 70 whereas ACTH dynamics are blunted-data from (22). After 12 weeks, ACTH dynamics normalize as well. In all panels, patient data are 71 denoted by thin gray line, and controls by a thicker black line.

73
The persistent blunting of ACTH responses after the resolution of hyper-cortisolemia 74 is surprising. The dynamics of ACTH in the classical HPA axis model is strongly associated 75 with the dynamics of cortisol, and so once cortisol normalizes, so should ACTH. It is also 76 unclear how a deficient ACTH response produces a normal cortisol response, given that ACTH 77 is the main regulator of cortisol secretion. Explaining this dysregulation thus requires a process 78 on the scale of weeks that decouples the dynamics of ACTH and cortisol. This timescale cannot 79 be readily explained by existing models of the HPA axis, where the relevant timescale is the 80 lifetime of hormones, which is minutes to hours (30). One important process with potentially 81 a timescale of weeks is epigenetic regulation of GR sensitivity (31-33). This process, however, 82 effects of mass changes, we assume that the per-unit-biomass secretion rates " and ' are 150

constant. 151
In order to describe the dynamics of the functional masses, we model the known major 152 factors that affect cell mass, namely the upstream hormones. CRH increases the growth of the 153 functional mass of pituitary corticotrophs, and ACTH increases the growth of the functional 154 mass of the adrenal cortex. After perturbations, the doubling time of the corticotrophs is on the 155 order of days in rats (47), while that of the adrenal mass is on the order of weeks in rats and 156 other animal models (34). 157 These hormone-driven mass changes are described by two additional equations, for the 158 masses of the corticotrophs and the adrenal cortex (Figure 2A and Methods). The important 159 parameters in these equations are the turnover times of the functional masses, which we 160 estimate to be on the order of days-weeks in order to explain the observed hormone dynamics 161 in Figure 1. Parameter values are given in Table 1. Note that the dynamics of the functional 162 masses is thus much slower than the dynamics of the hormones which have half-lives of 163 minutes to hours. 164

Model shows HPA axis dysregulation after prolonged activation 165
We simulated the response of the HPA model to a prolonged stress input. We model 166 the effect of a prolonged stressor as a pulse of input u to the HPA axis which lasts for several 167 weeks or longer. 168 169 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint of HPA axis dysregulation observed in experiments. At the end of a prolonged stress period, the adrenal mass is enlarged and the corticotroph 175 mass is slightly enlarged, which results in hypercortisolemia and blunted ACTH responses in the CRH test (B). After a few weeks, corticotroph 176 mass drops below baseline, while adrenal mass is slightly enlarged, causing normal cortisol dynamics with blunted ACTH responses to the 177 CRH test (C). Finally, after a few months both tissue masses return to normal, leading to normalization of both cortisol and ACTH dynamics 178 (D). In all panels, simulations are of a CRH test (Methods), where 'case' is after stress and 'control' is after basal HPA axis activation, as 179 described in Figure 3.

181
We then simulated a CRH test at several time-points after the end of the input stress 182 pulse. The CRH test is modelled by adding CRH and measuring ACTH and cortisol over 4 183 hours of simulated time (methods). We find three phases of recovery. In the first few weeks 184 after the end of the input pulse, cortisol is high and ACTH responses are blunted ( Figure 2B). 185 Then, for a period of several months, cortisol has returned to its original baseline but ACTH is 186 remains blunted ( Figure 2C). Finally, after several months, both cortisol and ACTH return to 187 their original baselines ( Figure 2D). Thus, the model recapitulates the observed dynamics of 188 The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint has several distinct phases of dysregulation during and after prolonged stress (Figure 3). These 193 phases are caused by changes in the functional masses of the pituitary corticotrophs and adrenal 194 cortex ( Figure 3AB). These functional mass changes lead to changes in the level of the 195 hormones ( Figure 3C). We also simulated a CRH test at each timepoint and calculated the 196 response as the peak hormone level after CRH administration relative to a control simulation 197 without the prolonged stress input ( Figure 3D). We plotted the CRH response as a function of 198 time for ACTH and cortisol in Figure 3E. A blunted response corresponds a response below 199 hormones over several hours. The response is defined as the maximum response to a CRH test (case) relative to the maximum response to a 206 CRH test in steady-state with basal input conditions (control). In E, the response to a CRH test given at time t is shown as a function of t. thus, 207 this is the predicted response to CRH tests done at different days during and after the stressor. The model (black lines) shows blunted (reduced) 208 responses to CRH tests after the stress similar to those observed in Figure 1, as well as a mismatch between cortisol and ACTH dynamics that 209 develops a few weeks after cessation of stress.

210
The initial phase occurs after the onset of the stressor and before adaptation to the 211 stressor (Marked ONSET in Figure 3). It lasts several weeks. In this phase, the increase in input 212 u causes elevated levels and responses of hypothalamic CRH, ACTH and cortisol. However, 213 over weeks of stress input, the corticotroph and adrenal masses grow. The gland masses thus 214 effectively adjust to the stressor, as example of the more general phenomenon of dynamical 215 compensation in physiological systems (61), discussed below. This functional mass growth 216 causes a return to baseline of hypothalamic CRH and ACTH, due to negative feedback by 217 cortisol. More precisely, a larger adrenal functional mass means that less ACTH is needed to 218 produce a concentration of cortisol that drives ACTH down to baseline. 219 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint Such a return to baseline is called exact adaptation. Exact adaptation is a robust feature 220 of this circuit due to a mathematical principle in the functional mass equations called integral 221 feedback (61) (Methods). Exact adaptation does not occur in models without the effects of 222 functional mass changes-the hormones do not adapt to the stressor ( Figure S1, Supplementary 223 Section 1). 224 The enlarged functional masses results in elevated cortisol levels during the stress 225 period, but in adapted (that is, baseline) levels of CRH, ACTH, and blunted responses of CRH 226 and ACTH to inputs ( Figure 3DE). During prolonged stress there is thus a transition from an 227 elevated to a blunted response of ACTH that occurs due to changes in functional masses, and 228 results from cortisol negative feedback. 229 At the end of the prolonged stress pulse, which we term Early Withdrawal (or EW), the 230 functional masses are abnormal, and take weeks to months to recover. This fundamental 231 process is the reason for the hormonal dysregulation that is the subject of this study. In the first 232 weeks after the stressor is removed, the adrenal and corticotroph functional masses shrink, 233 accompanied by dropping cortisol levels. ACTH responses are blunted, and blunting may even 234 worsen over time. 235 Then, cortisol and CRH levels and responses simultaneously normalize. This marks the 236 beginning of the next phase which we term Intermediate Withdrawal (IW, Figure 3). In this 237 phase, ACTH responses remain blunted, despite the fact that cortisol is back to baseline, 238 because adrenal functional mass is enlarged, and corticotroph functional mass is deficient.  We conclude that the model is sufficient to explain the dynamics of recovery from 248 chronic HPA activation in several conditions mentioned in the introduction -anorexia, alcohol 249 addiction, and pregnancy ( Figure 3). In order to explain the timescales of recovery, the only 250 model parameters that matter are the tissue turnover times. Good agreement is found with 251 turnover times on the scale of 1-3 weeks for the corticotrophs and adrenal cortex cells (see 252 Figure S2 for comparison of different turnover times, Supplementary Section 2). The model 253 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint therefore explains how ACTH responses remain blunted despite normalization of cortisol 254 baseline and dynamics. 255 A model in which functional masses remain constant does not show these phenomena, 256 but instead shows an immediate normalization of hormone levels and CRH-test responses after 257 stress ( Figure S1, red lines). We also tested several alternative mechanisms with a slow 258 timescale of weeks. We tested models with constant gland masses and the following processes 259 to which we assigned time constants on the order of one month: GR resistance following HPA 260 activation ( Figure S1, purple lines), slow changes in input signal ( Figure S1, blue lines), and 261 slow changes in cortisol removal rate ( Figure S1, gray lines). None of these models shows the 262 dysregulation that we consider. The reason is that these slow processes do not cause a mismatch 263 between ACTH and cortisol needed to capture the blunting of ACTH despite normal cortisol 264

responses. 265
Deficient GR feedback exacerbates HPA dysregulation following prolonged stress 266 We next considered the role of glucocorticoid receptor (GR) in recovery from 267 prolonged HPA activation. GR mediates the negative feedback of cortisol on CRH and ACTH 268 secretion. Impaired feedback by GR is observed in many cases of depression. The relation between depression and impaired GR function seems paradoxical, since 275 GR signaling mediates many of the detrimental effects associated with high cortisol levels such 276 as hippocampal atrophy (68). One explanation for the association between impaired GR 277 feedback and depression is that impaired GR feedback leads to failure of the HPA axis to 278 terminate the stress response on the timescale of hours, leading to excessive cortisol levels 279 (7,14). GR feedback also plays a role in ultradian and circadian rhythms in the HPA axis 280 (12,57,69). 281 Here we present an additional explanation for the association between GR feedback  The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint reduced when GR affinity is high. In other words, dysregulation is more severe the weaker the 287 feedback from GR. 288 289 Figure 4. GR provides resilience to the HPA axis against prolonged stressors. Here we show the dynamics of the HPA axis during and 290 after a prolonged pulse of stress, described in Figure 3, for +, =2 (strong feedback), +, =4 (moderate feedback), and +, =8 (weak feedback).

291
Stronger feedback from the GR attenuates the dysregulation of all HPA axis hormones.

293
The reason for this effect is as follows. After an increase in stress levels ( . → " ), the 294 adrenal gland mass increases ( Figure 4B). When GR feedback is weak ( +, ≫ " ), the adrenal 295 increases by about a factor of " / . . However, if GR feedback is strong ( +, ≪ " ), the 296 adrenal increases to a smaller extent, because less cortisol is required to inhibit ACTH to the 297 level required for precise adaptation. The smaller adrenal mass means a smaller dysregulation 298 of cortisol and ACTH after stress. Strong GR feedback therefore provides resilience to the HPA 299 axis against stress on the slow timescale. 300

Mass changes provide robustness and dynamic compensation to the HPA axis 301
Finally, we discuss several potential advantages provided by the functional mass 302 changes in the HPA axis. The first advantage is a natural way for the cell populations to 303 maintain their steady-state mass, by balancing growth and removal. Growth and removal must 304 be precisely balanced in order to avoid exponential growth or decay of the tissue. Thus, 305 although cells turn over, the feedback through the hormones couples with the hormone trophic 306 effects to set the masses at a functional level. 307 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint The second beneficial feature is the robustness of the steady-state hormone levels with 308 respect to physiological parameters. The simple form of the equations for the masses C(t) and 309 A(t) (given in Figure 2)  . This is remarkable because all other model parameters do not affect these 312 steady-states. Similarly, the cortisol baseline level, 45 ≈ = 9 > 9 8 7 > 7 8 ? ./@ > , depends on very 313 few parameters. As appropriate for a stress-response hormone, cortisol steady state depends 314 monotonously on the input to the hypothalamus (averaged over weeks), which corresponds 315 to physiological and psychological stresses. It is independent on almost all other parameters, 316 including production rate per cell ( " , ' ) or removal rates ( " , ' ) of CRH and ACTH, as 317 well as the rates of the proliferation and removal of the adrenal cortex cells D and D . This 318 robustness is due to the ability of the functional masses to grow or shrink to buffer changes in 319 these parameters. 320 This robustness also makes the steady-state level of cortisol and ACTH independent of 321 total blood volume. This is because blood volume only enters through the production 322 parameters " , ' . These parameters describe the secretion rate per cell, and since hormone 323 concentration is distributed throughout the circulation, these parameters go as one over the 324 blood volume (b . relates to the hypophyseal portal vein and not to total blood volume (70)). 325 The functional masses therefore adjust to stay proportional to blood volume. 326 The third feature of the present model is dynamical compensation, in which masses C 327 and A change to make the fast-timescale response of the system invariant to changes in the 328 production rates " , ' . A full analysis of dynamical compensation in the HPA axis is provided 329 in the SI (Supplementary Section 3). 330

Discussion 331
In this study we sought to understand the physiological mechanisms that underlies HPA 332 axis dysregulation after prolonged stress. In particular, we focused on a puzzling phenomenon 333 which recurs in several conditions: ACTH responses remain blunted for a few weeks after 334 cortisol dynamics return to baseline. This cannot be readily explained by existing models of 335 HPA axis activation. We show that a model which incorporates two known interactions that 336 govern functional masses -control of corticotroph growth by CRH and control of adrenal 337 growth by ACTH -is sufficient to explain this dysregulation. The dysregulation is a robust 338 property of this model, and explains clinical data with physiologically reasonable parameters. 339 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint This suggests that feedback interactions between masses and hormones may underlie the 340 dysregulation of the HPA axis after prolonged stress. 341 There may be several implications of the blunted ACTH response that lasts for weeks 342 after stress is removed. One set of psychopathological complications is due to the fact that 343 ACTH secretion is tightly linked with the secretion of endogenous opioids such as β-endorphin.  The present mass model can further provide insight into a physiological hallmark of 355 several psychiatric disorders: deficient GR feedback. Impaired GR feedback is implicated in 356 depression (7), and several genetic and environmental factors have been associated with this 357 impaired feedback (83,84). Deficient GR feedback has been suggested to affect ultradian (57) 358 and circadian (69) HPA rhythms. We find that strong GR feedback (high GR affinity to 359 cortisol) can also work via functional mass dynamics on the timescale of weeks: it protects the 360 HPA axis from dysregulation after prolonged activation. Reduced GR feedback causes larger 361 dysregulation after prolonged stress. This effect is directly mediated by changes in adrenal and 362 corticotroph masses: strong feedback allows smaller changes of adrenal mass during the stress 363 period, and hence to a smaller dysregulation after stress is removed. 364 Since chronic stress can by itself lead to reduced GR sensitivity ("glucocorticoid 365 resistance") (31), the present findings suggest that prolonged stress makes the HPA less 366 resilient to the next prolonged stress. Perhaps this reduced resilience is relevant to the 367 progression of depressive episodes (85). 368 If the present model is correct, functional masses are potential targets to address HPA 369 dysregulation. Measuring these masses and their dynamics using imaging may test this model 370 and provide clinically relevant information. Interventions that seek to normalize functional 371 masses can potentially reduce the extent of dysregulation during and after prolonged stresses. 372 One class of interventions may use control-engineering approaches, by periodically measuring 373 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2020.01.01.892596 doi: bioRxiv preprint masses or hormones and administering defined HPA agonists or antagonist doses in order to 374 reach desired functional masses and hormone levels (86). 375 The interactions between functional masses and hormones which underlie the present 376 HPA axis dysregulation also provide exact adaptation to CRH and ACTH (Figure 3) during 377 stress. Exact adaption has been extensively studied in the context of biochemical circuits (87), 378 but less is known about the role of exact adaptation in physiological circuits, and, in particular, 379 about its relevance to psychopathologies. One exception is mass changes in the insulin-glucose 380 system, in which exact adaptation can provide control of glucose levels and dynamics (61) 381 despite changes in insulin resistance, a form of dynamical compensation. It would be 382 fascinating to study what functional roles exact adaptation has in the HPA axis, in order to 383 understand the trade-offs that impact its dysregulation. The present modelling approach may 384 be able to address such questions. 385 386 387 author/funder. All rights reserved. No reuse allowed without permission.

HPA axis model. 389
We model HPA axis using the following equations, where cortisol provides negative 390 feedback on CRH secretion through MR, and negative feedback on both CRH and ACTH 391 secretion through GR. The functional mass of corticotrophs in the pituitary is C, and that of the 392 cortisol-secreting cells in the adrenal cortex is A. We assume a linear dependence of secretion 393 on the upstream hormones (55 We added to this classic description of the HPA axis the following equations for functional 400 mass dynamics, C(t) and A(t). The equations describe functional mass growth due to the 401 upstream hormones, and mass decline at turnover rates W and D , which are much smaller than 402 the hormone turnover rates: 403 In the simulations of Figure 2,  Where the gamma parameters are the turnover rates. All parameters are provided in Table 1 For high cortisol levels to induce strong resistance, we use the simple forms: f(x)=1, g(x)= 1 + 434 " . 435 The second alternative slow process is a slow decrease in input, u. We considered an 436 exponentially decreasing input signal ( Figure S1). 437 author/funder. All rights reserved. No reuse allowed without permission. The dose D and pulse width W were calibrated to provide the observed mean CRH test results 457 in non-stressed control subjects, providing D=20 and W=30 min. 458

Proof for dysregulation of ACTH after cortisol normalization. 459
We briefly show how the model Eq. 8-14 provide ACTH dysregulation even after cortisol 460 normalizes. Consider a prolonged stressor, like the one presented in Figure 3 (that is, a pulse 461 increase in the input u). Since ACTH and CRH are adapted to the stressor, within hours after 462 withdrawal of the stressor, CRH and ACTH levels drop to below their pre-stressor baseline, 463 whereas cortisol remains above its pre-stressor baseline. CRH and cortisol recover 464 simultaneously over a few weeks because CRH dynamics depend only on cortisol and the input 465 u. It therefore remains to be shown that ACTH does not recover before cortisol. Since ACTH 466 recovers from below baseline ACTH=1, it recovers with a positive temporal derivative 467 author/funder. All rights reserved. No reuse allowed without permission. > 0. If (by negation) this happens before CRH and cortisol recover, we will come to a 468 contradiction: since CRH is below its baseline CRH=1, C has a negative derivative due to Eq. 469 13; the derivative of A is zero due to Eq 7. Because ACTH levels are approximately 470 proportional to > >… † > b † > >… † > (Eq. 15) we find that ACTH has a negative temporal derivative 471 when it crosses its baseline. We come to a contradiction. We conclude that CRH and cortisol 472 return to baseline before ACTH does. 473 Cortisol shows normal fast timescale response after recovery despite ACTH blunting. 474 We now use the HPA model Eq. 8-14 to show that when CRH and cortisol levels return to 475 baseline after stress (point IW in Figure 3 proportional to the product of gland masses AC (Eq. 15-16), as noted above, one obtains D W = 481 ] is also independent of 491 D , W . We conclude that after cortisol and CRH return to baseline, the fast timescale dynamics 492 of cortisol and CRH to any input (such as a CRH test) becomes normalized and equal to the 493 dynamics before the stressor. In the simulations, after CRH and cortisol have normalized for 494 the first time, they only deviate from baseline slightly (at most 9%), and so this consideration 495 holds approximately for further times in the scenario of Figure 3. 496 497 author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was not peer-reviewed) is the .