Coordination of autonomic and endocrine stress responses to the Trier Social Stress Test in adolescence.

Dysregulations in autonomic and endocrine stress responses are linked to the emergence of psychopathology in adolescence. However, most studies fail to consider the interplay between these systems giving rise to conflicting findings and a gap in understanding adolescent stress response regulation. A multisystem framework-investigation of parasympathetic (PNS), sympathetic (SNS), and hypothalamic pituitary adrenal (HPA) axis components and their coordination-is necessary to understand individual differences in stress response coordination which contribute to stress vulnerabilities. As the first investigation to comprehensively evaluate these three systems in adolescence, the current study employed the Trier Social Stress Test in 72 typically developing adolescents (mean age = 13) to address how PNS, SNS, and HPA stress responses are coordinated in adolescence. Hypotheses tested key predictions of the Adaptive Calibration Model (ACM) of stress response coordination. PNS and SNS responses were assessed via heart rate variability (HRV) and salivary alpha amylase (sAA) respectively. HPA responses were indexed by salivary cortisol. Analyses utilized piecewise growth curve modeling to investigate these aims. Supporting the ACM theory, there was significant hierarchical coordination between the systems such that those with low HRV had higher sAA and cortisol reactivity and those with high HRV had low-to-moderate sAA and cortisol responsivity. Our novel results reveal the necessity of studying multisystem dynamics in an integrative fashion to uncover the true mechanisms of stress response and regulation during development. Additionally, our findings support the existence of characteristic stress response profiles as predicted by the ACM model.

INTRODUCTION

Stress response and regulation are integral to social, cognitive, and emotional functioning (Boyce & Ellis, 2005). The stress response system (SRS) consists of three main physiologic systems: the parasympathetic (PNS) and sympathetic (SNS) branches of the autonomic nervous system (ANS) and the hypothalamic–pituitary–adrenal (HPA) axis which function to mobilize resources to allow individuals to adapt to perceived threats (Russell & Lightman, 2019). Adolescence is a sensitive time in the development and plasticity of SRS with heightened susceptibility to genetic and environmental interactions (Dahl & Gunnar, 2009; Roberts & Lopez‐Duran, 2019; Rotenberg & McGrath, 2016). Dysregulations in individual components of the SRS have been linked to deficits in cognitive functioning and emotional regulation, and these stress‐related vulnerabilities are crucial contributors in the etiology and severity of adolescent psychopathology (McEwen & Akil, 2020; McEwen & Morrison, 2013; McEwen & Stellar, 1993). Uncovering the mechanisms of individual differences in stress responsivity across these systems and their coordination is integral to prevent and treat stress‐related diseases.

New evidence suggests that the coordination between these systems may underlie individual differences in stress responsivity, which may precede the onset of psychopathology in adolescence (Buss et al., 2018). Thus, a multisystem framework—integrating comprehensive investigation of PNS, SNS, and HPA components and their coordination—is necessary to understand differences in SRS coordination which contribute to stress vulnerabilities. Furthermore, theory suggests greater coordination of SRS components emerges across development (Alkon et al., 2003, 2011) with solidification of an individual's SRS profile in adolescence (Del Giudice et al., 2011); adolescence is therefore a crucial window of investigation for SRS coordination.

Stress response systems

Much of our knowledge on stress physiology derives from studies investigating an isolated SRS component without accounting for multisystem influences. The ANS and HPA axis act on different temporal scales, with the ANS responding to stress within milliseconds to minutes and the HPA axis responding over minutes to hours following stressor onset (Engert et al., 2011). As a fast response system, activation of the ANS in anticipation of and across stress leads to a quick, rapid physiologic response through the actions of the PNS and SNS. Within milliseconds of stressor anticipation, the PNS withdraws its “rest and digest” effects on various organs, particularly on the cardiovascular system via dampening of vagal tone. This has the effect of elevating heart rate and therefore cardiac output to mobilize resources to peripheral organs. This change in vagal tone is often indexed by high‐frequency heart rate variability (HRV: variation in time interval between heart beats) which decreases with loss of vagal input (Kim et al., 2018). In both adults and adolescents, studies show a rapid decrease in HRV in response to psychosocial stressors (Berntson & Cacioppo, 2004; Taelman et al., 2009).

The SNS arm of the response—often referred to as the “flight‐or‐fight” mechanism—acts through catecholamines released from the adrenal glands and locus coeruleus (Lovallo, 1975) to enable individuals to rapidly adjust systemic regulation mechanisms to prepare organs to meet high energy demands (Rotenberg & McGrath, 2016). The rapid effects include increases in blood pressure and heart rate. One index of the catecholamine response is salivary alpha amylase (sAA) which reflects plasma noradrenaline levels. Salivary glands secrete sAA in response to noradrenergic activation by catecholamines (Chatterton Jr. et al., 1996; Thoma et al., 2012); thus, secretion is controlled by direct sympathetic activation. Additionally, sAA levels during and following stress have been associated with plasma catecholamines and other markers of sympathetic activation including cardiovascular indices and basal skin conductance levels (Bosch et al., 2003; Nater & Rohleder, 2009; Nater et al., 2005; Rohleder et al., 2004, 2006; van Stegeren et al., 2006). For these reasons, sAA is a reliable marker of SNS activity across a stressor. In adults and adolescents, sAA increases following acute laboratory stressors (Nater et al., 2006) marking it as a useful index of the SNS response.

Adolescents' HPA responses, the slow arm of the stress response, is mediated through the release of cortisol which has long‐lasting effects including glucogenesis, volume regulation via renal mechanisms, lipolysis, bone resorption, immune suppression, and a host of other multiorgan system effects (Michaud et al., 2008). This response is quantified noninvasively by salivary cortisol which increases to a peak (reactivity phase) around 20–25 min (Narvaez Linares et al., 2020) following stressor onset and decreasing to baseline levels (recovery phase) around 60 min post‐exposure (Allen et al., 2014; Ji et al., 2016; Kirschbaum et al., 1993; Seddon et al., 2020). Cortisol reactivity is thought to index an individual's sensitivity to a stressor (Linden et al., 1997). Cortisol recovery is postulated to index an individual's capacity to withstand external threats and is linked to individual styles of coping (McEwen, 2004; Meuwly et al., 2012). Thus, consideration of an individual's entire cortisol response trajectory is necessary to understand the complexities of the HPA stress response (Lopez‐Duran et al., 2014).

These well‐documented single system responsivities fail to account for the integration between the ANS and HPA (Bauer et al., 2002) which is required to understand healthy, adaptive functioning in SRS physiology (Del Giudice et al., 2011). This work can then be extended to determine if dysregulation in these biologic mechanisms contributes to the etiology and maintenance of psychopathology.

Coordination of stress response systems

While these facets of the SRS are anatomically distinct, they are theorized to coordinate and integrate their dynamics in a hierarchical fashion to calibrate physiologic reactivity and regulation to the perceived demands (Del Giudice et al., 2011). Del Giudice et al. developed the “Adaptive Calibration Model” (ACM) to describe the mechanisms of integration across the SRS and how individual differences in SRS may inform our understanding of stress vulnerabilities, particularly during adolescence when coordination between SRS is thought to solidify and signs of psychopathology begin to emerge. The ACM framework posits that systems “come online” in a hierarchical fashion depending on the environmental context or perceived threat.

The PNS is described as the “gate‐keeper” of this hierarchy—it has the largest influence at baseline prior to stress (“rest and digest”) and is thought to calibrate the degree of SNS and HPA recruitment during a stress response (Del Giudice et al., 2011). Via changes in vagal tone in response to environmental demands, the PNS acts as the “brake‐or‐acceleration” of further SRS activation (Porges et al., 1994). If individuals sufficiently adapt to their perceived demands by withdrawal of PNS vagal tone (indexed by lower HRV during stressors), no other systems will be recruited in the response. The ACM theory designates this as a “buffered” response profile in typically developing individuals (Del Giudice et al., 2011). If the body's demands are not met at a certain threshold of vagal withdrawal, recruitment of SNS and HPA systems will be triggered. The ACM hypothesizes that the shape of the HPA response—the slowest system to respond but with the longest lasting effects—depends on the degree of acute SNS activation and PNS withdrawal during stress. This hierarchical recruitment of all SRS systems is coined a “sensitive” response profile. Thus, the HPA axis and SNS are theorized to work additively in healthy individuals while the withdrawal of the PNS synergistically promotes SNS and HPA effects only after reaching an individual's threshold. These thresholds are theorized to drastically change in clinical populations, particularly individuals suffering from trauma and stress‐related disorders (Beauchaine, 2001; Del Giudice et al., 2011; Juruena et al., 2020). In order to fully comprehend stress response profiles and individual differences in SRS coordination in clinical samples, we must first understand these dynamics in healthy populations.

The ACM model may further help explain contradictory findings in the literature examining bivariate relationships between SRS. In line with a sensitive response phenotype, evidence reveals linear associations between stress‐induced increases in cortisol and sAA (Cacioppo et al., 1998; Engert et al., 2011; Foretic et al., 2020; Grillon et al., 2007) and decreases in HRV (Smeets, 2010; Weber et al., 2010) in child and adult populations. Pulopulos et al. found that decreases in stress‐related anticipation HRV were related to cortisol increases, but not recovery (Pulopulos et al., 2018). Yet, many other studies find no relationship between HRV and cortisol (Altemus et al., 2001; Heilman et al., 2008; Looser et al., 2010) nor cortisol and sAA (Altamura et al., 2018; Karhula et al., 2017; Nater et al., 2006; Valentin et al., 2015) across stress. In studies utilizing the Trier Social Stress Test (TSST: a popular laboratory psychosocial stressor), some findings showed cortisol output significantly related to decreased HRV, while others found no relationship at all (Giles et al., 2014; Laurent et al., 2016; Marques et al., 2010; Rotenberg & McGrath, 2016). The dynamics of SRS coordination laid out in the ACM theory may explain conflicting findings: individuals may show buffered or sensitive response patterns; thus, the true coordination cannot be understood as these studies failed to capture and examine all three systems.

SRS dynamics during adolescence

Adolescence is a critical transition point during development in which there is a confluence of social and biological changes impacting SRS dynamics (Boyce & Ellis, 2005). The ACM model postulates that adolescence serves as a “switch point” for SRS; using an evolutionary framework, these systems become more reactive to psychosocial triggers as organisms are more focused on mating and reproductive behaviors rather than promoting bodily growth (Del Giudice et al., 2011). During this complex switch point, SRS undergoes rapid change biologically and is sensitive to environmental perturbations (Korte et al., 2005). Del Guidice and colleagues propose that adolescence defines the emergence of buffered and sensitive characteristic stress response phenotypes which are most prevalent in typically developing populations. The sensitive profile is seen in adolescents recruiting all three systems hierarchically: high PNS withdrawal (low HRV) and high‐to‐moderate SNS (sAA) and HPA (cortisol) activation. The buffered profile is seen in those with high‐to‐moderate PNS activation and low SNS and HPA activation.

A body of evidence in single or two SRS components supports the emergence of SRS coordination across development (Alkon et al., 2003, 2011). In infants and children, studies show SRS is dominated by a single physiologic system, typically the ANS (Ji et al., 2016; Parent et al., 2019; Perry et al., 2012; Quas et al., 2014). As individuals age, cross system coordination is thought to strengthen and solidify. TSST studies have found augmented responses (lower HRV and higher cortisol) with increasing age, but the strength of the relationship between systems across age groups is not well documented (Giles et al., 2014). Total sAA and cortisol have been shown to increase across childhood and into early adolescence (Ellis et al., 2005; Rohleder et al., 2006) alongside more pronounced decreases in HRV (Alkon et al., 2003). However, another study in healthy 8–10 and 15–17 year‐olds showed no age differences in either arm of the ANS response (Salomon et al., 2000). This is the first study to explore SRS coordination across three main physiologic systems (PNS, SNS, and HPA axis) during adolescence and findings may deepen our understanding of SRS physiology during development.

The current study

Despite advances in our understanding of stress regulation and developments in stressor paradigms, coordination of the SRS in adolescents remains unclear; this poses a barrier to understanding SRS regulation and its contributions to psychopathology. By investigating the coordination of autonomic (both PNS and SNS) and HPA stress responses, the current study advances conceptualizations of adolescent stress responsivity and regulation via a multisystem approach in a population of healthy adolescents subjected to the TSST for Children. This study achieves these aims by addressing key questions in an adolescent sample: in line with the ACM theory, do we find evidence supporting the existence of distinct response profiles (buffered and sensitive) as indicated by (1) high PNS (HRV) related to low‐to‐moderate SNS (sAA) and HPA (cortisol) responsivity, and (2) low PNS related to high SNS and HPA responsivity?

METHOD

Participants

As a part of a larger study, 160 adolescents were recruited through various means (word‐of‐mouth, advertisements, flyers, emails to community schools). In order to examine individual differences in stress responses across a laboratory stressor and minimize confounds, we used data from a subset of individuals in the larger study who completed the TSST and had no DSM‐IV axis disorder, chronic medical condition, not on any medications known to affect the stress response (stimulants, contraceptives, psychotropic medications, corticosteroids, etc.), smoking, drug use, alcohol use, nor any recent illness (i.e., flu). A trained clinician determined the presence of a DSM‐IV axis disorder using an abbreviated version of the Structured Clinical Interview for DSM‐IV and confirmed this via electronic health records when appropriate. And 72 typically developing adolescents (mean age: 12.5 ± 2.3 SD, range: 9–16 years; 39 males) met these criteria and were included in this study. Written, informed consent was provided by all parents and assent for minors was completed by all participants. The study protocol was approved by UNC Chapel Hill and Duke Institutional Review Boards.

Study design and procedure

Following consent, clinical interviews, and neurocognitive testing for the larger ongoing study, participants completed two stressor protocols on separate visits—the Trier Social Stress Test (TSST) and Montreal Imaging Stress Task (MIST)—as part of a larger, ongoing study. The focus of this study is the examination of stress system coordination in response to the TSST in adolescents; thus, responses to the MIST were not the focus of this study. Separate tests revealed no impact of stressor order on stress responses (Section S8). The TSST stressor session (Figure 1) included an acclimation period where participants practice working memory tasks and fill out questionnaires, electrophysiology (EEG) set‐up and capping, resting state periods, and TSST stressor administration (Kirschbaum et al., 1993) and repeat resting state and working memory tasks. EEG imaging results will be presented in a future publication. Participants were given instructions not to eat or drink 30 min prior to arrival and refrain from any activities which may alter salivary neuroendocrine markers (limit caffeine intake, gum chewing, strenuous exercise or activity, etc.).

FIGURE 1

Overall experimental set‐up and design for the stressor protocol. This depicts the order of consent, clinical interviews, and neurocognitive assessment as part of the larger overall study. The stressor session consisted of EEG set up with questionnaires, tasks prior to the stressor (TSST), and a repeat of the same tasks following the stressor. The line underneath corresponds to the time saliva samples were taken for alpha amylase and cortisol. The “t” represents minutes relative to TSST onset. Alpha amylase (fast SNS response) was not assessed at t60, but cortisol (slow HPA response) was

Psychosocial stressor: Trier Social Stress Test for children (TSST)

The TSST is the benchmark laboratory psychosocial stressor (Allen et al., 2014; Kudielka et al., 2007; Liu et al., 2017; Narvaez Linares et al., 2020) with excellent efficacy in eliciting a multisystem stress response across autonomic, endocrine, behavioral, and psychological domains in healthy children, adolescents, and adult populations (Seddon et al., 2020). The TSST relies heavily on stressors associated with social evaluation, public performance, and situational uncontrollability (Kirschbaum et al., 1993). Across a wide range of TSST protocols, the timing of the stress response is constant, even though different TSST paradigms differ in the degree of autonomic and HPA axis stimulation (Khoury et al., 2020; Narvaez Linares et al., 2020; Von Dawans et al., 2011).

This study administered the TSST similarly to other published protocols (Pruessner et al., 2013). The participant performed this task in front of a panel of researchers who were instructed to refrain from any verbal or non‐verbal feedback to the adolescent. The TSST consisted of a 5‐min preparation period (instruction was to prepare a story to a prompt in front of a panel of researchers and a video camera), a 5‐min story presentation, and a 5‐min mental arithmetic task (serial subtraction task).

For the TSST, participants were brought into a room next to the EEG suite where researchers introduced participants to the “behavioral panelists” who would be judging them on various series of tasks. There were 2 separate researchers in the panel, one male and one female, who were instructed to remain neutral throughout the duration of the TSST. The participants were unaware of the arithmetic task ahead of time, thus adding in an element of uncontrollability to the paradigm, in line with the original method for children (Buske‐Kirschbaum et al., 1997). EEG caps and electrodes were still connected, and the participants were seated throughout the duration of the instructions, preparatory, story, and math portions of the experiment. The seated position was especially important to minimize excessive movements for neural recordings and control for intrinsic cardiodynamic responses produced when transitioning from a seated to orthostatic position. Participants were seated directly in front of the behavioral panelists, a video camera, and a microphone. They were told that they would be recorded as part of the TSST deception. In line with the original TSST for children (Buske‐Kirschbaum et al., 1997), the participants received a story prompt (see Section S5 for full prompt) and were told that they had to finish telling the story in an exciting way for the committee and that they would be judged against all the other participants.

Unlike the original experimental design for children, participants were not allowed to keep their notes after the preparation period. This added an element of uncontrollability as they entered the story period as the behavioral panelists instructed them to discard their notes and begin. They were then asked to begin telling their story for 5 min and instructed to speak clearly so that the microphone and video camera could accurately record them. If adolescents finished their story prior to the 5‐min mark or paused for an extended period, the panelists instructed them to continue for the duration in a neutral tone. In line with the original methods, participants then completed a mental arithmetic task for 5 min. Adolescents had to serially subtract 7 from 758 (9–11 years) or 13 from 1023 (12–16 years). With each mistake, the panelist interfered to say “Stop, please start again.”

Measures

Pubertal staging

Participants completed the Pubertal Development Scale (PDS) in order to assess their relative pubertal stage. The PDS is a self‐report questionnaire composed of questions that assess gonadal, adrenal, and neuroendocrine changes throughout puberty. This scale has been shown to be equally reliable to the Tanner scale, which requires physical examinations by trained clinicians (Carskadon & Acebo, 1993). The PDS correlates with gonadal and adrenal hormone concentrations as well as bone age during pubertal maturation (Hibberd et al., 2015; Schmitz et al., 2004; Shirtcliff et al., 2009). Point values are averaged to give a possible range of scores from 1–12. Continuous scores indicate progression through 5 pubertal stages: prepubertal, early pubertal, midpubertal, late pubertal, and postpubertal. Scores were used in final analyses to control for a participant's relative pubertal status.

Salivary cortisol

Adolescents provided saliva samples via passive drool for a total of 2 min. Collection timing (Figure 1) was designed to catch typical peak response 20–25 min post‐stress and full recovery around 60 min post‐exposure. A total of 5 cortisol samples were collected over the course of the EEG session. The first sample was collected 30 min prior to the TSST (t = −30). The second sample was obtained at the start of the TSST (t = 0). Saliva samples were then collected at 20, 35, and 60 min following onset of the TSST.

While many studies conduct their testing during a limited time window due to circadian effects, this often is not feasible for large studies with adolescents. Individuals were scheduled in the morning and afternoon, and time was controlled for statistically. The exact time of day was logged for each collected saliva sample and used in all cortisol analyses to correct for circadian patterns in cortisol release (Kudielka et al., 2004), mitigate individual variation in timing, and improve fidelity of analyses.

Saliva samples for cortisol and alpha amylase analyses were stored at −80°C in a secured lab at UNC until they were analyzed through the UNC Psychiatry Biobehavioral Core. Salivary cortisol levels were determined using a commercially available competitive enzyme immunoassay (EIA) kit and protocol available from Salimetrics, State College, PA as described previously (Raff et al., 2003). Saliva samples from the same participant were assayed in the same batch. The sensitivity of the cortisol assay is <.007 μg/dl with a standard range of .007 to 1.8 μg/dl. The intra‐ and inter‐ assay variation are 3.88% and 6.69%, respectively.

Salivary alpha amylase

Salivary alpha amylase (sAA) samples were collected and stored similarly to cortisol samples; however, given the fast response dynamics of the SNS compared to the HPA axis, we only obtained sAA from time points −30, 0, 20, and 35 min. We assayed for sAA using a commercially available kinetic reaction assay (Salimetrics, State College, PA). Re‐runs were conducted when the activities from the initial runs were too high (≥600 U/ml) or too low (≤2 U/ml) by adjusting dilution folds. On average, the intra‐ and inter‐assay coefficients of variance were less than 10%.

Heart rate variability

HRV consists of changes in the time intervals between consecutive heartbeats and serves as an index of vagal tone on the cardiovascular dynamics; HRV decreases with vagal withdrawal and can thus provide insights into PNS input to cardiac pacemaker cells. We collected cardiovascular data via two ECG electrodes which were attached to the left and right chest (Biosemi ActiveTwo System). Raw data were digitized at a rate of 1024 Hz. HR data were recorded during the rest session prior to the stressor and continuously during the stressor task. HRV data were processed following our previously published protocol (Corr et al., 2020). Individual's raw data were converted to interbeat interval (IBI) format and passed to Kubios HRV software for automated artifact correction (Tarvainen et al., 2014). The accuracy of this automated artifact correction was inspected by a single trained rater, and any files with abnormal or biologically implausible peaks that were inadequately corrected by the automated process were manually edited in CardioEdit Software (Brain Body Center, 2007; Porges & Bohrer, 1990). Files with excessive artifacts (greater than 5% of analyzed beats) were excluded from analysis. To ensure that the continuous HR data collected during the TSST were correctly segmented, start times for preparation, story, and math periods were recorded during task administration. Edited IBI files were split according to these recorded times. For the preparation, story, and math periods, to ensure that only data collected during each specific period were included, files were segmented beginning 45 s after the recorded start time and ending 45 s before the full 5‐min task completion. Resting state analysis duration was the same as TSST task analysis duration (3.5 min). To index the autonomic response to stress, we extracted the high frequency HRV in the 0.12–0.4 Hz range, in accordance with prior adolescent HRV research (Cui et al., 2015; McLaughlin et al., 2015).

Average HRV was extracted for the pre‐stress resting state, the TSST preparatory period (indexing anticipation of threat) (Pulopulos et al., 2018), and the TSST tasks (math and story periods representing within‐stressor adaptation). Given the high correlation between average HRV and age due to the nature of cardiovascular development across early childhood to adulthood (Finley & Nugent, 1995), these were then regressed against age and residuals served as the age‐corrected mean HRV indexing PNS function via vagal tone on the heart. Additionally, a relative change score from each individual's baseline was used to index relative changes from rest for each individual. These final values were incorporated in models and served as age‐corrected (removed model estimation difficulties arising from highly collinear covariates with a predictor of interest) and relative‐change during anticipation and stress HRV scores.

Salivary affect ratings

At each time of saliva collection, individuals were asked to complete a short 5‐question affect questionnaire to assess changes in mood over the course of the sessions. The questions asked participants how (1) stressed, worried, or nervous, (2) happy, relaxed, or comfortable, (3) irritated, annoyed, or mad, (4) sad, down, or unhappy, (5) overwhelmed, unable to control things, or discouraged they felt on a scale of 1–5 (not at all to very much).

Data analyses

Data preparation

Missing data examination showed that data met criteria for missing at random (see Section S4) and thus all individuals (n = 72) in the study were included in the final analyses as robust maximum likelihood allows for missing data—each individual contributes their available observations in the analysis (Bhat, 1974; Ivanova et al., 2016). We used the Box‐Cox power transformation for time series data to normalize cortisol concentrations (Miller & Plessow, 2013) which has been shown to outperform typical log transformations. As distributions of predictors are known and fixed in linear regression estimations, we did not need to normalize other measures prior to modeling (Lindstrom & Bates, 1990).

There is much debate concerning the most useful measurement index for sAA including: point estimates, percent changes from baseline, slope increase across stress, and area under the curve (AUC) indices. Given the complexity of measuring a fast response system via saliva, we chose AUC with respect to the increase (AUCi) as it is less sensitive to timing as an aggregate index and is associated with other autonomic and HPA measures (Ali & Pruessner, 2012; Balodis et al., 2010; Nater & Rohleder, 2009; Rohleder & Nater, 2009). AUCi represents the total increase in sAA above baseline levels in response to stress. In order to compute this measure, data need to be complete. The breakdown of missing sAA values for each timepoint (minutes −30, 0, 20, and 35) is as follows: 4, 4, 2, and 4 respectively. As Little's missing completely at random tests were non‐significant (p > .05) and less than 5% of data was missing at each time point, we performed multiple imputation for these missing data using the “multiple imputation by chained equations” (mice) package in R software (version 3.5.2) (Zhang, 2016). We calculated AUCi following the trapezoidal formula for sAA (Pruessner et al., 2003) to index SNS reactivity over the course of the stressor. AUCi values were used in all final analyses.

Cortisol response was modeled using two‐piece growth curve modeling with person‐centered knot points (also known as landmark registration) following the procedure for modeling of neuroendocrine data (Lopez‐Duran et al., 2014) (see Section S1 and Figure S1 for detailed description of this analytic approach). This piecewise procedure allows for examination of the entire cortisol trajectory while examining predictors of the specific phases of the HPA response (reactivity and recovery). In line with the original procedure, we used a person‐centered approach which allowed the knot point for each trajectory to fall on its natural post‐stress peak (landmark registration). Thus, the spline of the piecewise model captured the rise (reactivity) and fall (recovery) of each individual's cortisol trajectory. If the individual showed a plateau, the highest concentration within the plateau that was at least 10% greater than baseline (Lopez‐Duran et al., 2014) was used as the individual's peak or trajectory knot point. This approach aligns with the current literature determining a salient HPA response to stress (Ji et al., 2016). Thus, these individuals were classified as stress responders (identifiable post‐stress peak) and this time‐at‐peak was used to place the trajectory knot point. In line with Lopez‐Duran and colleagues, we used the mode peak time (20 min post‐TSST) as the trajectory knot point for non‐responders (those without identifiable peaks). This knot point does not alter the overall trajectory, only where we estimate the spline along the curve (Rahal et al., 2020). This estimation procedure allows for examination of the entire cortisol response curve and allows for phase specific predictors prior to and after the knot point (reactivity and recovery respectively) (Lopez‐Duran et al., 2014). To account for diurnal influences, we used sample time of day as the time growth predictor in the model. Time of day was converted to decimal time and mean centered.

Piecewise multilevel growth curve models were estimated within R software (version 3.5.2) lme4 package with restricted maximum likelihood estimation (Bates et al., 2015). These models, unlike traditional ANOVA approaches, can handle missing data. We used the RePsychLing package with principal components analysis to test for overfitting of random effects structures (Barr et al., 2013). p‐Values were estimated using Satterthwaite approximations to degrees of freedom and heteroskedasticity‐consistent (HC3) robust standard errors were reported. Continuous predictors were mean centered for interpretability. Standardized regression coefficients were reported to allow for comparison of effect sizes.

According to the original published piecewise growth curve procedure (Lopez‐Duran et al., 2014) for neuroendocrine data, model parameters include: the intercept set to each individual's peak time (peak activation or knot point of the curve), baseline cortisol, minutes to peak (reactivity slope), and minutes after peak (recovery slope) as fixed effects with random intercepts and slopes to account for individual variability in these parameters. The outcome variable in these models is the repeated cortisol concentrations. An individual's “peak activation” indexes cortisol level at the knot point of the response curve (Lopez‐Duran et al., 2014). Cortisol “reactivity” and “recovery” index the slope of cortisol as it approaches and moves away from the peak activation (Section S1).

Primary statistical analyses

Descriptive analyses were initially performed to describe mood ratings, cortisol, alpha amylase, and HRV. Repeated measures ANOVA using time as a within‐subject factor was used to assess the isolated physiological (cortisol, sAA, and HRV) and psychological (perceived “stress, worry, nervousness”) stress responses. Greenhouse–Geisser corrections were used when assumptions of sphericity were violated. Spearman correlations were performed to assess collinearity between model predictors. We also examined correlations between baseline SRS measures to determine if systems were related at rest. As the larger study design may have influenced cortisol responses to the TSST, we ran a linear model to examine whether stressor order (MIST vs. TSST) and days between stressor protocol completion significantly affected results (see Section S8). Ruling out potential confounds, we found no effect due to participation in multiple stressor protocols and thus, these covariates were not included in subsequent models.

In order to address hypotheses 1 and 2 (examine buffered and sensitive response types), we performed two separate piecewise growth curve models with main autonomic predictors (sAA and HRV) added hierarchically to examine their individual and then interactive effects on cortisol response curves: (1) test whether HRV during the preparatory period (representing anticipation of threat by the PNS), sAA, and their interactive effect were associated with cortisol response to the TSST, (2) test whether HRV during the stressor tasks (representing acute PNS adjustment to concurrent stress), sAA, and their interactive effect were associated with cortisol response to TSST. Models 1 and 2 allow for examination of SRS coordination following acute psychosocial stressor exposure while distinguishing between two conceptually different aspects of the PNS response (Pulopulos et al., 2018). Biological sex (Men = 0, Women = 1), pubertal status, age, and time of day were controlled for in final analyses as these covariates are known to influence physiologic stress responses (Liu et al., 2017).

RESULTS

Examining single systems across stress

Descriptive information for key study variables is presented in Table 1. Mean cortisol, sAA, HRV and affect ratings across the TSST are presented in Figure 2 (see Figure S2 for additional affect ratings). Correlational analyses did not reveal any bivariate relationships among baseline SRS indices including baseline cortisol, sAA, and HRV (p > .05) Table (S3).

TABLE 1

Descriptive statistics for key study variables

	Sample size	Mean (SD) or frequency (%)
Biologic sex	72
Males	39	54.2%
Females	33	45.8%
Race/ethnicity	72
Black or African American, non‐Hispanic	20	27.7%
White or Caucasian, non‐Hispanic	48	66.7%
Other	4	5.6%
Mother's education status	72
High School Graduate or below	4	6.2%
College and post‐graduate education	61	84.1%
No response	7	9.7%
Father's education status	72
High School Graduate or below	12	19%
College and post‐graduate education	51	68.5%
No response	9	12.5%
Age in years	72	12.46 (2.31)
Raw cortisol concentrations (μg/dl)
Sample 1 (−30 min pre‐TSST)	72	0.13 (0.06)
Sample 2 (0 min pre‐TSST)	72	0.11 (0.05)
Sample 3 (20 min post‐TSST)	72	0.15 (0.10)
Sample 4 (35 min post‐TSST)	72	0.12 (0.07)
Sample 5 (60 min post‐TSST)	72	0.09 (0.05)
Raw alpha amylase concentrations (U/ml)
Sample 1 (−30 min pre‐TSST)	72	163.4 (107.0)
Sample 2 (0 min pre‐TSST)	72	185.9 (136.0)
Sample 3 (20 min post‐TSST)	72	237.8 (160.8)
Sample 4 (35 min post‐TSST)	72	209.0 (166.4)
Alpha amylase AUC i	72	2205 (5306)
Heart rate variability (age‐detrended)
Pre‐TSST resting state	65	−0.049 (0.13)
Preparation (TSST)	67	−0.59 (0.16)
Story (TSST)	67	−0.18 (0.15)
Math (TSST)	57	−0.14 (0.13)
Post‐TSST resting state	60	−0.055 (0.15)
Pubertal status (PDS score)	72	7.83 (3.04)

Note: Adolescent sex coded male = 0, female = 1.

Abbreviation: TSST, Trier Social Stress Test.

FIGURE 2

Mean and standard errors for affect ratings, salivary alpha amylase (sAA), heart rate variability (HRV), and salivary cortisol across the stressor session. Stressor (TSST) onset at time = 0. (a) Indicates on average how”stressed, worried, or nervous” individuals felt. (b) Average sAA (measuring SNS response) across the experiment. Note that saliva was not analyzed at the 60‐min timepoint. (c) Average HRV change from baseline (indexing PNS response) across the TSST periods. The x‐axis indicates the recording period starting with pre‐stress rest. (d) Average cortisol response across stress

Mode cortisol peak time was 20 min post‐TSST onset and over 80% of the participant population were responders (those with identifiable cortisol peaks post‐stress); this is considered within an acceptable range for a stressor task (above 70%) (Kirschbaum et al., 1993; Miller et al., 2016). During the TSST session, individuals showed typical cortisol responses (Table 2). Specifically, individuals' cortisol levels significantly increased (β = 0.32, SE = 0.03, p < .001) prior to their peak (reactivity phase) and significantly declined (β = −0.49, SE = 0.04, p < .001) following their peak when controlling for all covariates. In line with expected autonomic stressor responses, HRV significantly decreased from baseline (Figure 2) during stress [F(4.1, 148.3) = 14.37, p < .001] while sAA significantly increased from baseline across the TSST [F(3.0, 284) = 15.84, p = .016]. Aligning with descriptive cortisol findings, mode alpha amylase peak time was also 20 min post‐TSST onset. After quantifying changes in isolated systems, we then examined coordination across systems.

TABLE 2

Estimates for growth curve models with landmark registration of the cortisol response to psychosocial stress predicted by autonomic stress response

	Baseline piecewise growth curve models without autonomic interaction terms
	Model 2.1: Alpha amylase (SNS)			Model 2.2: Preparatory HRV (PNS anticipation)			Model 2.3: Task (story + math) HRV (PNS acute stressor adaptation)
	β	SE	t‐Value	β	SE	t‐Value	β	SE	t‐Value
Intercept	0.032	0.14	0.23	−0.050	0.17	0.75	−0.061	0.17	−0.36
Time before peak	0.32	0.038	8.24***	0.34	0.042	8.08***	0.33	0.041	8.15***
Time after peak	−0.46	0.043	−10.63***	−0.44	0.046	−9.48***	−0.45	0.046	−9.81***
sAA	0.16	0.099	1.64	–	–	–	–	–	–
Time before peak × sAA	0.040	0.036	1.09	–	–	–	–	–	–
Time after peak × sAA	−0.15	0.054	−2.83**	–	–	–	–	–	–
HRV	–	–	–	−0.20	0.11	−1.75	−0.23	0.11	−2.14*
Time before peak × HRV	–	–	–	−0.15	0.044	−3.36***	−0.21	0.045	−4.59***
Time after peak × HRV	–	–	–	0.087	0.045	1.95	0.10	0.044	2.34**

Note: Models 2.1 through 2.3 examine coordination between isolated autonomic systems and the cortisol response trajectory. 2.1 examines the effects of the sympathetic (SNS) system or sAA (using AUCi as a predictor). 2.2 and 2.3 examine the effects of different aspects of the parasympathetic (PNS) system: HRV during the preparatory period and HRV during the TSST story and math tasks. Intercept represents peak activation (when all other predictors are 0 and time is 0). Time before peak indexes cortisol reactivity slope (x = time, y = cortisol concentration), time after peak reflects cortisol recovery slope. Models control for baseline cortisol, biologic sex, time of day, and age in years.

p < .05

p < .01

p < .001.

SRS coordination across stress

We started by examining relationships between two response systems and hierarchically built on these baseline models by including their interactive effects to eventually assess coordinated dynamics between PNS, SNS, and HPA systems. A correlation analysis examining the relation between ANS system components (PNS with SNS) revealed no significant associations between HRV during preparation period nor HRV averaged during the stressor tasks (story and math) with sAA (p > .05) Table (S3). Initial piecewise growth curve models to identify the interaction of PNS and SNS with the HPA response revealed that a greater increase in sAA across the TSST was associated with a flatter cortisol recovery slope (β = −0.15, SE = 0.054, p = .005) (Figure 3, panel A), but sAA had no effect on the cortisol peak nor the reactivity slope (p > .05, see Table 2, model 2.1). Conversely, an increased HRV during the preparatory period was associated with a flatter cortisol reactivity slope (β = −0.147, SE = 0.044, p = .001) (Figure 3, panel B), but not cortisol peak nor recovery slope (p > .05, see Table 2, model 2.2). Increased HRV during the stressor tasks (see Table 2, model 2.3) was associated with a flatter cortisol reactivity (β = −0.21, SE = 0.045, p < .001), flatter recovery slope (β = 0.103, SE = 0.044, p = .020) (Figure 3, panel C), and a reduced cortisol peak (β = −0.232, SE = 0.11, p = .039).

FIGURE 3

Plots (a–c) show the bivariate relationship between autonomic indices and cortisol response trajectories. X axes are cortisol concentration plotted against time either before (time before peak) or after (time after peak) the knot point in the piecewise growth curve model. Thus, the lines represent cortisol reactivity or cortisol recovery (x = time, y = cortisol, slope = reactivity or recovery) respectively. Shading around lines represent 90% confidence interval of the slope of the line. (a) Line colors depict mean sAA ± 1 standard deviation (SD) used only for visualization purposes (continuous variable in model). Individuals with higher sAA had steeper cortisol recovery slopes. (b) Line colors depict mean HRV during the preparatory period ±1 SD used only for visualization purposes. Individuals with lower preparatory HRV had a sharper increase in their cortisol reactivity slope. (c) Line colors depict mean HRV during the TSST tasks (story and math) ± 1 SD used only for visualization purposes. Participants with lower HRV had a steeper recovery slope

To test the combined, coordinated effects of the fast PNS and SNS responses on the slower HPA response, we hierarchically added their interactive effects to the above initial piecewise model (Table 3). The inclusion of these interactions in the piecewise model allowed for investigation of sensitive and buffered responses: high PNS (HRV) related to low‐to‐moderate SNS (sAA) and HPA (cortisol) responsivity, and low PNS related to high SNS and HPA responsivity respectively. Our findings revealed that HRV during stress (story and math tasks) significantly moderated the relationship between sAA and cortisol reactivity slope (β = −0.115, SE = 0.048, p = .019) (Table 3, model 3.2). Specifically, individuals with low HRV differed in their cortisol reactivity based on level of sAA—those with high sAA showed steeper cortisol reactivities (Figure 4). However, HRV during stress did not moderate the relationship between sAA and cortisol peak nor recovery slope (p > .05). Unlike the HRV during stress, there was no significant interaction between the preparatory HRV and sAA on cortisol response trajectories (p > .05) (Table 3, model 3.1).

TABLE 3

Estimates for growth curve models with landmark registration of the cortisol response to psychosocial stress predicted by autonomic stress responses

	Hierarchical piecewise growth curve models with autonomic interaction terms
	Model 3.1: Alpha amylase (SNS) × preparatory HRV (PNS anticipation)			Model 3.2: Alpha amylase (SNS) × task HRV (PNS acute stressor adaptation)
	β	SE	t‐Value	β	SE	t‐Value
Intercept	0.002	0.17	0.001	0.001	0.17	0.009
Time before peak	0.34	0.043	7.86***	0.33	0.041	7.99***
Time after peak	−0.47	0.048	−9.87***	−0.48	0.047	−10.14***
sAA	0.20	0.12	1.67	0.21	0.12	1.79
Time before peak × sAA	0.068	0.041	1.67	0.089	0.041	2.15*
Time after peak × sAA	−0.16	0.059	−2.61**	−0.19	0.060	−3.08**
HRV	−0.17	0.11	−1.55	−0.20	0.11	−1.88
Time before peak × HRV	−0.15	0.044	−3.35***	−0.19	0.046	−4.17***
Time after peak × HRV	0.16	0.046	2.52*	0.13	0.045	2.88**
sAA × HRV	−0.11	0.12	−0.93	−0.13	0.13	−0.98
Time before peak × sAA × HRV	−0.083	0.044	−1.87	−0.12	0.048	−2.38*
Time after peak × sAA × HRV	0.12	0.061	1.96	0.11	0.059	1.77

Note: Results from proposed models examining ANS dynamics on the slow HPA response trajectories are combined above. Interaction terms were added hierarchically to baseline models presented in Table 2. Model 4.1 examines the interactive effects of the SNS and the anticipatory stressor response of the PNS. 4.2 examines the interactive effects of the SNS and the acute acclimation of the PNS during the stressor. Intercept represents peak activation (when all other predictors are 0 and time is 0). Time before peak indexes cortisol reactivity slope (x = time, y = cortisol concentration), time after peak reflects cortisol recovery slope. Models control for baseline cortisol, biologic sex, time of day, age in years, and pubertal scores.

p < .05

p < .01

p < .001.

FIGURE 4

HRV moderates the relationship between sAA and cortisol reactivity slope. Plotting effects from the piecewise growth curve model examining cortisol trajectory before the knot point. Panel depicts individuals with HRV during the TSST tasks (story and math) 1 standard deviation (SD) below the mean. Y axis is the cortisol concentration. X axis is the time before peak (representing reactivity time or growth prior to the knot point). Line colors depict mean sAA ± 1 SD used only for visualization purposes (continuous variable in GCM). Shading around lines represent 90% confidence interval of the slope of the line. In those with lower HRV, cortisol reactivity (x = time, y = cortisol, slope = reactivity) varied with levels of sAA: Higher sAA related to steeper cortisol increases

DISCUSSION

The present study is the first to report coordination of autonomic (PNS and SNS) and HPA stress responses across the TSST in an adolescent sample. In accordance with the ACM theory of SRS dynamics, findings supported the existence of buffered and sensitive stressor response profiles (Del Giudice et al., 2011). This work highlights the necessity of examining these three main physiologic systems to understand the mechanisms of stress response coordination.

Buffered and sensitive stress responses

Analyses revealed significant coordination among SRS systems: higher SNS activation (increasing sAA) was associated with greater HPA activation (steeper increases in cortisol reactivity), but this co‐activation was only evident in those individuals with greater PNS withdrawal (decreasing HRV in response to stress). In those with less PNS withdrawal (or greater vagal tone), there was no significant increase in SNS and HPA activation. This is in line with our hypotheses of hierarchical SRS coordination suggesting the presence of both sensitive (high PNS withdrawal with high SNS and HPA activities) and buffered (low‐to‐moderate PNS withdrawal with low‐to‐moderate SNS and HPA activities) stressor phenotypes in typically developing adolescent populations. Our findings provide support for the presence of characteristic SRS dynamics in adolescence and synchronization between the systems which may be governed at the outset by PNS withdrawal or vagal tone (Porges et al., 1994); but at a threshold, the combined influences of SRS activation determine the course of an individual's stress response.

Buffered responses are characterized by withdrawal of PNS activation (vagal withdrawal indexed via decreased HRV in response to stress) without subsequent activation of SNS or HPA systems. The immediate PNS response is thought to represent an individual's threshold for further stress response activation—those who successfully accommodate to demands of the threat via PNS withdrawal do not trigger robust SNS and HPA responses (Porges, 2009). This buffered profile is reflected in our findings showing that less PNS withdrawal (higher vagal tone) across stress significantly associated with less activation in SNS and HPA systems. However, once a threshold of PNS withdrawal is reached for an individual without successful calibration of physiologic systems to meet the acute stressor demands, SNS and HPA systems become activated.

Sensitive responses are characterized by hierarchical recruitment of all SRS systems—PNS activation (vagal withdrawal indexed via decreased HRV in response to stress) and SNS and HPA activation (increased sAA and cortisol, respectively, in response to stress) (Del Giudice et al., 2011). Anticipatory PNS withdrawal to stress (indexed by decreasing HRV during TSST preparation) was significantly associated with increased HPA activation (indexed by steeper cortisol reactivity). Greater PNS withdrawal and regulation during stress (indexed by decreasing HRV within the math and story stressor tasks) significantly related to greater HPA activation (increased cortisol reactivity slope) and recovery (steeper cortisol recovery slopes). Furthermore, greater SNS activation in response to stress (indexed by increased sAA output) significantly related to prolonged HPA activation (flatter cortisol recovery slope).

Sensitive and buffered SRS profiles not only reflect the susceptibility of an individual to psychosocial threats, but the level of top‐down control over whole‐system activation (Del Giudice et al., 2011). Adolescents with increased PNS activity (vagal tone), a hallmark of a buffered response, are less reactive physiologically to psychosocial stressors, and this autonomic stability maintained by the PNS across stress may in turn attenuate HPA axis reactivity (conceptualized as an individual's sensitivity to a stressor) (Wolff et al., 2012). This could represent a potential resiliency factor in the overactivation of these systems and mitigate the long‐term effects of chronic, repeated elevations in cortisol; more longitudinal work is needed to determine if this response pattern is protective against the negative effects of chronic stress.

Ultimately, stress response profiles must be considered in a context dependent manner. A buffered phenotype may only be advantageous in certain scenarios. If an individual encounters an extreme threat that requires a heightened and sustained physiologic stress response, a sensitive phenotype would be the appropriate response. In this context, failure to activate a sensitive response would be maladaptive (lack the mobilization of appropriate resources in response to the degree of the environmental threat).

Multi‐system stress coordination

Hierarchical addition of autonomic interactions in the piecewise modeling procedure allowed for identification of the unique effects of the PNS and SNS on the HPA response above and beyond their combined dynamics. Additionally, analyses revealed the nature of multisystem stress coordination by accounting for shared (PNS × SNS interaction) and unique effects on HPA response trajectories. Some isolated effects of single ANS systems no longer remained significant when accounting for multisystem coordination (Table 3). For example, the effect of an individual's PNS withdrawal during stress on HPA peak activation did not remain significant when accounting for SNS activation. When the interaction between both arms of the ANS are incorporated, the sole predictor of HPA responsivity is the interaction between the PNS and SNS; individual effects of the PNS and SNS were no longer significant. This reveals a critical gap in existing literature—failure to examine all SRS elements may misattribute true regulation mechanisms to single components and fail to reveal dynamics that are dependent on coordination between systems.

Pulopulos and colleagues' examination of preparatory HRV—reflecting anticipatory PNS engagement in response to impending stress—revealed a significant effect of this PNS measure on HPA response to stress (Pulopulos et al., 2018, 2020). While our findings suggest preparatory HRV had a small isolated significant effect on cortisol reactivity, we failed to find a significant moderation effect of preparatory HRV on HPA response to stress when controlling for SNS activation. These discrepant findings may be due to the inclusion of SNS effects which were absent in the Pulopulos study. Furthermore, the relative effect size of preparatory HRV on cortisol reactivity compared to that of within stressor task HRV and SNS activation (sAA) was small (Table 3); this suggests PNS anticipation of stress may play a small contributory role in HPA activation, but is less influential than SNS activation or the within‐stressor PNS response. Our findings exemplify the importance of examining all three systems—PNS, SNS, and HPA—when investigating SRS coordination.

Discrepancies in findings may be accounted for by several methodological factors. A majority of studies which find relations between HRV and cortisol reactivities or peaks do not account for the covariance between an individual's reactivity, peak, and recovery—the entire profile or shape of the HPA response. While cortisol reactivities and recoveries are conceptually different, they are dependent on one another and can be confounded when extracting isolated slopes, area under the curve, or point estimates (Lopez‐Duran et al., 2014). For instance, a significant effect of HRV on a cortisol reactivity estimate may be driven by an individual's peak threshold or a slow, prolonged recovery if the maxima used in the computed slope is truly a part of the recovery period. Area under the curve estimates may equivocate vastly different HPA response curves and fail to capture distinct reactivity and recovery dynamics driving significant effects. Failure to incorporate an individual's entire HPA response profile and opting for single estimates represent misspecifications and over‐attribution of significance to the extracted parameters included in the model.

From a theoretical perspective, our findings may reveal why numerous inconsistencies in SRS literature exist (Altemus et al., 2001; Balodis et al., 2010; Cacioppo et al., 1998; Engert et al., 2011; Karhula et al., 2017; Khoury et al., 2020; Kudielka et al., 2004; Laurent et al., 2016; Looser et al., 2010; Maruyama et al., 2012; Myers et al., 2017; Nater et al., 2006; Quas et al., 2014; Rahal et al., 2020; Smeets, 2010; Valentin et al., 2015; Weber et al., 2010). Not only is it vital to consider the entire HPA response trajectory (to account for significant covariation in cortisol reactivity, peak, and recovery), but a multisystem approach using indices across SRS is needed as significant relationships between two systems disappear when accounting for the combined interactions across all three systems. Without taking into account the interconnected nature of the PNS, SNS, and HPA systems in determining the course of an acute stress response, studies may be misattributing the underlying mechanisms to the wrong system. Failure to account for these coordinated dynamics creates an incomplete framework for understanding stress response regulation and compounds inconsistencies across studies.

Limitations and future directions

There are a few limitations of the current study which are important to consider alongside our findings. Wake time was not collected in the current study, which prevented us from taking into account its potential impact on cortisol as done in prior research (Morgan et al., 2017). While we controlled for exact time of day in all analyses to remove diurnal confounds, further studies should conduct testing between a limited afternoon window to remove this potential source of variability in the design. Second, while the sampling in the study was designed to capture the distinct phases of HPA response (reactivity, peak, recovery), it is possible we missed the exact time at peak which may have reflected a point truly in the individual's reactivity or recovery phase; we mitigate and account for this potential confound by examining the entire trajectory of an individual's response profile, thus simultaneously controlling for each respective phase.

Increased sAA sampling frequency during the stressor should be considered to assess SNS responsivity on the appropriate time scale. Other stressor studies measuring sAA have implemented a longer sampling protocol with saliva collection following stressor completion and demonstrate an increase in sAA values at 15–20 min following stressor onset (Anesiadou et al., 2021; Kang, 2010; Ma et al., 2018; Nater et al., 2005; Thoma et al., 2012). For example, Gordis et al demonstrated sAA peak 15–20 min post‐stressor onset (Gordis et al., 2006). We attempted to mitigate our sampling frequency by using an area under the curve approach to estimate total sAA output (Nater et al., 2006; Pruessner et al., 2003).

While age and pubertal status failed to show a significant effect on stress responses, these two covariates are highly collinear and confounded, making their unique contributions difficult to isolate in cross‐sectional studies. Longitudinal investigations should examine how age and pubertal status affect SRS dynamics over time to determine their distinct influences. Additionally, researchers should also control for phase of menstrual cycle in older adolescents as this may impact between and within‐sex differences in SRS coordination, particularly the HPA response to acute stress (Dahl & Gunnar, 2009). Future studies should test the role of parental influences, rearing environment, and past experiences on SRS dynamics in line with the ACM theory to investigate whether sensitive or buffered types are determined by early childhood environments in typically developing adolescents (Del Giudice et al., 2011).

Conclusions

In line with the ACM model of stress response dynamics, we found evidence of hierarchical SRS coordination and two stress response profiles in typically developing adolescents: buffered and sensitive. Those with high PNS tone (evidenced by higher HRV across stress) had attenuated SNS (lower sAA) and HPA (flatter cortisol reactivity and recovery) responses as postulated for those with buffered phenotypes. With more PNS withdrawal (evidenced by decreasing HRV), SNS and HPA response systems collinearly activated as theorized for sensitive response profiles. This is the first time these profiles and patterns of SRS dynamics have been examined in adolescence. Our findings highlight the need for future studies to consider the interactive effects across all three systems of response: PNS, SNS, and HPA axis.

CONFLICT OF INTEREST

Authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent‐licensing arrangements), or non‐financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.

AUTHOR CONTRIBUTIONS

Sarah Glier: Conceptualization; data curation; formal analysis; investigation; methodology; project administration; software; validation; visualization; writing – original draft; writing – review and editing. Alana Campbell: Conceptualization; investigation; project administration; software; supervision; writing – review and editing. Rachel Corr: Data curation; investigation; writing – review and editing. Andrea Pelletier‐Baldelli: Investigation; writing – review and editing. Mae Yefimov: Investigation; project administration. Carina Guerra: Data curation; investigation; writing – review and editing. Kathryn Scott: Data curation; investigation. Louis Murphy: Investigation. Joshua Bizzell: Investigation. Aysenil Belger: Funding acquisition; investigation; project administration; resources.

TABLE S1 Estimates for the growth curve model with landmark registration of the cortisol response to psychosocial stress and the effects of known cortisol confounders

TABLE S2 Correlations between baseline measures

TABLE S3 Correlations between growth curve model predictors

TABLE S4 Testing for differences in cortisol responses based on session order. Individuals' reactivity slopes, recovery slopes, and peak activations were extracted from the GCM‐LR and used in Welch's t‐tests to determine if there were underlying group differences based on session order. Estimates from Welch's t‐tests are presented below

TABLE S5 Testing for differences in cortisol responses based on timing differences (in days) between stressor sessions. Individuals' reactivity slopes, recovery slopes, and peak activations were extracted from the GCM‐LR and used in separate linear models to determine if there were underlying effects of length of time between stressor sessions. Estimates from separate regressions are below. Time in days between sessions was the predictor in all models and individual parameter estimates from the GCM‐LRs were the dependent variables in separate regression models

FIGURE S1 Simulated cortisol trajectories for 2 individuals. (a) Demonstrates timing confounds when examining reactivity and recovery slopes in traditional models. (b) Landmark registration (LR) procedure centers individual curves on time 0 (person‐centered knot points) so individuals' reactivity, peak, and recovery slopes can be investigated in a single piecewise growth curve model

FIGURE S2 Mean and standard errors for affect ratings across TSST. Stressor onset at time = 0. Each graph depicts the scores for a different affect question seen in the title above each individual graph

Click here for additional data file.