Dynamic changes in the central autonomic network of patients with anorexia nervosa

Autonomic cardiac dysfunction is a common complication in patients with anorexia nervosa (AN). Despite its high prevalence, physicians often overlook this clinical condition, and little research has been dedicated so far. To probe the functional role of the neurocircuitry underpinning the poorly understood autonomic cardiac dysfunction, we examined dynamic functional differences in the central autonomic network (CAN) between 21 acute AN individuals and 24 age, sex and heart rate‐matched healthy controls (HC). We assessed functional connectivity (FC) changes in CAN using seeds in the ventromedial prefrontal cortex, left and right anterior insular cortex, left and right amygdala and dorsal anterior cingulate cortex. The overall FC between the six investigated seeds is reduced in AN individuals compared to HC, although no changes were observed for single connections. Moreover, AN exhibited higher complexity in the FC time series of such CAN regions. Contrary to HC, we found that the degree of complexity between FC and heart rate (HR) series did not correlate in AN, suggesting a shift from central to peripheral control of the heart in AN individuals. Using dynamic FC analysis, we showed that the CAN transits across five functional states with no preference for any. Strikingly, at the state of weakest connectivity, the entropy significantly diverges between healthy and AN individuals, reaching its minimum and maximum values, respectively. Overall, our findings provide evidence that core regions of the CAN engaged in cardiac regulation are functionally affected in acute AN.


| INTRODUCTION
Anorexia nervosa (AN) is a severe mental disorder mainly affecting girls and young women, showing the highest mortality rate among all mental disorders (Arcelus et al., 2011;Papadopoulos et al., 2009). The most common characteristics of AN are an excessive restriction of food intake and an intense fear of gaining weight. The constant state of energy deprivation and malnutrition may eventually lead to severe cardiovascular complications (Sachs et al., 2016). From a physiological perspective, the adaptation to prolonged low energy state, presumably mediated by the hormonal system (Muñoz-Calvo & Argente, 2016), triggers a series of physiological events to maintain body functioning, being slow heart rate (HR) and increased HR variability the most noticeables (Bär et al., 2006;Mazurak et al., 2011;Peyser et al., 2021). How the sympathovagal balance changes during the disease is still under debate, as indicators of HR variability do not unequivocally indicate an increase in parasympathetic outflow (Peyser et al., 2021).
Although there are some inconsistencies in the literature, the majority of research suggests that AN patients have increased HR variability and reduced HR compared to healthy individuals (Bär et al., 2006;Melanson et al., 2004;Spaulding-Barclay et al., 2016). A few studies have reported no differences in frequency or time domain cardiac indices between AN patients and healthy subjects (Melanson et al., 2004;Vigo et al., 2008), while others found reduced HR variability in patients (Lachish et al., 2009). Such discrepancies across studies are likely attributable to sample heterogeneity such as illness state (i.e., acute or recovered AN) or duration, medication, psychiatric comorbidities or methodological issues in HR recordings. Importantly, studies investigating different eating disorders found no differences in autonomic cardiac measures, suggesting that central regulatory mechanisms are directly responsible for changes in HR and its variability and not weight fluctuations or malnutrition (Mazurak et al., 2011).
Investigations of brain-heart interaction has been mainly oriented towards identifying the pathways by which the central control of the heart is achieved (Valenza et al., 2017(Valenza et al., , 2019. Benarroch (1993), in his pioneer work, described the so-called central autonomic network (CAN) as responsible for monitoring and regulating the peripheral organs via sympathetic and parasympathetic outflow. The central CAN regions involved in cardiac regulation include cortical and subcortical limbic regions, such as the ventromedial prefrontal cortex (vmPFC), anterior cingulate cortex (ACC), amygdala and anterior insula (Beissner et al., 2013;Thayer & Lane, 2000). Several neuroimaging studies have shown that variations in cardiac cycle covariate with restingstate functional connectivity (FC) changes in these regions (de la Cruz et al., 2019;Iacovella & Hasson, 2011;Valenza et al., 2019). In the context of AN, although a large body of literature exists on FC, none has paid attention to the possible relation between autonomic imbalance and alteration in CAN connectivity. The lack of CAN studies can be attributed to the diffuse anatomical boundaries of the CAN and the challenges of assessing brain-heart interaction in fMRI due to the influence of physiological confounds. In addition, the high heterogeneity of FC methods mainly focusing on whole-brain approaches (Gaudio et al., 2016(Gaudio et al., , 2018Geisler et al., 2020;Lotter et al., 2021) has limited the formulation of specific hypotheses regarding biological mechanisms underlying FC alterations in AN, such as whether there exists an association between connectivity of brain regions involved in HR control and the cardiac autonomic dysfunctions reported in AN.
Neuroimaging studies on dynamic functional brain connectivity in AN are even more scarce. Contrary to classical FC analysis, dynamic FC allows to monitor dynamic fluctuations of brain connections instead of a static measure of them (Hutchison et al., 2013). There is an overall consensus that dynamic changes in biological systems (e.g., brain and heart) are complex, characterized by nonlinearity and intrincate interactions (Francesco et al., 2012;Goldberger et al., 1990;Isaeva, 2012;Motta & Pappalardo, 2013). A large body of research studies suggests that complexity of both brain and heart signals changes with ageing or disease (Blasco-Lafarga et al., 2010;Jia et al., 2017;Lin et al., 2019;Lipsitz & Goldberger, 1992;Sokunbi et al., 2014;Vaillancourt & Newell, 2002). In particular, the resting-state brain signal in psychiatric disorders, such as in schizophrenia or depression, shows higher complexity compared to that of healthy subjects likely attributable to a dysregulation of the nonlinear dynamics of underlying neuronal systems (Sokunbi et al., 2014). The quantification of biosignals' complexity can be achieved through sample entropy (SampEn), a statistical measure of the irregularity in data sequence (Richman & Moorman, 2000). Unlike other indices assessing complexity, SampEn requires a shorter data length, which makes this measure very suitable for analysing functional magnetic resonance imaging (fMRI) data. The value of SampEn of short time series cardiac data appears primarily determined by parasympathetic modulation (Beckers et al., 2006;Jelinek et al., 2017; and reflects the level of adaptability of the cardiac control mechanism to internal and external stress (Ucak et al., 2021). Given the interaction of brain and heart functions (Azzalini et al., 2019), we should expect some degree of correlation in the degree of complexity of physiological outputs characterizing both organs, such as between FC and HR time series.
In this exploratory study, we aim to investigate FC alterations in the CAN in patients with AN. To this end, we collected resting-state fMRI data in 21 acute AN patients and 24 healthy controls (HCs) with similar HR in order to discard potential differences that may arise due to different HRs (de la Cruz et al., 2019). Based on previous reports of autonomic imbalance (Peyser et al., 2021) and FC abnormalities in brain areas putatively involved in cardiac regulation, such as insula, amygdala or ACC (Gaudio et al., 2016), we hypothesized reduced FC of CAN regions in AN individuals. Consistent with the theory that physiological signals in individuals with psychiatric disorders are more complex than those in HC (Lin et al., 2019;Sokunbi et al., 2014;Vaillancourt & Newell, 2002), we also hypothesized that the degree of complexity of CAN connectivity dynamics differs between HC and AN subjects, with higher complexity in AN patients. In addition, given the close link between HR or its variability and FC of CAN regions de la Cruz et al., 2019), we examined the association between the degree of complexity in CAN connectivity dynamic and HR time series. Finally, we aimed to identify and characterize discrete FC patterns or CAN states and their relation to HR complexity.

| Participants
We recruited 19 female patients with AN and 24 age, HRand gender-matched HCs. All AN patients met the DSM-IV criteria for AN according to the Structured Clinical Interview for DSM-IV Axis I disorders (SCID-I; First et al., 2002). General psychopathology was assessed using the Eating Disorder Inventory-2 (EDI-2; Thiel et al., 1997), the State-Trait Anxiety Inventory (STAI trait, STAI state; Laux et al., 1981) and the Beck Depression Inventory (BDI-2) self-report questionnaire (Beck et al., 1996). Patient's study started within the 3rd and 7th days after admission to hospital, allowing patients to acclimate to the hospital routine and reduces interference from any type of interventions.
All participants were right-handed persons as evaluated by the modified version of Annett's handedness inventory (Briggs & Nebes, 1975) with age range of 18 to 40 years (one patients with 40 years and the rest of participants 18 to 32 years). Participants were thoroughly examined by means of a clinical examination, standard electrocardiogram, routine laboratory investigations and screened for a history of other disorders such as major depression, personality disorders or obsessivecompulsive disorder. According to the SCID-I interview, none of the HC had a current episode or a history of mental disorder.
For AN, an expert clinician confirmed from medical records no comorbid psychiatric diagnoses other than eating disorders. Six AN participants were taking antidepressant medication at the time of the study. Informed written consent was obtained from all participants, and the study was carried out following the protocols approved by the local Ethics Committee of the University Hospital Jena and according to the guidelines of the Helsinki Declaration (from 2013). Table 1 shows a description of participants' demography, physiology and questionnaire variables.

| Image acquisition and physiological recordings
MRI data were collected on a 3T whole body-system equipped with a 12-element head matrix coil (MAGNETOM Prisma, Siemens Healthcare, Erlangen, Germany). Participants were instructed to keep their eyes open during the entire measurement. We acquired T 2 *weighted images using a multiband multislice GE-EPI sequence with the following parameters: repetition time (TR) = 484 ms, echo time (TE) = 30 ms, flip angle (FA) = 90 , multiband factor = 8, matrix size = 78 Â 78 with in-plane resolution of 2.5 Â 2.5 mm 2 and 56 contiguous transverse slices of 2.5 mm thickness covering the entire brain. The session lasted approximately 15 min and contained a series of 1,900 whole-brain volume sets.
Cardiac and respiratory activities were recorded (500 Hz) during rsfMRI data acquisition using an MRcompatible BIOPAC MP150 polygraph (BIOPAC Systems Inc., Goleta, CA, USA). Respiratory activity was measured by a strain gauge transducer incorporated in a belt tied around the chest at the level of the processus xiphoideus. The cardiac signal was recorded using a standard photoplethysmograph (PPG) attached to the proximal phalanx of the index finger of the subject's left hand. The PPG can be used as a surrogate technique for the electrocardiogram and is usually preferred for cardiac recordings since electrocardiogram-derived signals exhibit greater susceptibility to electromagnetic and biologic interference (Birkholz et al., 2004). To remove MRI-related or movement artefacts PPG and respiratory signals were digitally filtered at .05-3 Hz and 0-10 Hz, respectively.

| rsfMRI preprocessing
Data preprocessing was performed using the 'afni_proc. py' script in the AFNI software package (https://afni. nimh.nih.gov/; Cox, 1996). After discarding the first 20 volumes, time courses were despiked and artefacts time-locked to the cardiac/respiratory cycles, and slow blood oxygenation level fluctuations modelled via RET-ROICOR (Glover et al., 2000) and respiration volumes per time (RVT) regressors (Birn et al., 2008;Chang et al., 2009). The model consisted of eight RETROICOR regressors (four respiratory and four cardiac) and five RVT functions delayed at 0, 5, 10, 15 and 20 s (Birn et al., 2008), generated on a slice-wise basis by AFNI's RetroTS.m script (Jo et al., 2010). Next, data were motion corrected, aligned to the anatomical scan, warped to the Montreal Neurological Institute (MNI) template and smoothed with a 6-mm full-width half-maximum (FWHM) Gaussian kernel. Finally, we filtered the data to retain frequencies in the .01-.1 Hz frequency band and reduced non-neural sources by regressing the following nuisance variables: (1) 12 motion regressors (6 realignment parameters and their derivatives), (2) voxelwise local white matter regressors and (3) 3 principal components of ventricle signals (ANATICOR; Jo et al., 2010). Masks for white matter and ventricles were generated from each participant's anatomical scan using Freesurfer 7.1.0 (http://surfer.nmr.mgh.harvard.edu). There were no differences in head motion between AN patients and T A B L E 1 Demographics and clinical data. healthy individuals as measured by mean framewise displacement ( p > .6 Mann-Whitney U-test).

| Static FC
We focused on six core regions of the CAN, namely, vmPFC, dACC, left and right amygdala and left and right anterior insula (aINS). Details on the definition of these seed regions are given in our previous study (de la Cruz et al., 2019). Briefly, we drawn the vmPFC and dACC seeds as a sphere of 10 mm radius, centred at MNI-coordinates, x = 0, y = 44, z = 14 and x = 5, y = 32, z = 36, respectively. For definition of left and right aINS and amygdala seed regions, we used the WFU Pick Atlas tool for SPM (Maldjian et al., 2003). The MNI-coordinates of the centre of gravity for left and right aINS were x = À33, y = 15, z = À2 and x = 38, y = 16, z = À4 and for left and right amygdala were x = À23, y = À2, z = À19 and x = 28, y = À1, z = À19, respectively. Figure 1 shows the locations of all six seed regions.
For FC analysis, we first extracted the average time series of each seed region and correlated via Pearson's correlation with all seed regions to generate a 6 Â 6 symetric FC matrix for each subject. Entries below or above the main diagonal represent the correlation strength between each pair of seeds. FC values were finally transformed to z score using a Fisher's Z transformation to yield a more normally distributed data.

| Dynamic FC
We adopt a sliding-window approach to estimate temporal changes in FC. For this analysis, scans were split into 60 s sliding time windows, with the onset of each window progressively shifted by 20 s (40 TR) from that of the previous window, resulting in a total of 43 windows. This window size of 60 s and overlapping of 40 s  results in a good trade-off between the sensitivity of detecting potentially interesting transients in FC and the signal-to-noise ratio of the estimated FC (Hutchison et al., 2013). Similar to the static FC, for each window, we calculated a 6 Â 6 FC matrix to reflect the connectivities between all pairs of regions. This analysis yielded dynamic FC matrices for each subject. From these FC matrices, we then extracted FC time series between each pair of regions.
We applied k-means clustering in order to identify a set of CAN states or recurring dynamic FC patterns. All dynamic FC matrices were concatenated across all participants, and k-means clustering was performed to classify each FC matrix using squared Euclidean distance as the cost function. In order to determine the optimal number of clusters, we used the elbow method by varying the cluster number from 2 to 20. This procedure consistently yielded five cluster centroids or CAN states, which were then used to classify the FC matrices for each subject. Next, we computed the number of state transitions and time spent in each state during the resting-state scan of each participant.

| SampEn of dynamic FC
We computed SampEn of FC time series to answer the hypothesis that the degree of complexity of CAN connectivity dynamics differs between healthy and AN patients. SampEn is a modification of approximate entropy, used for assessing complexity or irregularity of physiological time series. Contrary to approximate entropy, SampEn does not have self-counting and is mostly independent of the length of the series (Delgado-Bonal & Marshak, 2019). SampEn is the negative logarithm of the conditional probability that two similar sequences of m points remain similar at the next point m + 1, where F I G U R E 1 Locations of seed regions used to define the central autonomic network. Red, ventromedial prefrontal cortex (vmPFC); yellow, dorsal anterior cingulate cortex (dACC); green, anterior insula (aINS); blue, amygdala (AMYG).
larger SampEn values correspond to more complexity or irregularity in the time series. For a time series {x 1 , x 2 , … x N }, SampEn is defined by where C m r ð Þ indicates the probability of satisfying jx m i À x m j j < r for i ≠ j i, j ¼ 1, 2, … N-m. The quantities x m i , x m j represent the elements of an m-dimensional template vector x i , x iþ1 , …, x iþmÀ1 f g , and r is a distance threshold or tolerance value. In line with previous works (Miši c et al., 2010;Omidvarnia et al., 2021;Wang et al., 2018), we set m = 2 and r = .5. Along with SampEn of FC time series (SampEn FC ), we also computed the SampEn of HR time series (SampEn HR ) using the same window size of 60 and 40 s overlapping between windows. HR time series was defined by computing one average HR value for each window. We also computed a SampEn value at each window of the raw HR signal. In this way, we obtained time series of SampEn HR and CAN states during the resting-state scan for each participant. All SampEn computations were performed using the pracma R package.

| Statistical analysis
We conducted group-comparisons of FC and SampEn FC in two ways. First, individual entries below the main diagonal were statistically compared using a Student's t test to reveal specific differences between pair of regions in the CAN. Each p value assigned to each entry in the FC or SampEn FC matrix was corrected (Bonferroni correction) for multiple comparison at α < .05/15 (15 is the number of pair of regions). Second, for each group, we built a 1D-array with all entries below the main diagonal from the individual FC matrix. The two 1D-arrays, one for AN (15 Â 21 rows) and one for HC (15 Â 24 rows), were then compared using a Student's t test. In contrast to the first approach, the second approach aims to estimate group-differences in CAN connectivity as a whole and not on specific pair of seed regions. Before using Student's t test, we tested the normality assumption with the Shapiro-Wilk's test, and when violation, we used the Mann-Whitney U-test instead.
We performed an analysis of covariance (ANCOVA) to assess SampEn HR by diagnosis group (AN, HC) interaction (SampEn FC $ SampEn HR * diagnosis). This analysis is motivated by previous findings that static and dynamic CAN connectivity correlates with HR or its variability (Chang et al., 2013; de la Cruz et al., 2019), and the fact that the autonomic balance appears disrupted in AN (Bär et al., 2006), which may suggest an alteration in the brain-heart axis in AN.
To assess the influence of CAN states on SampEn HR , data were entered into a two-way mixed ANOVA model with state and group as within-and between-subject factors, respectively, and SampEn HR as the dependent variable. In order to simplify the ANOVA model, we did not enter the time series of SampEn HR and CAN states but instead computed an average SampEn HR for each CAN state. To achieve this, we averaged SampEn HR across all windows classified in a given CAN state for each participant.  Figure 2c. An interesting point to note was that connectivities among the investigated regions were always positive, strengthening the notion that they make up a functional network. The uncorrected p-value results in Figure 2b indicate important trends in the data for potential future investigation.

| Static FC
After verification of linear model assumptions, we conducted a post hoc linear regression analysis between the subject's average CAN connectivity in AN and HR variability (RMSSD) in order to explain the reduced FC in AN in terms of their autonomic imbalance. However, this analysis did not yield a significant relationship between these two variables (r = À.26, p > .2).

| Dynamic FC complexity
For each pair of seed regions depicted in Figure 1, we quantified the SampEn of FC time series (SampEn FC ). Figure 3 shows the average SampEn FC for control and patient groups (a), uncorrected p values for each F I G U R E 2 (a) Group average functional connectivity matrices showing pairwise correlations for all channels. (b) Uncorrected p value for between-group differences in functional connectivity in each channel and (c) between-group difference in functional connectivity for all channels together. Error bars represent 95% confidence intervals. Abbreviations: aINS_L, left anterior insula; aINS_R, right anterior insula; AMYG_L, left amygdala; AMYG_R, right amygdala; AN, anorexia nervosa; dACC, dorsolateral prefrontal cortex; FC, functional connectivity; HC, healthy control; vmPFC, ventromedial prefrontal cortex.
F I G U R E 3 (a) Group average matrices showing pairwise SampEn FC for each channels. (b) Uncorrected p value for between-group differences in SampEn FC in each channel and (c) between-group difference in SampEn FC for all channels together. * in (b) indicates survival of FDR correction. Error bars represent 95% confidence intervals. Abbreviations: aINS_L, left anterior insula; aINS_R, right anterior insula; AMYG_L, left amygdala; AMYG_R, right amygdala; AN, anorexia nervosa; dACC, dorsolateral prefrontal cortex; FC, functional connectivity; HC, healthy control; vmPFC, ventromedial prefrontal cortex.
F I G U R E 4 Group differences in SampEn HR (a), linear regression between SampEn FC and SampEn HR (b) and Spearman's correlation coefficients between FC and HR time series (c). The p values in (c) refer to one-sample t test comparing Spearman's correlation coefficients to 0. In a and c, violin plots represent the data distribution, the box indicates the interquartile range (from 25% to 75%), the line within the box represents the median and the triangle indicates the mean. Abbreviations: AN, anorexia nervosa; HC, healthy control; HR, heart rate; SampEn FC and SampEn HR , sample entropy of FC and HR time series. channels (b) and between-groups comparison for all channels together (c). Similar to static FC results, we found a significant difference between groups in Sam-pEn FC in the pair dACC-right anterior insula (SampEn FC-HC = .80 ± .22, SampEn FC-AN = 1.00 ± .19; t[41] = 3.10, p = .003) but now surviving the FDR correction. The group comparison showed that the overall SampEn FC in CAN is significantly higher (t[643] = 2.64, p = .009) in AN patients compared to HC, as shown in Figure 3c.

| Association of FC complexity with HR complexity
As shown in Figure 4a, we found a trend for group differences in SampEn HR (median HC = .65, median AN = .85, Mann-Whitney U = 299.0, p = .08 two-sided), a result that is surprising because groups were previously matched by HR. The ANCOVA analysis yielded a significant interaction effect (F[1, 39] = 5.7, p = .02) between SampEn HR and diagnosis groups. To visualize this interaction, we linearly fit the raw data as shown in Figure 4b. This analysis indicated a significant positive linear relationship between SampEn FC and SampEn HR in the control group (r = .62, p = .001) but not in AN patients (r = À.08, p > .7). To further examine differences in the brain-heart axis between groups, we performed a post hoc analysis correlating FC and HR time series of each subject via Spearman correlation. We then performed a one-sample t test to examine how the correlation values in each group differed from 0 and a twosample t test to compare the groups. As shown in Figure 4c, the correlation between dynamic changes in FC and HR did not differ between groups (t[41] = 1.7, p > .1 two-sided), yet there is a clear trend for higher coupling in HCs (mean Spearman's r = À.11; t[23] = À1.58, p = .064 one-sided) not present in AN patients (mean Spearman's r = .04; t[18] = .83, p > .2 one-sided).
The k-means clustering yielded five stables CAN states, as shown in Figure 5. States 1-4 show moderate to high correlation among the sub-cortical CAN regions, that is, anterior insula and amygdala. In particular, states 2 and 4 distinguish themselves from the other states in the positive connectivity among all channels. State 1 is similar to state 3 except that the connectivity of dACC with subcortical CAN regions is moderate or high. State 4 closely resembles the average FC matrices (see Figure 3a) of healthy (r = .96) and AN (r = .94) groups.
F I G U R E 5 Representative connectivity matrices (cluster centroids) of CAN states 1-5. The numbers in parentheses are the time spent in the state for all participants in the healthy (HC) and AN groups. The plot in the right bottom row indicates the interaction between states and groups on SampEn HR (F[4, 80] = 2.74, p = .034). Post hoc analysis of simple main effects revealed that state 5 drives the state x group interaction. When the CAN is in state 5, the SampEn HR of healthy subjects reaches its minimum value, while the SampEn HR of AN patients reaches its maximum value. *p < .01 corrected.
State 5 shows overall the weakest connectivity among all CAN states; specifically, both amygdalae and vmPFC regions appear weakly connected with the anterior insula and dACC. Most strikingly, we found an interaction effect between state and group on SampEn HR (F[4, 80] = 2.74, p = .034) driven by state 5, as revealed by a post hoc analysis of simple main effects (see Figure 5, right bottom row). This analysis showed that SampEn HR reaches its minimum and maximum values in HCs and AN patients, respectively, when the CAN is in state 5. Concerning the number of transitions and time spent in a given state, we found that groups did not differ in these metrics. In addition, subjects consumed a comparable time across CAN states irrespective of the group membership, implying no favourite state existed.

| DISCUSSION
To the best of our knowledge, this is the first study exploring the degree of complexity of FC and HR time series in acute AN and age-and HR-matched HCs. AN patients showed FC alterations in core regions of the CAN engaged in cardiac regulation. The overall FC between the six investigated seeds, namely, vmPFC, dACC, left and right insula and left and right amygdala, is reduced in AN patients compared to HCs. Regarding SampEn of FC time series (SampEn FC ), patients exhibited higher values at the overall level and in the channel dACC-right anterior insula, indicating higher complexity of the CAN connectivity dynamic. The SampEn of HR time series (SampEn HR ) did not differ between AN patients and controls, but it shows a strong relationship with SampEn FC in HCs. Interestingly, this relationship is not present in AN. Finally, we found that the CAN transits between five functional states with no particular preference for any of them. However, in state 5, the SampEn HR significantly differed between groups, reaching its minimum value in controls and maximum in AN patients.
Our results of reduced FC in AN are consistent with previous studies investigating FC of large-scale resting state networks. For example, Gaudio et al. (2015) reported reduced FC between the executive control network and the ACC, likely attributable to the impaired cognitive flexibility in relation to body image and appetite in AN individuals (Gaudio et al., 2015). In another functional study, Phillipou et al. (2016) found decreased FC between the sensorimotor and visual networks in AN individuals compared to controls (Phillipou et al., 2016). They similarly explained their findings in terms of the altered visuospatial processing and the body image distortion experienced by AN individuals. Previous neuroimaging studies have reported an association of increased parasympathetic modulation with increased FC between default-mode and sensorimotor networks with CAN regions (Dhond et al., 2008;Schumann et al., 2021). In this study, we did not find a significant association between FC and the parasympathetic index RMSSD. A possible explanation is a shift from central to more peripheral influences on the heart in AN, which may lead to uncoupling of the central regions involved in cardiac regulation. An argument favouring this hypothesis was the uncoupling found in AN patients between the degree of complexity of FC and HR time series (see Figure 4b).
The use of our seed regions, that is, vmPFC, amygdala, anterior insula and dACC, is justified by neuroimaging studies examining neural correlates of autonomic control and by their role in the disorder (Beissner et al., 2013;Gaudio et al., 2016;Kaye et al., 2013;Valenza et al., 2019). Valenza et al. (2019) previously revealed that they make up the core of the neural substrate sustaining cardiovascular oscillations at rest. In a recent study of our group, we showed that FC between these CAN regions covaries with HR (de la Cruz et al., 2019), giving more evidence for their involvement in cardiac regulation. We should bear in mind that CAN regions are also implicated in emotion regulation (Mather & Thayer, 2018), a domain that is particularly affected in AN (Harrison et al., 2009). Thus, our findings of disrupted connectivity and entropy may be possibly reflecting at the same time the emotional processing deficit and speculated brain-heart uncoupling present in AN. On the other hand, the anterior insula has an essential role in integrating interoceptive information (e.g., internal physical sensations like taste and hunger). Many of the symptoms of AN, such as distorted body image, failure to appropriately respond to hunger and diminished motivation to change, are related to disturbed interoceptive awareness due to altered insula activity (Pollatos et al., 2008).
In line with our hypothesis, AN patients have higher complexity in FC time courses than HC despite having overall reduced CAN connectivity. This finding has also been observed in other psychiatric disorders and suggests a dysregulation of the nonlinear dynamics of underlying neuronal systems (Sokunbi et al., 2014). Indeed, an intriguing hypothesis poses that system complexity increases with illness (Vaillancourt & Newell, 2002). Thinking in the neurobiology of AN, the abnormally higher complexity may be related to an abnormal neurochemical mechanism. Accumulating evidence suggests that alteration in serotonin and dopamine systems in AN causes a skew towards aversive or inhibitory responses, which may underlie behavioural traits observed in the disorder, such as inhibition, inflexibility and anxiety (Kaye et al., 2013). Dopamine (DA) and serotonin (5-HT) are two monoamine neurotransmitters produced primarily in brainstem regions that help regulate many bodily functions. These brainstem neurotransmitters nuclei have widespread projections to several subcortical and cortical areas, including vmPFC, amygdala, insula and dACC (Conio et al., 2020). A very recent study using a mouse model of AN discovered that the dopamineserotonin circuit is super activated in AN, when compared to controls, and responsible for the lack of appetite and excessive exercising seen in the animals (Cai et al., 2022). Although very speculative, it is possible that the higher complexity of the FC time series observed in AN patients may be a consequence of a dysregulation in the serotonin and dopamine activities, which does not allow to keep the system stable. Another interesting hypothesis is that metabolic changes as the result of malnutrition drive higher entropy in AN patients. The welldocumented alterations in glucose metabolism in AN (Delvenne et al., 1995;Phillipou et al., 2014) may effect the BOLD signal and, henceforth, FC dynamic.
We did not find significant differences in SampEn HR between AN and HC. This finding contrasts with previous studies showing altered complexity in HR series in AN patients (Mazurak et al., 2011;Vigo et al., 2008). Yet we should note that some studies used another measure of entropy, that is, approximated entropy, which is inherently biassed due to self-matches, heavily dependent on the data length and less consistent than SampEn (Chen et al., 2005). Additionally, since we matched groups by HR, it is expected not to find differences in SampEn HR . However, we must stress that the control group has an unusual lower HR (62 bpm) than the average healthy population of about 70 bpm, and therefore, we cannot discard differences when compared to another sample of healthy individuals. Moreover, the degree of complexity of HR series in AN depends on the stage of the disease. It has been shown that in anorexic patients, the change varies from increased complexity in the acute phase to loss of complexity in the chronic phase .
We found a strong relationship between the Sam-pEn FC and SampEn HR in HC that is not present in AN. As stated above, a possible explanation is that HR in anorexia is more influenced by peripheral signals than by central commands. This leads to the interesting speculation that the uncoupling brain-heart assessed by Sam-pEn FC and SampEn HR can be used as a biomarker of diagnosis and disease status. This uncoupling as a biomarker may be more robust than HR or its variability, which findings are inconsistent across the AN literature (Peyser et al., 2021) and do not directly account for the brain-heart axis. Yet we acknowledge that the cross-sectional nature of our study does not allow us to differentiate between state or trait markers and make a call for further investigation of SampEn in longitudinal studies and large data sets. In addition, we suggest that future research focuses on assessing correlation between brain and cardiac measures in AN individuals because the absence of a group effect found here (see Figure 4c) cannot be taken as final conclusion given the low power of this study.
We identified and characterized recurrent patterns of CAN connectivity dynamics for the first time. Interestingly, despite the CAN transits between five discrete states, there is no preference for any of them, regardless of group memberships. Even more strikingly, the HR complexity significantly differed between groups in the state characterized by the weakest connectivity among the CAN regions, despite groups being matched by HR. An outcome of this weaker coupling is a diminished tonic inhibition of the amygdala by the vmPFC, which might result in changes in brain-heart integration (Thayer et al., 2012) and, with this, altered HR entropy. However, why the HR entropy of healthy and AN patients diverges at the state of weakest connectivity is intriguing and deserves further investigation. A possible explanation is that the already disrupted autonomic regulation in AN patients becomes less stable at this state, leading to higher complexity in the HR dynamic.

| LIMITATIONS
Although this study is the first to examine the degree of complexity of CAN connectivity and its relation with HR complexity in AN patients, it is not without limitations. First, our study sample is relatively small. A larger sample size would be particularly important to power the statistical analyses conducted to determine group differences in CAN connectivity and SampEn. Second, the cross-sectional nature of the study did not allow for assessment of FC abnormalities occurring with changes in disease severity and, therefore, whether the observed findings are an effect of acute starvation or rather something 'trait' like or 'scar' like (Lotter et al., 2021). Third, the exclusive presence of female participants raises the question whether the results are generalizable to male AN patients. Moreover, although we controlled for age and BMI, we cannot exclude the influence of other potential confounders not accounted for in our results, such as the menstrual cycle (Bai et al., 2009;Petersen et al., 2014) and physical exercise (Tulppo et al., 2001). Furthermore, the amount and type of seed regions used to define the CAN may also influence the findings. While we focused on six core CAN regions, there are other brain areas involved in cardiovascular control. For instance, Valenza et al. (2020) recently reported a set of brain regions functionally linked to complex cardiovascular oscillations that did not include the classical CAN regions involved in sympathovagal control. On the other hand, it seems that FC abnormalities in AN are distributed across many brain areas, which raises questions regarding the specificity of FC findings in small predefined networks (Lotter et al., 2021). Finally, the techniques used for dynamic analysis, that is, sliding window correlation and SampEn, are sensitive to their input parameters (Shakil et al., 2018;Yentes et al., 2013). Thus, we recommend interpreting these results with caution and suggest further research to assess the effects of different parameter settings and methodological approaches (Cabral et al., 2017;Eavani et al., 2013).

| CONCLUSIONS
In sum, our study shows for the first time that core regions of the CAN engaged in cardiac regulation are functionally affected in acute AN individuals. On average, AN patients have reduced FC in these regions, possibly due to a shift from central to more peripheral modulation of the cardiovascular system. The CAN transits between five discrete states with no preference for any, though at one, the HR complexity significantly differs between healthy and AN individuals. As observed in patients with other psychiatric disorders, AN individuals have higher entropy in surrogates of brain activity than healthy subjects. Moreover, AN patients show uncoupling of the brain-heart axis, as assessed by the correlation in the degree of complexity between FC and HR series, which could also be attributed to the shift in the underlying mechanisms subserving HR regulation. Whether the nonlinear FC dynamics alterations described herein have an adverse prognostic significance similar to that seen in other psychiatric diseases remains to be elucidated.

CONFLICT OF INTEREST STATEMENT
The authors report no biomedical financial interests or potential conflict of interest.