Nervous-System-Wise Functional Estimation of Directed Brain–Heart Interplay Through Microstate Occurrences

Background: The quantification of functional brain–heart interplay (BHI) through analysis of the dynamics of the central and autonomic nervous systems provides effective biomarkers for cognitive, emotional, and autonomic state changes. Several computational models have been proposed to estimate BHI, focusing on a single sensor, brain region, or frequency activity. However, no models currently provide a directional estimation of such interplay at the organ level. Objective: This study proposes an analysis framework to estimate BHI that quantifies the directional information flow between whole–brain and heartbeat dynamics. Methods: System–wise directed functional estimation is performed through an ad-hoc symbolic transfer entropy implementation, which leverages on EEG-derived microstate series and on partition of heart rate variability series. The proposed framework is validated on two different experimental datasets: the first investigates the cognitive workload performed through mental arithmetic and the second focuses on an autonomic maneuver using a cold pressor test (CPT). Results: The experimental results highlight a significant bidirectional increase in BHI during cognitive workload with respect to the preceding resting phase and a higher descending interplay during a CPT compared to the preceding rest and following recovery phases. These changes are not detected by the intrinsic self entropy of isolated cortical and heartbeat dynamics. Conclusion: This study corroborates the literature on the BHI phenomenon under these experimental conditions and the new perspective provides novel insights from an organ–level viewpoint. Significance: A system–wise perspective of the BHI phenomenon may provide new insights into physiological and pathological processes that may not be completely understood at a lower level/scale of analysis.

Within a brain-body framework, the set of complex interactions comprising anatomical, functional, biochemical, and bioelectrical CNS-ANS links are commonly referred to as functional brain-heart interplay (BHI). Modulations of the BHI have been correlated with physiological conditions, such as emotion perception [9], sleep [10], human body allostatic responses to external stimuli (e.g., autonomic maneuvers) [11], deep breathing [12], and cognitive load [13]). Moreover, dysfunctional BHI has been reported in several pathophysiological conditions [3], [6], [14], [15], including also schizophrenia [16], epilepsy [17], and mild depression [8]. An important aspect of BHI assessment is directionality. Indeed, while CNS actively influences the ANS, a different influence is concurrently exerted by the ANS to the CNS. For example, there is a higher neural modulation on heartbeat dynamics in subclinical depression [8]. Moreover, a cardiac sympathovagal initiation of a functional brain responds to emotional [9] or somatosensory [18] stimuli.
The quantification of functional BHI faces a number of technical issues from a purely methodological viewpoint. These issues include being intrinsically multimodal and multivariate, and being diffuse over the CNS and not specifically localized. Additionally, there is a directionality issue mentioned above and physiological plausibility that need to be considered when applying classical signal processing tools. Notwithstanding, several techniques have been applied or specifically developed to estimate the BHI. A synthetic data generation model was developed and exploited in different contexts, stressing the directionality of the BHI phenomenon [8], [9]. This model has also been investigated in conjunction with other methods [19], [20]. Information theory quantifiers were also developed to disentangle linear and nonlinear interactions [21], [22]. Moreover, measures such as point-process-based transfer entropy focused on instantaneous heartbeat response to scalp activity [11]. Furthermore, heartbeat-evoked potentials investigated the grandaverage of scalp response to heartbeats perceived as interoceptive stimuli [18]. Other studies have found that BHI extended to the multifractal domain [23].
Despite the significant results achieved by these methodologies, some aspects remain unclear owing to technical details. Specifically, functional BHI has mainly been identified at a single EEG-electrode level, thus quantifying the interplay between single-channel EEG oscillations and heartbeat oscillations, neglecting a nervous-system-wise level of interaction.
Indeed, a recent scientific consensus highlights brain activity as a whole, as well as embedded networks, as crucial way of functioning that sustains high-level neural processing [24], [25]. To this extent, in this study we focus on a whole-brain activity estimation through identification of EEG microstates [26], [27].
EEG microstate analysis has been widely employed to investigate the spatial and temporal properties of the brain in a non-invasive manner [26], [27], [28], [29]. It is possible to reconstruct the EEG amplitude topographical distribution on the scalp surface with a high resolution in time and space using multi-channel recording arrays and volume conductance models. Consequently, scalp activity can be visualized as a series of transient scalp electric-field maps [26], [28]. Each instantaneous map reflects the sum of momentarily active brain processes. Therefore, changes in the spatial configuration of the map imply that different neural elements performed an activation shift (i.e., passing from inactive to active, or vice-versa). However, this continuous scalp electrical activity is instead the concatenation of many building blocks (EEG microstates), which are quasistable topographic maps of electric potentials. On average, each segment remains stable for approximately 60−250 ms before rapidly transitioning to a different topography that remains stable again [26], [28]. Many reports have found correlations between microstate features that are consistent with the idea that EEG microstates may reflect the activity of brain networks, including frequency of occurrence, transition probabilities, average duration, various cognitive activities, behavioral states, and neuropsychiatric diseases. For example, microstates of certain topographies have a shorter average lifespan in schizophrenia [30], and a longer lifespan in panic disorder [31]. Microstate series vary with cognitive and behavioral states, such as drowsiness [32] and sleep stages [33].
Activity of several cortical and subcortical brain regions have been linked to microstate dynamics, with particular reference to the insular cortex, thalamus, amygdala, anterior cingulate cortex, and others [34], [35]. Remarkably, those regions are reportedly known belonging to the CAN [2], [4].
Here, we hypothesize that a causal, bidirectional, functional link occurs between brain microstates and heartbeat dynamics. Accordingly, we propose a novel methodology to quantify BHI at a nervous-system-wise level aiming to: i.) investigate the causal link between microstate series and cardiovascular dynamics; ii.) investigate functional BHI at a wholebrain level, going beyond the standard EEG-channel-specific approach. To this end, a directional symbolic transfer entropy (STE) analysis is developed to estimate the information transfer between a symbolic representation of heartbeat dynamics and a symbolic representation of whole-brain activity through EEG-microstates.

II. SIGNAL PREPROCESSING AND EXPERIMENTAL DATA
Two different experimental conditions eliciting concurrent CNS and ANS changes were analyzed to validate the proposed system-wise BHI analysis framework. The first experimental condition is associated with the cognitive workload performed through consecutive mental arithmetic calculations, and the second is associated with the cold pressor test (CPT), an autonomic maneuver that causes strong sympathovagal elicitation through thermal stress. This study received formal approval from the qualified ethical committee of the University of Pisa.

A. Mental Arithmetic
Experimental mental arithmetic approaches can activate the sympathetic nervous system through CNS manipulation. Participants are usually required to complete various cognitive tasks by frequently clicking a button or performing algebraic calculations within a certain amount of time [36]. Although mental arithmetic tasks have frequently been investigated at the CNS [37], [38] and ANS [39] levels, only a few studies have focused on their functional BHI correlates. In response to stress, variations in cardiac output correlate with neural activity in the left temporal and lateral frontal lobes [40]. Additionally, the BHI appears to increase in magnitude, and the information flows from the scalp's post-central and central regions to the heart appear to increase during mental arithmetic [41].
The first dataset (D1) was EEG During Mental Arithmetic Tasks [42], which is publicly available from the Physionet.org data repository (https://physionet.org/content/eegmat/1.0.0/). This dataset is comprised of EEG and ECG gatherings from 36 healthy volunteers undergoing a 180 s resting phase and a 60 s mental cognitive workload task (i.e., performing mental arithmetic). Recordings from four subjects were rejected after visual inspection owing to gross artifacts. Thus, data from 32 persons (24 females) with ages of 18 ± 2.01 years on average were retained for further processing. The electrophysiological signals were sampled at 500 Hz. The eligibility criteria were normal or corrected-to-normal visual acuity, normal color vision, no clinical manifestations of mental or cognitive impairment, and no learning disabilities. The use of psychoactive medication, drug or alcohol addiction, and psychiatric or neurological complaints were considered as additional exclusion criteria.
Power line notch (50 Hz) and [0.5 Hz−45 Hz] band-pass filters were applied in the EEG series before independent component analysis, and were used to identify artifacts (i.e., eyes, muscles, and cardiac pulsation) that were subsequently rejected. Further details on signal acquisition and preprocessing can be found in [42]. EEG signals were re-referenced in accordance with the REST method [43], as per recommendation from a previous microstate investigation [44]. To have series with equal length, we took the first segment of 60 s for each experimental condition.

B. Cold Pressor Test
Several studies have extensively investigated BHI changes during a CPT [11], [23], [45]. The CPT is a test for examining the body's autonomic functioning and CNS response to intense thermal and sub-pain-threshold stimuli [46], [47], [48]. A CPT activates physiological systems such as the baroreflex to maintain the body in a homeostatic state through the enhanced sympathetic activity of the ANS [49]. Pragmatically, this process typically entails submerging a distal limb (hand or foot) in cold water for 1 to 5 min. Numerous cortical and subcortical brain regions are involved in the brain correlates of a CPT, such as frontal regions in various frequency bands, posterior-parietal areas in the alpha band, and peripheral-bilateral temporal regions in the beta band [11], [48], [50]. Therefore, it is hypothesized that BHI estimates employing an organ-level whole-brain approach may be suitable for studying physiological responses to CPT. Furthermore, BHI findings have highlighted diffuse bidirectional interplay with stronger intervention of brain dynamics into the heartbeat activity [11], [45], [51].
The second dataset (D2) used in this study was collected from 30 healthy right-handed subjects (26.7 yrs on average; 15 males) who volunteered to participate in the experiment. The subjects were seated on a comfortable chair and performed an initial 3 min resting state followed by the actual cold pressor stimulation, which consisted of submerging their non-dominant hand into cold water that was maintained at approximately 4 • C.
Participants were asked to hold the position for up to 3 min, which is considered a time threshold that will not elicit pain perception [47]. However, participants were free to remove their hand if they felt uncomfortable. As a result, six participants did not reach 2 min of cold pressor stimulation and were excluded from further analysis (i.e., 24 subjects remaining). Physiological signals were recorded using a 128-electrode EEG and 1-lead ECG with a sampling frequency of 500 Hz. Researchers may obtain raw data through reasonable mail requests if ethical requirements are met.
An R-beat detection from the ECG series was applied using the well-known Pan-Tompkins algorithm [52], followed by an automated and visually inspected artifact rejection, all implemented in Kubios Software [53].
EEG signals were preprocessed through the Harvard Automated Processing Pipeline for Electroencephalography (HAPPE), which is extensively described in [54]. Then, the EEG signals were implemented through an EEGLAB toolbox in MATLAB software (MathWorks Inc.) [55]. Briefly, 57 more external channels were rejected, then a bandpass filter (between 1 and 100Hz) and notch filter (50 Hz) were applied. Furthermore, bad channels were removed after being identified as the most external 1% tail of a joint distribution based on high-order statistical moments. HAPPE implements a wavelet-enhanced independent component analysis to detect and remove periodic artifacts (e.g., eye blink, heartbeats, and respiration). Subsequently, an automated independent component analysis based algorithm was used to remove the remaining artifacts (e.g., muscular and motor) [54]. Rejected EEG channels were spherically interpolated and the final 71 channels were re-referenced in accordance with the REST method [43], as per recommendation from a previous microstate investigation [44]. To have series with equal length, we took the first segment of 2 min for each experimental condition.

III. SYSTEM-WISE BHI ANALYSIS FRAMEWORK
A comprehensive diagram is presented in Fig. 1 and a detailed description of the proposed system-wise BHI analysis framework is provided below.

A. EEG Microstate Analysis
Extensive details on the microstates identification from EEG signals are reported in [26], [29], and a MATLAB toolbox may also be used for their estimation [56].
Briefly, microstates are defined as EEG scalp topographies that do not change over a short time window (about 60− 250 ms) [26]. Consequently, microstates associated with different topographies alternate over time and may be task-and subject-specific. Identification of microstates from EEG signals is performed through non-parametric clustering, and each cluster denotes a microstate. A cluster is primarily assigned to peaks (i.e., local maxima) of the global field potential (GFP) [29], which is the standard deviation across all EEG electrodes calculated at each time sample [57]. In this study, the so-called modified k-means clustering algorithm was used [29], [56]. The dataset-specific number of microstates k was identified according to a metacriterion that accounts for different fit measures such as global explained variance, cross-validation criterion, Krzanovsky-Ly criterion, and variance [26]. A final visual inspection analysis is then performed to verify the physiological plausibility of the identified topographies, also known as prototypes. Once identified, each microstate prototype is assigned to each EEG sample based on the similarity between EEG scalp topography and prototypes. The reliability of this procedure, also known as back-fitting, is evaluated through a goodness-of-fit analysis based on the calculation of global explained variance [26]. Finally, microstate series is smoothed through the calculation of the mode in non-overlapping windows with 250 ms length to ensure continuity, i.e., minimizing rapid changes [56], and to match the timing of cardiovascular variability dynamics.

B. Heart Rate Variability Analysis
Heart Rate Variability (HRV) series were interpolated through a spline function to obtain a uniform sampling rate of 4 Hz that matches the microstates sampling rate. Then, the series were symbolized through a partition operation. A max − min transformation was applied to all the datasets, experimental conditions, and subjects to obtain sequences of symbols. The partition was defined as a quantization through five amplitude levels after the maximum and minimum values were determined within each signal. Each level was associated with a symbol between 1 and 5.

C. Symbolic Transfer Entropy
In order to estimate the system-wise functional BHI, a robust formulation of Symbolic Transfer Entropy (STE) based on permutation entropy [58], [59] was implemented.
For a given, but otherwise arbitrary i, m amplitude values X i = x(i), x(i + l), . . . , x(i + (m − 1)l) were arranged in ascending order x(i + (k i1 − 1)l) = x(i + (k i2 − 1)l) ≤ · · · ≤ x(i + (k im − 1)l), where l denotes the time delay, m is the embedding dimension, and {k i1 . . . k im } ∈ N represent the rearrangement indexes [58]. Consequently, a symbol was formalized asx i ≡ (k i1 , k i2 , . . . , k im ). After that, it was possible to estimate the joint and conditional probabilities of the sequence of permutation indices using the relative frequency of the symbols. In this study, the embedded dimension is set to m = 3.
Using symbol sequencesx i andŷ i , ST E Y →X is expressed as: where the sum runs over the entire symbolic series, δ indicates the time step (in this case δ = l = 250 ms), and the log has a base of 2. Thus, ST E Y →X is provided in bits. Moreover, ST E X→Y is specularly defined [58]. As suggested by [58], the directionality index was defined as D X Y = ST E Y →X − ST E X→Y . D X Y quantified the preferred information flow direction, assuming positive and negative values for net couplings where Y was mainly driving X and X was mainly driving Y , respectively. D X Y 0 is expected for symmetric bidirectional couplings. Note that D X Y is an index of directionality and not an index of coupling strength.
Additionally, ST E X→Y represents the strength of coupling from X to Y , which is not related to coupling from Y to X.
For comparison reasons, the information associated with each system (central nervous system and cardiovascular system) was estimated through the calculation of self entropy (SE) [59], [60], which quantifies the average reduction in uncertainty onx i+δ resulting from the knowledge ofx i [60]. Formally SE is defined as follows: where the sum runs over the entire symbolic series, δ indicates the time step (in this case δ = l = 250 ms), and the log has a base of 2.
The framework was implemented in MATLAB (Mathworks, Inc.), and the source code is available on a GitHub repository at https://github.com/CatramboneVincenzo/BHI_ SymbolicTransferEntropy.

D. Statistical Assessment and Comparisons
Functional brain-heart interplay estimates were employed to separately investigate significant changes between experimental conditions for the two datasets.
To this end, non-parametric statistics was employed, which did not make assumptions on the original distribution nature and was recognized to be robust against outliers [61]. D1 experimental sessions (resting state and arithmetic mental workload) were compared using a signed-rank test for paired samples, and D2 experimental sessions (resting state, CPT, and recovery) were compared groupwise using the nonparametric Friedman test for multiple paired samples [61], and pair-wise through a signed-rank test. Significance was set at 1 %.
Reliability of STE estimates was performed through a surrogate data analysis, with a null hypothesis of no causal interactions between systems [62].
Specifically, to test the reliability of an estimate, e.g., ST E Y →X , where Y is the driver system and X is the target system, we generated 500 synthetic series by random shuffling of the driver series. Consequently, while the surrogate driver series has the same statistical distribution as the original series, the target series is unaltered and so is its autocorrelation function. The distribution ofŜT E from the surrogate series constitutes the reference distribution of the null hypothesis, allowing for the calculation of p-value for the original statistic ST E Y →X . If the original estimate was more external than the 5% tail of the null distribution, the estimate was deemed significant and retained for further analyses. In our study, X and Y refer to brain and cardiovascular system through symbolic EEG-microstate and RR series, respectively, on both D1 and D2 datasets.

IV. RESULTS
The same system-wise analysis framework was applied for datasets D1 and D2.

A. Mental Arithmetic
The application of the microstate analysis to dataset D1 (mental arithmetic) led to the identification of three microstates for all subjects and experimental conditions. These microstates are illustrated in Fig. 2, also in terms of their temporal dynamics, along with an exemplary EEG series from a random subject (specifically, subject 31, first 5 seconds of the mental arithmetic task). Their topographical representations appear smooth over the scalp and can be associated with three main 'gradient directions'. The first represented an occipital-to-frontal gradient, whereas the second was in the opposite direction, introducing a slight left-to-right shift. The third microstate individuated a clear left-to-right gradient that was symmetric to the medial axis.
The experimental results for the two BHI directions (from heart-to-brain and brain-to-heart) and SE estimations (SE H and SE B ) are shown as boxplots with highlighted statistical differences in Fig. 3.
Sub-panels (a) and (b) represent SE boxplots for the EEGmicrostate (SE B ) and heartbeat (SE H ) series, respectively. Subpanels (c) and (d) represent ST E boxplots for heart-tobrain (ST E H→B ) and brain-to-heart (ST E B→H ) information transfer, respectively. Note that 24 out of 32 estimates resulted statistically significant in the surrogate data analysis on both ST E H→B and ST E B→H . Notably, SE measures did not significantly changed between the two experimental phases in terms of central tendency.
Evident differences have been found between the two experimental conditions in both directions considering BHI estimations, as shown in Fig. 3(c) and (d). More specifically, the statistical comparison between mental workload and the ST E H→B calculated during the resting phase resulted in a p-value of p = 2.7016e − 05. However, estimations resulted in a p-value of p = 1.8215e − 05 when applied to ST E B→H . Focusing on the absolute values of ST E estimations, ST E H→B achieved higher values than ST E B→H , indicating that the directional index D B H was significantly positive (Fig. 4). This quantifies the higher level of interplay in the heart-to-brain direction for both experimental conditions.
The differences between the two experimental conditions can be inferred from the results presented in Fig. 3. Specifically, the  elicitation provided by the cognitive workload implied appreciable difference in terms of the amount of information exchanged, or better interchanged, between the two systems bidirectionally. Such interchange, during the mental arithmetics task, become more predominantly heart-to-brain directed (Fig. 4).

B. Cold Pressor Test
The application of the microstate analysis to dataset D2 (CPT) led to the identification of five microstates for all subjects and experimental conditions. These microstates are illustrated in Fig. 5, also in terms of their temporal dynamics, along with an exemplary EEG series from a random subject (specifically, subject 30, first 5 seconds of resting state).
In this case, two additional microstates prototypes were identified with respect to dataset D1, namely the second and the fifth. This may be due to the higher level of variability introduced by the CPT dataset in terms of brain dynamics, particularly in the central area of the scalp, and on the right temporal lobe.
The experimental results for the heart-to-brain and brainto-heart BHI directions, as well as for SE H and SE B are summarized as boxplots with highlighted statistical differences in Fig. 6.
First, sub-panels (a) and (b) represent SE boxplots for the EEG microstate (SE B ) and heartbeat (SE H ) series, respectively. Second, sub-panels (c) and (d) display ST E information transfer boxplots for heart-to-brain (ST E H→B ) on the bottom left, and brain-to-heart (ST E B→H ) on the bottom right, respectively. Note that 20 out of 24 estimates resulted statistically significant in the surrogate data analysis on both ST E H→B and ST E B→H . Specifically, SE H and SE B did not report any significant group-wise differences between the experimental conditions (initial rest, CPT, and recovery) and the same was reported by the ST E H→B analysis. A single pair-wise significant comparison between CPT and recovery was detected in the SE H analysis, with SE H in CPT being significantly lower than during recovery.
Significant statistical differences were found among ST E B→H values extracted in the three experimental conditions (Friedman test p-values p = 0.0004) when focusing on BHI estimates through information transfer, as shown in Fig. 6(c). Interestingly, no differences were detected in ST E H→B , even if it reached higher values on average than ST E B→H . This was also reported for D B H values, which were generally positive, as shown in Fig. 7. This quantifies the higher level of interplay in the heart-to-brain direction for all the experimental conditions compared to the brain-to-heart direction. These results align with those found in dataset D1.

V. DISCUSSION
This study proposed a novel framework to estimate nervoussystem-wise directional functional BHI, which exploited microstates to summarize information for a whole-scalp EEG time series. The framework includes a uniform partition to represent HRV dynamics, and ST E to quantify directional information transfer between the brain and heart. Accordingly, particular attention was paid to the directionality of BHI, discovering that heart-to-brain interaction was a driving phenomenon compared to the brain-to-heart direction, in terms of absolute values. Indeed, previous studies on functional BHI highlights the importance of directionality assessment and may indicate a task-and state-dependent dominance for the brain-to-heart and heart-to-brain information flow [8], [9], [19], [21], [63], [64]. To illustrate, the functional heart-to-brain direction shows early changes during emotional processing in healthy subjects, therefore emotions may be associated with brain responses to peripheral changes to emotional stimuli [9]. Moreover, the functional brain-to-heart direction is dominant in resting state in subclinical depression [8].
Two different datasets were analyzed to validate the proposed model. The first dataset (D1) consisted of publicly available electrophysiological signals recorded during a cognitive workload task. The second dataset (D2) consists of strong sympathovagal and cerebral changes elicited through a CPT. A proof of concept of this research has been published in [65]. While findings of this study corroborate earlier research that highlights BHI occurring at single EEG channel level, they build further knowledge based on analysis at a whole-brain level. From a physiological viewpoint, results on mental arithmetics (dataset D1) are in agreement with previous findings [40], [66], [67], highlighting strong changes in directional functional BHI. In fact, no significant differences were found in EEG microstate or partitioned HRV series SE. However, directed ST E H→B and ST E B→H were significantly enhanced by the mental arithmetic task with respect to the preceding resting state. Furthermore, the directionality index D B H showed a positive sign, suggesting that heart-to-brain information transfer mainly drives such experimental conditions. First, scalp and heartbeat activity are known to be strongly affected by mental arithmetic in terms of power spectral density with respect to the resting state [40], [67], [68]. However, this does not necessarily imply a variation in the entropy of their activity. Second, mental arithmetics significantly enhance communication between the brain and heart in both directions, quantified here as ST E, which increases during mental calculation with respect to resting state, both in the ST E H→B and ST E B→H estimations. The authors believe that this may be due to the concurrent action of the cognitive and affective elicitation induced by the experimental task, which therefore might be associated with changes in heart-to-brain communications [9]. Note that cognitive stress reverberated together with sympathovagal activation in the limbic regions (e.g., amygdala), cortical regions (e.g., ventro-medial-prefrontal cortex), salience network, and brain areas and networks involved in CAN [68], [69], [70].
The experimental results on dataset D2 showed insignificant differences between experimental conditions (resting state, CPT, and recovery) in SE H and SE B , except for the CPT vs recovery SE H comparison. This was not surprising because previous studies on the same dataset needed to employ non-Gaussian expansion to find differences that were not detected in the linear analysis of HRV [71] and EEG [72] series. Differences in the experimental phases were detected in ST E B→H , which is consistent with previous studies [11], [45]. It should be noted that literature on BHI changes elicited by a CPT has also highlighted variations in the heart-to-brain direction, which was not found in the present study. This is not surprising because of the peculiarities of the methodologies introduced. Specifically, the ST E techniques allow for an organ-level quantification of the BHI and reasonably loses some partial region-or frequency-specific dependencies that are not spread over the scalp and spectrum. Focusing on the directionality of BHI in experiment D2, the index D B H showed a positive sign, even if no differences were enhanced among the experimental conditions. This confirmed a main heart-to-brain directional interplay.
Results referred to an embedding dimension m = 3. This dimension limits the curse of dimensionality and thus may be the best compromise between possible embedding dimensions, especially considering the different dynamics associated with heartbeat and EEG series.
There were some limitations in the present study, which focused on brain-heart interplay and neglected important autonomic covariates, such as respiratory activity and blood pressure. Moreover, ST E only measures part of the BHI dynamics; to illustrate, it does not reach the multifractal domain as the BHI is known to do [23]. The proposed surrogate data analysis employing random permutation of the driver series may constitute a sub-optimal evaluation; generation of time-shift surrogates or iterative-adjusted amplitude Fourier transform surrogates [73] should also be investigated used in future studies. Future research should also focus on applying the developed framework to other datasets that include healthy persons in various experimental settings (such as emotional elicitation) and subjects with pathological conditions, particularly those with mental or mood disorders.

VI. CONCLUSION
This study proposed a novel analysis framework to quantify the directional information flow between brain and heartbeat dynamics and investigated functional BHI. The proposed framework relied on information theory for STE estimation, microstate analysis to summarize EEG dynamics in a symbolic series of quasi-stable states, and statistical assessment through synthetic data generation. Two different datasets were utilized to validate the methodology. The first dataset focused on cognitive workload and the second investigated an autonomic maneuver as a CPT. The experimental results aligned with past research, highlighting the strong changes elicited by the experimental conditions with respect to the related baseline. From a purely methodological viewpoint, the proposed estimation of BHI provided novel insights into functional brain-heart communication. This includes: i) Employing microstate series to quantify BHI; ii) Providing directional estimates, thus quantifying information transfer from-heart-to-brain and from-brain-toheart separately; iii) Providing a directionality index, quantifying the main direction of information flow; iv) Being the first attempt estimating functional BHI at the organ level, thus going beyond the usual estimates at the single-sensor level; v) Statistically assessing the calculated measures against the null BHI hypothesis using surrogate data generation. The proposed technique provides new organ-level insights and perspectives on functional directional BHI and the framework can be used to represent factors related to various organs and systems in networks, which could be extremely useful in high-level studies such as network physiology.