In this study, two types of saccade task were used: the reflexive saccade task and voluntary saccade task. The reflexive saccade task contains Pro Saccade task (PS), and the voluntary saccade task contains Anti Saccade task (AS) and the Memory Guided Saccade task (MGS).

The driving factor for the PS task is primary exogenous, requiring no complex cognitive processes to be involved [32]. Specifically, in the PS task, a white cross appears at the center of the screen, and the participant must immediately fixate on the cross. When the cross disappears, a white dot randomly appears to the left or right of the fixation point. The participant must quickly and accurately direct their gaze towards the location of the white dot. The schematic of PS is shown in Fig. 1A.

Unlike the PS task, both the AS and MGS tasks are typically driven primarily by endogenous factors, involving the participation of higher cortical areas [33]. The AS task requires participants to suppress their automatic response to exogenous stimuli, thus effectively controlling their attention [34]. The MGS task, on the other hand, focuses on examining the memory function of participants [35].

In the AS task, the sequence of the white cross and white dot appearing on the screen is the same as in the PS task, but the participant must look at the opposite, mirrored location of the white dot after it appears (shown as Fig. 1B).

In the MGS task, when the white cross appears at the center of the screen, the participant must focus on the cross. During this fixation, a white dot randomly appears and disappears on the left or right side of the screen. The participant is asked to maintain focus on the central cross and remember the location of the appeared white dot. Once the cross disappears, the participant must immediately look at the location where the white dot appeared (shown as Fig. 1C).

A total of 30 young volunteers were recruited under three saccade tasks, during which both electroencephalogram (EEG), electrocardiogram (ECG) and eye movement signals were recorded simultaneously. Prior to the experiment, each participant was required to thoroughly clean their scalp with shampoo to remove any oils or residues. After this, the EEG cap (32-channel Quik-Cap Neo Net, Neuroscan, Compumedics Limited., Abbotsford, Victoria, Australia) was fitted, and conductive gel was applied to ensure that the impedance at each electrode site was below 10 k$\mathrm{\Omega }$. Subjects were seated in a dark, closed room, with their jaws placed on a special bracket for head immobilization, and their eye movements were recorded by an infrared vision-based video eye-tracker (EyeMind 2000, Jasmine Science and Technology Ltd., Chengdu, Sichuan, China) with a sampling rate of 1000 Hz. All visual stimuli were displayed on a 27-inch liquid crystal display (LCD, Dell LCD display, Dell Inc, Round Rock, TX, USA) monitor located 57 cm in front of the subject with a resolution of 1920 $\times$ 1080 pixels and a refresh rate of 100 Hz. The presentation of the visual stimuli and the acquisition of eye movement data were based on the Psychtoolbox toolkit (Psychtoolbox-3.0.19.10, https://psychtoolbox.org) for MATLAB (R2022b MathWorks, Inc., Natick, MA, USA), which was run on another computer. EEG and ECG signals were recorded simultaneously using NeuroScan EEG equipment (Neuroscan SynAmps RT, Neuroscan, Compumedics Limited.). Due to noise in the physiological data of four participants, their data were excluded from the analysis. Ultimately, data from 26 participants (aged 23–26 years, including 11 females and 15 males) were included.

This study was approved by the Ethics Committee of Hebei University Affiliated Hospital. Informed consent was obtained from all participants, who also confirmed that they had no history of neurological, cardiovascular, or other related diseases.

Before the official experiment began, participants read the experimental instructions and completed practice trials. Only after confirming their clear understanding of the instructions did the formal experiment commence. The task of the formal experiment was performed in different sessions, each consisting of 60 trials, 30 left and 30 right, with the order of appearance of the trials being randomized. Each session consisted of only one eye-tracking task. Before the start of each task, subjects were required to stare at nine visual stimulus points that appeared sequentially at different locations on the screen to complete the calibration and validation. Subjects were allowed to perform the task only when the error in eye position was less than that allowed by the validation; otherwise, they were required to perform the calibration and validation again. At the end of each session, subjects were required to take a 5-minute rest with their eyes closed before proceeding to the next session.

EEG data preprocessing was conducted using the EEGLAB toolbox (EEGLAB 2022.1, the Swartz Center for Computational Neuroscience (SCCN) at the University of California, San Diego, CA, USA, https://sccn.ucsd.edu/eeglab/index.php) and FieldTrip toolbox (FieldTrip toolbox 20240701, the Donders Institute for Brain, Cognition and Behaviour, Radboud University, Nijmegen, https://www.fieldtriptoolbox.org) in MATLAB. The preprocessing of EEG data involved the following steps:

(1) Channel Rejection: Bad channels were removed by inspecting the log power spectral density plot of all channels generated by EEGLAB.

(2) Filtering: A bandpass filter from 0.5 to 45 Hz was applied using the built-in finite impulse response (FIR) filter in EEGLAB.

(3) Re-referencing: The re-reference method used was the average of all electrodes across the brain [36].

(4) Artifact Removal: Independent Component Analysis (ICA), a function provided by EEGLAB, was used to remove non-brain-related components.

For power spectral analysis, we applied Short-Time Fourier Transform (STFT) with a Hanning window. A sliding window of 2 seconds (corresponding to 2 $\times$ 1000) was used with 50% overlap. The number of samples per window and overlap was chosen to ensure adequate spectral resolution, and to enable sufficient temporal coverage of task-related changes. Prior to spectral decomposition, the stationarity of EEG signals within each 2-second window was assessed by inspecting time-domain trends and visually inspecting spectrograms for abrupt shifts or trends. Non-stationary segments were excluded from further analysis. The Fast Fourier Transform (FFT) length for EEG data was set to the next power of two greater than the window length to enhance frequency resolution while avoiding excessive zero-padding that could distort time sensitivity. Subsequently, the time series were integrated across three frequency bands: below $\alpha$ (0.5–12 Hz), $\beta$ (12–30 Hz), and $\gamma$ (30–45 Hz). Furthermore, we evaluated the stationarity of the resulting EEG power spectral density (PSD) time series, which serve as inputs to the model using the Augmented Dickey-Fuller (ADF) test (Fig. 2A). The ADF test is a widely used statistical method for determining whether a time series is stationary by testing for the presence of a unit root.

The ECG time series were band-pass filtered between 0.5 and 45 Hz. A peak detection algorithm was used to identify the R-peaks of the ECG signal [36]. The detected R-peaks were visually inspected to identify any potential errors, which were then corrected accordingly. After detecting the heartbeats, we constructed the IHI sequence. To enable spectral analysis, the IHI sequence was uniformly resampled using spline interpolation to achieve a sampling rate of 4 Hz [37], corresponding to a temporal resolution of 0.25 seconds. This resampling facilitates time-frequency analysis. We then applied STFT with a Hanning window of 2 seconds (8 samples at 4 Hz) and 50% overlap. To verify the stationarity of the IHI signal within each 2-second window, we visually inspected the time series and spectrograms for abrupt shifts or drifts. Segments that clearly deviated from stationarity were excluded from analysis. For spectral decomposition of the ECG, the FFT length was set to 256 [38], which was selected to provide sufficient frequency resolution for observing LF (0.04–0.15 Hz) and HF (0.15–0.4 Hz) components. Furthermore, we evaluated the stationarity of the resulting ECG PSD time series, which serve as inputs to the model using the ADF test (Fig. 2B).

The SDG model evaluates the mutual interplay between brain and IHI dynamics across specific frequency bands (below alpha band (0.5–12 Hz), beta band (12–30 Hz), and gamma band (30–45 Hz)) in both low and high-frequency ranges [39].

To address the directional influence from the brain to heart activity, we adopted the Integral Pulse Frequency Modulation (IPFM) model, which conceptualizes heartbeat generation as a threshold-based integration process, wherein heartbeats are emitted once an integrated signal reaches a preset threshold. Specifically, this model was parameterized using the properties of the Poincare diagram [40], where the heartbeat was considered as a collection of impulse signals, and each impulse corresponds to a specific moment $t_k$. The heartbeat formation process builds up a characterization of the heartbeat signals by integrating between each heartbeat interval $t_k$ and $t_{k+1}$ until the result of this integration reaches one.

(1)$$X(t)=\sum _{K=1}^N\delta \left(t-t_k\right)$$

In this equation, $\delta \left(t-t_k\right)$represents the Dirac delta function marking the occurrence of heartbeats.

(2)$$1={\int }_{t_k}^{t_{k+1}}\left[{\mu }_{HR}+m(t)\right]𝑑t$$

Where ${\mu }_{HR}$ is the reference heart rate (Hz) and $m\left(t\right)$ is the interplay function. The interplay function $m\left(t\right)$ describes the variation of the inter-heartbeat period, which contains two oscillators related to sympathetic-vagal and parasympathetic activity, which influence heart rate variability (HRV). These oscillatory terms simulate how ANS inputs dynamically affect the generation of heartbeats through modulation of instantaneous heart rate.

(3)$$\mathrm{m}(\mathrm{t})={\mathrm{C}}_{\mathrm{LF}}(\mathrm{t})\times \sin \left({\omega }_{\mathrm{LF}}\times \mathrm{t}\right)+{\mathrm{C}}_{\mathrm{HF}}(\mathrm{t})\times \sin \left({\omega }_{\mathrm{HF}}\times \mathrm{t}\right)$$

Here, $C_{LF}$ and $C_{HF}$ are time-varying interplay constants calculated based on Poincare diagram features [40].

The angular frequencies ${\omega }_{LF}$ = 2$\pi \cdot$ 0.1 rad/s, ${\omega }_{HF}$ = 2$\pi \cdot$ 0.25 rad/s in which 0.1 Hz and 0.25 Hz correspond to LF and HF band central frequencies, respectively.

The interplay coefficient between the brain and the heart was calculated by the following equation:

(4)$$C_{LF}(t)=SD{G}_{\text{Brain}\to LF(t)}\times P_f(t-1)$$

(5)$${\mathrm{C}}_{\mathrm{HF}}(\mathrm{t})={\mathrm{SDG}}_{\text{Brain}\to \mathrm{HF}(\mathrm{t})}\times {\mathrm{P}}_{\mathrm{f}}(\mathrm{t}-1)$$

$P_f(t-1)$ is the EEG power at frequency f in the previous time window, and $SD{G}_{Brain\to LF\left(t\right)}$ is the interplay coefficient of brain to LF of heart. The negative value of $SD{G}_{Brain\to LF\left(t\right)}$ indicates that the activity of LF of heart decreases linearly with increasing brain, and vice versa. The level of $SD{G}_{Brain\to LF\left(t\right)}$ represents the strength of brain interplay of LF of heart. $SD{G}_{Brain\to HF\left(t\right)}$ is the interplay coefficient of brain to heart HF. A negative value of $SD{G}_{Brain\to HF\left(t\right)}$ indicates that the activity of heart HF decreases linearly with the increase of brain, and vice versa. The level of $SD{G}_{Brain\to HF\left(t\right)}$ represents the strength of brain interplay of heart HF.

Interplay from the heart to the brain was modeled by an adaptive Markov process to generate synthetic EEG sequences [41]. In this model, the upward interplay of the heart to the brain was realized through an autoregressive process and least squares. To be specific, this autoregressive process was designed to simulate how short-term HRV components in the LF and HF bands influence EEG power fluctuations across time. The model used the band and time window information of a specific EEG channel to integrate previous neural activity and current heartbeat dynamics to construct a Markovian neural activity sequence. In this process, the EEG signal EEG(t) can be expressed as a superposition of a series of frequency components with the following equation:

(6)$$EEG(t)=\sum _{f=f_1}^{f_n}{P}_f(t)\times \sin \left(2\pi ft+{\theta }_f\right)$$

$P_f(t)$ denotes the power at frequency $f$ at time $t$, while ${\theta }_f$ is the corresponding phase parameter.

The estimation of$P_f\left(t\right)$ follows the following autoregressive model:

(7)$$P_f(t)=k_f\times P_f(t-1)+{\mathrm{\Psi }}_f(t)+{\epsilon }_f$$

$k_f$ is a constant coefficient, reflect the temporal inertia of EEG power in each band. ${\mathrm{\Psi }}_f\left(t\right)$denotes the interplay of heart dynamics on brain-specific frequency band components, and ${\epsilon }_f$ is an error term used to adjust for deviations between the model and the actual data.

The heart to brain interplay coefficient was calculated by analyzing the effect of LF or HF on the generation of EEG data with the following equation:

(8)$$SD{G}_{LF\to Brain}(t)=\frac{{\mathrm{\Psi }}_f(t)}{IH{I}_{LF}(t)}$$

(9)$$SD{G}_{HF\to Brain}(t)=\frac{{\mathrm{\Psi }}_f(t)}{IH{I}_{HF}(t)}$$

$IH{I}_{LF}\left(t\right)$ and $IH{I}_{HF}\left(t\right)$ denote LF and HF IHI, respectively, which were obtained by analysing the time series of RR intervals in the ECG signals and calculating the frequency components using the STFT.

The rationale behind this modeling approach is to infer the directional influence of ANS activity (reflected in LF and HF components) on cortical dynamics. By quantifying this influence through normalized coupling coefficients, we aimed to capture time-resolved regulatory patterns.

$SD{G}_{LF\to Brain}\left(t\right)$ is the interplay coefficient of LF of heart to brain. A positive value of $SD{G}_{LF\to Brain}\left(t\right)$ indicates that brain increases linearly with increasing LF of heart, and vice versa. The level of $SD{G}_{LF\to Brain}\left(t\right)$ indicates the strength of LF of heart interplay of brain. $SD{G}_{HF\to Brain}\left(t\right)$ is the interplay coefficient of heart HF to brain. The value of $SD{G}_{HF\to Brain}\left(t\right)$ is positive. It shows that with the increase of heart HF, brain increases linearly, and vice versa. The level of $SD{G}_{HF\to Brain}\left(t\right)$ indicates the strength of heart HF interplay of brain.

Interplay coefficients proposed by the SDG model were averaged for each subject (N = 26) at the time points, both AS and MGS belong to the voluntary saccade tasks, so we averaged the interplay coefficients of the two saccade tasks. To reduce the effect of individual differences, the median was calculated between all participants for each EEG channel, and Friedman’s test, followed by Bonferroni correction ($\alpha$ = 0.05/3), was used to assess the differences in the interplay coefficients across different frequency bands, while the Wilcoxon test ($\alpha$ = 0.05), was applied for between-group comparisons. No correction method was applied for the between-group comparisons because the group comparisons in this study did not involve multiple comparisons. Both tests were conducted using a two-tailed approach and were suitable for data that did not follow a normal distribution. For the analysis of the saccade task, we did not perform a trial-by-trial analysis because the time of saccade is very fast, and one trial is not enough to calculate the brain-heart interplay coefficient. The resultant plots were all plotted by MATLAB.

The interplay coefficients of brain to LF of heart under the two groups saccade tasks were shown in Fig. 3A. Only gamma band showed significant interplay coefficients in the two groups of saccade tasks, especially in the parietal region, whereas the interplay coefficients in other frequency bands (below alpha and beta) were not significant (Table 1). The values of the interplay coefficients from the brain to the LF of the heart were all negative. The negative interplay coefficients suggest that during the performance of the saccade tasks, the increase of brain activity led to a decrease in the LF of heart. This may stem from the cognitive load triggered by the demands of the task [42]. We observed that the absolute mean interplay coefficients across different frequency bands were higher for the voluntary saccade task compared to the reflexive saccade task. This difference may be attributed to the voluntary saccade task engaging higher cortical regions, thereby increasing cognitive load and resulting in stronger negative interplay relationships.

The pattern of brain to HF of the heart interplay was shown in Fig. 3B. As can be seen, the below alpha band of EEG demonstrated limited interplay with HF of the heart. However, the interplay coefficient rises as the frequency band increases (Table 1). The most pronounced interplay influence was observed in the gamma band within the parietal areas. All frequency bands showed positive interplay coefficients. The positive interplay coefficients suggested that the increased brain leads to increased heart HF during the performance of a saccade task [43]. No distinct pattern was observed in the brain average interplay coefficients across different frequency bands between the two groups of saccade tasks. This requirement may have driven the brain to regulate heart HF in a consistent manner across all tasks, thereby effectively mitigating excessive physiological stress [44].

Average interplay coefficients from the different frequency bands of the brain to the LF and HF of the heart proved that the interplay coefficients of the gamma band were the strongest. To further prove this, we performed a p-value test on the interplay coefficients in different frequency bands for two groups of saccade tasks (shown in Fig. 4). p-values were obtained by Friedman’s test and Bonferroni correction. As can be seen from the figure, there are significant differences in the gamma band.

Subsequently, comparison of the interplay coefficients from brain to heart among two groups of saccade tasks was conducted, significant difference was found in the interplay coefficients from brain to LF of heart (shown in Fig. 5). And there were no significant differences in the interplay coefficients from brain to HF of heart. The main significant differences were found in the prefrontal lobes that closely related to cognition and memory (p = 0.0163).

The interplay patterns from LF and HF of the heart to the brain under the two groups saccade tasks were shown in Fig. 6A,B. As can be seen, the interplay coefficients from LF and HF of the heart to EEG different bands were the most significant below alpha band under the two groups saccade tasks (Table 2). The average brain interplay coefficients were all positive and the values of average brain interplay coefficients of reflexive saccade task was greater than voluntary saccade task.

The positive interplay coefficients suggest that heart is positively correlated to brain [43]. And both the interplay coefficients from LF and HF of the heart to EEG different bands exhibit reflexive saccade task is greater than voluntary saccade task. This may be due to the limited attentional and cognitive resources required for the brain in the reflexive saccade task, and the brain has more capacity to contact the heart. In contrast, voluntary saccade tasks require a higher degree of attention and cognitive resources, and the higher cortex of the brain is more involved, requiring more resources to be devoted to processing information and performing tasks. Although the heart is still in contact with brain, the interplay from heart to brain may be affected to some extent with the increase in cognitive resources and the complexity of the task [45].

Similarly, to verify that the interplay coefficients from the LF and HF of the heart to the brain are most significant below alpha band, we also did a p-value test for the different frequency bands under the two groups of saccade task (Fig. 7). p-values were obtained by Friedman’s test and Bonferroni correction. As can be seen in the figure, there are significant differences below alpha band.

Comparison of the interplay coefficients from heart to brain among different groups of saccade tasks was conducted, significant difference was found in the interplay coefficients from LF and HF of heart to brain (shown in Fig. 8). Significant differences were observed in the interplay coefficients from both the heart LF and HF to the brain below alpha band in the parietal (p = 0.0203) and occipital regions (p = 0.0097).