Autonomic nervous activity analysis based on visibility graph complex networks and skin sympathetic nerve activity

Background: Autonomic nerve system (ANS) plays an important role in regulating cardiovascular function and cerebrovascular function. Traditional heart rate variation (HRV) and emerging skin sympathetic nerve activity (SKNA) analyses from ultra-short-time (UST) data cannot fully reveal neural activity, thereby quantitatively reflect ANS intensity.

Methods: Electrocardiogram and SKNA from sixteen patients (seven cerebral hemorrhage (CH) patients and nine control group (CO) patients) were recorded using a portable device. Ten derived HRV (mean, standard deviation and root mean square difference of sinus RR intervals (NNmean, SDNN and RMSSD), ultra-low frequency (<0.003 Hz, uLF), very low frequency ([0.003 Hz, 0.04 Hz), vLF), low frequency ([0.04 Hz, 0.15 Hz), LF) and high frequency power ([0.15 Hz, 0.4 Hz), HF), ratio of LF to HF (LF/HF), the standard deviation of instantaneous beat-to-beat R-R interval variability (SD1), and approximate entropy (ApEn)) and ten visibility graph (VG) features (diameter (Dia), average node degree (aND), average shortest-path length (aSPL), clustering coefficient (CC), average closeness centrality (aCC), transitivity (Trans), average degree centrality (aDC), link density (LD), sMetric (sM) and graph energy (GE) of the constructed complex network) were compared on 5-min and UST segments to verify their validity and robustness in discriminating CH and CO under different data lengths. Besides, their potential for quantifying ANS-Load were also investigated.

Results: The validation results of HRV and VG features in discriminating CH from CO showed that VG features were more clearly distinguishable between the two groups than HRV features. For effectiveness evaluation of analyzing ANS on UST segment, the NNmean, SDNN, RMSSD, LF, HF and LF/HF in HRV features and the CC, Trans, Dia and GE of VG features remained stable in both activated and inactivated segments across all data lengths. The capability of HRV and VG features in quantifying ANS-Load were evaluated and compared under different ANS-Load, the results showed that most HRV features (SDNN, LFHF, RMSSD, vLF, LF and HF) and almost all VG features were correlated to sympathetic nerve activity intensity.

Conclusions: The proposed autonomic nervous activity analysis method based on VG and SKNA offers a new insight into ANS assessment in UST segments and ANS-Load quantification.

1 Introduction

The autonomic nerve system (ANS), composed of the sympathetic nervous system and the parasympathetic nervous system, plays an important role in regulating cardiovascular function and cerebrovascular function. Dysregulation of the ANS can affect the brain’s perception of various stressors, disrupt the adaptive capacity of homeostasis restoration, and ultimately increase the risk of stress-related disorders such as cardiac arrhythmia, hypertension, atherosclerosis, and stroke (Sternberg and Schaller, 2020). More recently, ANS modulation has been proposed as a promising therapeutic strategy for the management of autonomic dysfunction-related stroke (Mo et al., 2019). Therefore, ANS monitoring and analysis before the development of stress-related disorders is of paramount importance for improving the prognosis of patients with dysfunction-related stroke.

The most widely used clinical assessment method of sympathetic nerve activity (SNA) is evaluating end-organ responses to ANS physiological stimuli, such as tilt table testing, valsalva maneuver, plasma catecholamines, baroreflex sensitivity, thermoregulatory sweat test, and heart rate variation (HRV) (Thomas et al., 2019). Among these tests, HRV analysis is a widely accepted and implemented method to non-invasively and conveniently assess sympatho-vagal balance (Thomas et al., 2019). In general, HRV is generated and analyzed from the long-term electrocardiogram (ECG) waves, such as 24-h Holter, and its change can reflect the dynamic/trend of ANS activity over time (Bodapati et al., 2017; Schneider et al., 2018). It is reported that poststroke patients with raised SNA and low HRV are at higher risk for arrhythmias (atrial fibrillation, ventricular tachyarrhythmia) or other ECG changes (prolonged QT, inversed T wave) (Constantinescu et al., 2018). In addition, HRV is also used as a biomarker for classifying acquired brain injury patients and healthy controls (Galea et al., 2018). Meanwhile, multiple functional outcomes (cognitive functions, physical activity, and emotional expression) can be manifested in HRV (Forte et al., 2019; Sharma et al., 2019). Thus, HRV can not only serve as an indicator of cardiac function, but also reflect the central modulation capacity to stress (Yperzeele et al., 2015; Fyfe-Johnson et al., 2016; Kim et al., 2018). However, HRV quantifies ANS modulation at the sinoatrial level, which is difficult to generalize to cardiac patients with abnormal rhythms (atrial fibrillation, premature beats, etc.) (Zhao et al., 2020).

As a non-invasive and versatile SNA assessment method, skin sympathetic nerve activity (SKNA) has been applied to many clinical events (Doytchinova et al., 2017; Kusayama et al., 2020) and been proven to have the potential to predict sympathetic tone in many applications (i.e., acute myocardial infarction (He et al., 2020), neurologic recovery patients (Liu et al., 2021a), and sleep apnea (Kutkut et al., 2021)). To this juncture, several parameters have been derived from SKNA to quantify SNA. The average voltage of SKNA (aSKNA) is validated to be correlated with heart rate, and can be used as a biomarker for fitness level and efficacy of exercise training (Liu et al., 2021b). The burst numbers of SKNA (bSKNA) and variable value of SKNA (vSKNA) (Zhang et al., 2019) are higher in ventricular arrhythmia patients than in control groups, indicating SKNA can be used to predict the ventricular arrhythmogenesis recurrence. The envelope of SKNA (eSKNA) is extracted to depict the temporal pattern of SKNA, and the cross-comparison results between SKNA clustering groups and non-SKNA clustering groups demonstrate that eSKNA can act as a valid surrogate marker to classify ANS regulation ability in acute myocardial infarction patients (Liu et al., 2021a). Although these parameters can reflect the ANS changes by empirical threshold-based nerve bursts detection, the low signal-to-noise ratio of SKNA will lead to misjudgments (Xing et al., 2022a). In addition, the low amplitude SKNA signal (0.5–80 µV) is susceptible to noise, increasing the difficulty of extracting sympathetic-related information (Zhang et al., 2022). Therefore, more work is still needed to effectively analyze the autonomic nervous activity from SKNA, especially in real-time application.

Complex network is an emerging nonlinear dynamics analysis method for complex systems. It has been employed in a variety of physical and engineering systems: weather conditions (Fang et al., 2017), finance (Zhao et al., 2018), biomedical applications (Gao et al., 2021). Recently, several network-based approaches have been proposed to map time series into complex networks, such as visibility graph (Xu et al., 2018), recurrence plot (Eroglu et al., 2018), ordinal partition network (Santos et al., 2022). In particular, visibility graph (VG) is a simple and fast computational framework for us to bridge the gap between time series and complex networks, and it has been successfully implemented in different fields. Bhaduri and Ghosh (2016) studied cardiac dynamics during meditation through multi-fractal detrended fluctuation and RR interval-based VG, and they found that VG was superior to multi-fractal detrended fluctuation in reflecting physiological effects on subjects undergoing meditation. Gao et al. (2017) developed a time-dependent limited penetrable VG, and applied it to RR intervals for classifying heart states of healthy, congestive heart failure and atrial fibrillation. León et al. (2020) used HRV features and VG features derived from the heart rate time series for the prediction of late onset sepsis in preterm infants, the results showed that the VG features in HRV analysis were useful for sepsis prediction in newborns. From these studies, VG complex networks are often constructed from RR intervals for heart rate-related applications, while no work has focused on the application of VG analysis in evaluating ANS with SKNA.

In this study, an autonomic nervous activity analysis method was proposed based on VG complex network and SKNA. Based on previous studies (Naredi et al., 2000; Chun-jing et al., 2013), we hypothesized that SNA was elevated in patients with cerebral hemorrhage (CH). Therefore, we collected ECG and SKNA from CH patients and control group (CO), and compared the derived HRV and VG features to evaluate their effectiveness in distinguishing CH from CO. In addition, the ANS analysis performance of HRV and VG features on ultra-short-time (UST) segments were evaluated to verify their robustness under different data lengths. Finally, the correlations between HRV and VG features and ANS-Load were investigated under different data lengths to explore their potential for quantifying the intensity of SNA.

2 Methods

2.1 Data acquisition

The ECG and SKNA were recorded by a portable data acquisition device designed in our previous work (Xing et al., 2022b). It consists of low-noise analog-front-end (ADS1299, Texas Instruments, Dallas, TX) for bio-potential signal acquisition, a microcontroller (STM32L476, STMicroelectronics) for the management of the whole system, and a power management circuits (powered by a 3.7 V rechargeable lithium polymer battery). In order to reduce the system noise floor, a low-noise first-stage amplifier (INA128) was implemented with the ADS1299 chip. The clinical signals were measured at 4 kHz sampling frequency using conventional disposable silver/silver-chloride (Ag/AgCl) electrodes attached to the users’ chest. The signal measurements were carried out in a noise-free sound insulation room. After an adjustment period of at least 10-min, the 10-min signal of each subject was acquired in a supine position. The recorded signals were stored on a local trans-flash card, and processed off-line by MATLAB.

2.2 Patients

Patients with spontaneous CH who had a history of hypertension were recruited. All patients were male and had no definite cardiovascular and cerebrovascular events other than hypertension. The location of cerebral hemorrhage in all patients was located in the basal ganglia, and the hemorrhage did not break into the ventricle. The course of cerebral hemorrhage had passed through the acute phase and was in the subacute phase. Age- and sex-matched normal volunteers, no other obvious cardiovascular and cerebrovascular diseases except hypertension, were recruited as CO from the hypertension clinic. All the patients were enrolled from the Department of Neurosurgery, First Affiliated Hospital of Nanjing Medical University from October 2021 to December 2021. Exclusion criteria included: 1) patients with traumatic cerebral hemorrhage, ischemic stroke or hemorrhagic conversion; 2) cerebral hemorrhage patients underwent the unstable phase (with shock, large fluctuations in heart rate or blood pressure); 3) patients with thyroid disease, diabetes, cardiac arrhythmia, and other disorders that may affect ANS.

Sixteen patients were enrolled in this study, including seven CH patients and nine CO patients. A 10-min single-lead ECG and SKNA were recorded in a supine position for each patient, and they were asked to avoid unnecessary movement during the recording. Three Ag/AgCl electrodes were placed in the left subclavian, right subclavian, and right lower abdomen, and the sampling rate was 4 kHz.

2.3 Data process

2.3.1 Signal pre-processing

Due to the small amplitudes of ECG and SKNA, the signal is easily contaminated by various noises. Therefore, the signal quality is visually assessed before signal processing. Those episodes that are corrupted by severe background noise and cannot distinguish QRS complexes are eliminated. Afterwards, only 5-min segments with more than 90% high signal quality are reserved, and the ECG and SKNA are extracted from these segments by 10th-order Butterworth bandpass filters with cutoff frequencies of 0.5–150 Hz and 500–1,000 Hz, respectively. For further HRV analysis, the QRS complexes are identified by P&T method (Pan and Tompkins, 1985), and false and missing detection are calibrated artificially. To clearly label neural clusters, eSKNA was extracted by performing moving average (MA) and root mean square (RMS) on SKNA (Eqs. 1, 2). Referring to (Liu et al., 2021a), the window size and sliding step of MA are 100-ms and 2-ms, respectively.

${\mathit{x}}_{\mathit{MAk}}\mathrm{=}\frac{{\displaystyle {\sum }_{i\mathrm{=}s\mathrm{\ast }k}^{s\mathrm{\ast }k\mathrm{+}w}{\mathit{x}}_{\mathit{i}}}}{\mathit{m}}(1)$

$\mathit{j}\mathrm{=}\frac{(Fs\mathrm{\ast }T\mathrm{-}w)\mathrm{\ast }s}{Fs}(2)$

where ${\mathit{X}}_{\mathit{MA}}$ is the array of input signal after MA; ${\mathit{x}}_{\mathit{MAk}}$ is the kth sample of ${\mathit{X}}_{\mathit{MA}}$; j is the number of ${\mathit{X}}_{\mathit{MA}}$; w is the window size; s is the sliding step; ${\mathit{x}}_{\mathit{i}}$ is the ith sample of the input signal; Fs is the sampling frequency; T is the duration in second of selected data.

In RMS calculator (Eq. 3), the ${\mathit{X}}_{\mathit{RMS}}$ is extracted from ${\mathit{X}}_{\mathit{MA}}$ with a window size of 100-ms and a sliding step of 2-ms:

${\mathit{X}}_{\mathit{RMS}}\mathrm{=}\sqrt{\frac{{\displaystyle {\sum }_{i\mathrm{=}k}^{k\mathrm{+}j}{\mathit{x}}_{\mathit{MAi}}^2}}{\mathit{n}}}(3)$

where n is the number of samples in a window; ${\mathit{x}}_{\mathit{RMS}}$ is the kth sample value of RMS; j is the number of samples of ${\mathit{X}}_{\mathit{RMS}}$; ${\mathit{X}}_{\mathit{MAi}}$ is the ith sample value of the array ${\mathit{X}}_{\mathit{RMS}}$.

The ${\mathit{X}}_{\mathit{RMS}}$ is defined as eSKNA, and a threshold-based method is performed on it for SKNA bursts determination. The threshold is calculated as follows:

$\mathit{Threshold}\mathrm{=}\left(Baseline\mathrm{-}Min\right)\mathrm{\ast }5\mathrm{+}Min(4)$

where Baseline is the average of the lower 20% samples in the selected window; Min is the minimum of the selected window.

In order to analyze the effectiveness of VG features in quantifying ANS from UST segments, the 5-min signals were split into 10-s, 20-s, 30-s, 40-s, 50-s and 60-s segments, respectively. The burst load of each segment was calculated as the ratio of burst time to total time, and partitioned to 5 equal intervals from 0 to 1 ([0, 0.2), [0.2, 0.4), [0.4, 0.6), [0.6, 0.8), [0.8, 1.0]). Then, the segments were marked as activated (burst load > 0) and inactivated (burst load = 0) according to the burst load. Thereafter, the HRV and VG analysis were conducted on these data. The flowchart of this paper is illustrated in Figure 1.

FIGURE 1

FIGURE 1. The flowchart of data process in this paper. MA and RMS mean the operation of moving average and root mean square, respectively.

2.3.2 Heart rate variation analysis

The time-domain features, frequency-domain features, and nonlinear features were extracted by the PhysioNet Cardiovascular Signal Toolbox (Vest et al., 2018). The time-domain analysis included sinus RR intervals-related features (mean (NNmean), standard deviation (SDNN), and root mean square difference (RMSSD)). The frequency-domain features consisted of the power in different frequency range (ultra-low frequency power (<0.003 Hz, uLF), very low frequency power ([0.003 Hz, 0.04 Hz), vLF), low frequency power ([0.04 Hz, 0.15 Hz), LF), high frequency power ([0.15 Hz, 0.4 Hz), HF), ratio of LF to HF (LFHF)). The nonlinear features consisted of the standard deviation of projection of the Poincaré Plot on the line perpendicular to the line of identity (y = −x, SD1), and approximate entropy (ApEn). The standard deviation of the projection of the Poincaré Plot on the line of identity (y = x, SD2) was not included in this study because it was not suitable for UST HRV analysis.

2.3.3 Visibility graph features extraction of envelope of SKNA

As a natural graph-theoretical description of nonlinear systems, VG can simply convert a time series into a scale-free network. The statistical measures from the constructed network can demonstrate the dynamic behaviors of the nonlinear systems, and have been proven to be related to the self-similarity and complexity of the time series (Bhaduri and Ghosh, 2016). Every data point of time series is mapped to a node in its associated VG network, and an edge between two nodes is connected if the corresponding time samples can view each other. Suppose the original time series is X, and the ith point of X is X_i. Two data points (t_a, X_a) and (t_c, X_c), at time t_a and t_c, are connected if and only if any other data (t_b, X_b) between them (t_a < t_b < t_c) satisfies the following criterion:

${\mathit{X}}_{\mathit{b}}\mathrm{<}{\mathit{X}}_{\mathit{c}}\mathrm{+}\left({\mathit{X}}_{\mathit{a}}\mathrm{-}{\mathit{X}}_{\mathit{c}}\right)\mathrm{\ast }\frac{t_c\mathrm{-}{t}_b}{t_c\mathrm{-}{t}_a}(5)$

The VG network are applied on eSKNA, and ten network measures (Pan and Tompkins, 1985; Hou et al., 2016; Vest et al., 2018; Xing et al., 2022b; Santos et al., 2022) are extracted for further analysis in this study:

2.3.3.1 Diameter

The longest shortest path between any two nodes in the network (Eq. 6).

$\mathit{Dia}\mathrm{=}\mathit{ma}{\mathit{x}}_{\mathit{i}\mathrm{,}\mathit{j}}{\mathit{D}}_{\mathit{ij}}(6)$

where D_ij is the length of the shortest path between node i and node j.

2.3.3.2 Average node degree

The degree of a node in a graph is defined as the number of connected edges to this node, and the mean degree (Eq. 7) is the average value of all node’s degree in this graph (León et al., 2020).

$\mathit{aND}\mathrm{=}\frac{\mathrm{1}}{\mathrm{{\rm N}}}{\displaystyle \sum _{n\mathrm{=}1}^N{d}_{\mathit{n}}}(7)$

where N is the total number of nodes, and d_n is the degree of node n.

2.3.3.3 Average shortest-path length

The shortest path is a reflection of transmission and communication in the graph, the average shortest path length (Eq. 8) is the average of shortest path over all couples of nodes (Hou et al., 2016).

$\mathit{aSPL}\mathrm{=}\frac{\mathrm{1}}{\mathrm{{\rm N}\ast }\left(\mathrm{{\rm N}-1}\right)}{\displaystyle {\sum }_{\mathit{i}\mathrm{,}\mathit{j}\in \mathit{V}\mathrm{,}\mathit{i}\ne \mathit{j}}{D}_{\mathit{ij}}}(8)$

where N is the total number of nodes, and V is the set of N nodes.

2.3.3.4 Clustering coefficient

The cluster coefficient of a node in a graph is the ratio of all triangles involving that node to the number of connected triples centered on that node, and the cluster coefficient of a graph (Eq. 9) is the average of the cluster coefficient of all nodes (León et al., 2020).

$\mathit{CC}\mathrm{=}\frac{\mathrm{1}}{\mathrm{{\rm N}}}\mathrm{\ast }{\displaystyle {\sum }_{\mathit{i}\in \mathit{V}}\frac{{\mathit{e}}_{\mathit{i}}}{{\mathit{k}}_{\mathit{i}}\mathrm{\ast }\left({\mathit{k}}_{\mathit{i}}\mathrm{-}1\right)}}(9)$

where e_i is the actual number of edges between all the couples of neighbors of node i, and k_i*(k_i-1) is the maximum possible number of edges between all the k_i neighbors of node i.

2.3.3.5 Average closeness centrality

Closeness centrality is the sum of the distances from a node to other nodes, representing the convenience and ease of connection between the focal node and other nodes (Zhang and Luo, 2017).

$\mathit{aCC}\mathrm{=}\frac{\mathrm{1}}{\mathrm{{\rm N}}}{\displaystyle \sum _{i\mathrm{=}1}^N\frac{N\mathrm{-}1}{{\displaystyle {\sum }_{j\mathrm{=}1}^N{D}_{ij}}}}(10)$

2.3.3.6 Transitivity

The transitivity (Eq. 11) is the ratio between the triangle numbers and the connected triple numbers in a graph to obtain the global information of CC (León et al., 2020).

$\mathit{Trans=}\frac{\mathit{3\ast number}\mathit{of}\mathit{trangles}\mathit{in}\mathit{the}\mathit{graph}}{\mathit{number}\mathit{of}\mathit{connected}\mathit{triples}\mathit{in}\mathit{the}\mathit{graph}}(11)$

2.3.3.7 Average degree centrality

Degree centrality is defined as the ratio between the number of nodes connected to the current node, and the total number of all nodes in the network (Zhang and Luo, 2017).

$\mathit{aDC}\mathrm{=}\frac{\mathrm{1}}{\mathrm{{\rm N}}}{\displaystyle \sum _{i\mathrm{=}1}^N\frac{{\displaystyle {\sum }_{j\mathrm{=}1}^N{e}_{ij}}}{N\mathrm{-}1}}(12)$

2.3.3.8 Link density

Link density (Eq. 13) is the ratio between the number of edges and the maximum possible number of edges (N*(N-1)/2) (Liu et al., 2015).

$\mathit{LD}\mathrm{=}\frac{{\displaystyle {\sum }_{i\mathrm{=}1}^{\mathrm{{\rm N}}}{\displaystyle {\sum }_{j\mathrm{=}1}^N{e}_{ij}}}}{N\ast \left(N\mathrm{-}1\right)/2}(13)$

2.3.3.9 sMetric

The sMetric (Eq. 14) is the sum of products of degrees across all edges (Li et al., 2005).

$sM\mathrm{=}{\displaystyle \sum _{i\mathrm{,}j\mathrm{=}1}^{\mathrm{{\rm N}}}{d}_i\ast d_j}(14)$

2.3.3.10 Graph energy

Graph Energy (Eq. 15) is defined as the sum of the absolute values of the real components of the eigenvalues (${\lambda }_i$) of the graph (Balakrishnan, 2004).

$GE\mathrm{=}{\displaystyle {\sum }_i\left|{\lambda }_i\right|}(15)$

3 Experiments and results

3.1 Comparison of heart rate variation and visibility graph on autonomic nerve system analysis

The eSKNA segments are converted into a scale-free graph by natural VG method. The typical 10-s eSKNA of CH and CO segments are illustrated in Figure 2A, and their corresponding VG complex networks are shown in Figures 2B,C, respectively. For clear demonstration, the communities of these complex networks are colorized according to their modularity classes by Gephi software. It can be seen that the amplitudes of CH eSKNA fluctuate smoothly, while there is a clear burst in CO eSKNA around about 6-s time point (Figure 2A). The communities of the CH complex network are dispersed as all intermediate peaks obstruct the visible range between the front and rear peaks (Figure 2B). Conversely, the CO complex network consists of several small communities crowded with a central community, since the burst can view almost all other nodes in the network (Figure 2C).

FIGURE 2

FIGURE 2. The typical eSKNA signals and their corresponding VG complex networks of CH and CO segments. (A) The typical eSKNA and their VG, (B) and (C) are the complex networks of CH and CO with colored communities.

To evaluate the validity of HRV and VG features from the 5-min signal, these features are normalized to [0, 1], and the comparisons between CH and CO for these features are illustrated in Figure 3. To further quantify the distribution differences between CH and CO, the WRS test is carried out for each feature. Significant difference (p < 0.05) between two groups is marked with red “*”, and extremely significant difference (p < 0.01) is marked with red “**“. In this paper, p values less than 0.05 were regarded as statistically significant for each test. Statistical analyses were performed using MATLAB (R2022a) on a PC with Intel® Core™ i7-7700 3.6 GHz processor and 32 GB RAM.

FIGURE 3

FIGURE 3. HRV and VG features between CH and CO. “*” stands for significant difference between two groups (p < 0.05) and “**” stands for extremely significant difference (p < 0.01). ANS Analysis Performance of VG on UST Segments.

It can be seen from Figure 3, the distributions of most HRV features between CH and CO are overlapped, indicating the difficulty to distinguish CH from CO by these features. In contrast, the distribution differences between CH and CO are evident in most VG features, validating the effectiveness of VG features in classifying CH from CO. In addition, almost all p-values calculated from VG features are < 0.01 (except Trans and Dia), while only SDNN, vLF and ApEn show significant difference in HRV features.

The ANS function assessment by HRV analysis is typically performed on either 5-min ECG recordings or nominal 24-h recordings, which limits its application in dynamic conditions, such as dynamic sympathetic assessment in athletes. To investigate the ANS analysis performance of these parameters on UST signals, the 5-min signals are split to 10-s, 20-s, 30-s, 40-s, 50-s and 60-s segments, respectively. In addition, these segments are marked as ANS activated and inactivated depending to whether they contain ANS bursts labeled from eSKNA. Note that segments without valid QRS complexes are removed. The final number of each data length is illustrated in Table 1. The distribution of each feature among different data length is compared and quantified by KW test. Furthermore, the total runtime of VG feature extraction for each data length is compared.

TABLE 1

TABLE 1. The number of data segments according to different data length, ANS status and burst load.

Figure 4 depicts the distribution differences of each feature (HRV and VG) for activated and inactivated segments under different data lengths. For almost all HRV features (except ApEn), the distribution for activated segments seldom changes with the data length increase, and the distribution for inactivated segments varies sightly in uLF, vLF and SD1. The ApEn for both activated and inactivated segments increases with data length expands. For VG features, the distributions for both activated and inactivated segments remain stable in CC, Trans, Dia and GE, and change slightly in aSPL and aCC. However, they decrease (increase) sharply with the data length increases in aND and aDC (sM).

FIGURE 4

FIGURE 4. The distribution of features (HRV and VG) for different data length. Effectiveness of HRV and VG on ANS-Load Determination.

To quantitatively characterize the stability of features (HRV and VG) in UST segments, their distribution differences (for activated and inactivated segments, respectively) in different data lengths (10-s, 20-s, 30-s, 40-s, 50-s and 60-s) are compared by KW test. The results are shown in Table 2. It is clear that the p-values of NNmean, SDNN, RMSSD, LF, HF, LFHF and SD1 are > 0.05 in both activated and inactivated segments, indicating that these features are not distributed differently across segments with different length. Similarly, the VG features of CC, Trans, Dia and GE also show no distribution differences in both activated and inactivated segments, as they all had p-values > 0.05. However, the p-values for vLF and ApEn of HRV features and aND, aDC, LD and sM of VG features are all < 0.01 for both groups, implying that their distribution varies significantly during data length increase.

TABLE 2

TABLE 2. The KW test results of feature distribution differences in different data lengths.

Although several features (HRV and VG) show stable performance in short-term segments, their efficiency still needs to be investigated to ensure their practical application. As we all known that the computational complexity of HRV features is very low, therefore, we only compare the running time of VG features under different data lengths. Figure 5 shows the histogram based on empirical cumulative distribution function and kernel density estimation of the running time for VG features extraction from different data lengths. The distribution is heavy-tailed in 10-s and 30-s segments, but appears approximately normal distribution in the remaining segments. The average time for each data length is around 15-s, 80-s, 210-s, 360-s, 565-s and 1,420-s, respectively. Obviously, the average runtime increases rapidly with data length expands and shows an exponential growth trend. The reason is that the nodes of the VG complex network increase with data length, resulting in a rapid growth of computational complexity for extracting features from the constructed adjacent matrix.

FIGURE 5

FIGURE 5. The running times of VG feature extraction under different data lengths.

In order to evaluate the capability of HRV and VG features in quantifying ANS-Load, the burst of each segment is determined by a threshold-based method. Furthermore, the burst load (ANS-Load) is extracted by calculating the ratio of the burst time to total time of the segment. The burst load is coarse-grained to 5 equal partitions from 0 to 1, the segment number of each partition under different data lengths is listed in Table 1. Then, the correlations between features and segment length under different burst load are quantified by Kendall rank correlation coefficient. As there are not enough ANS-Load in partitions [0.6, 0.8) and [0.8, 1.0], we only compare the correlation coefficient in ANS-Load among (0, 0.6).

The distributions between each HRV and VG feature and data lengths under different ANS-Load, associated with their mean values, are shown in Figure 6. In HRV features, the SDNN and LFHF (RMSSD, vLF, LF and HF) decrease (increase) with the ANS-Load increase in different data lengths, implying that the variation of ANS-Load would influence the time-domain and frequency-domain features of HRV. Meanwhile, almost all the VG features present an increasing or decreasing trend with the ANS-Load increase. The reason may be that the increased autonomic activity is reflected in increased bursts in eSKNA, resulting in the variation of connections between two nodes.

FIGURE 6

FIGURE 6. The distributions between each HRV and VG features and data lengths under different ANS-Load and their corresponding mean values.

The Kendall rank correlation coefficients for HRV and VG features are illustrated in Figure 7, and the red “*” and “**” represent p-values < 0.05 and < 0.01, respectively. In this study, we only focus on the degree of correlation other than its direction, therefore the positive correlation and negative correlation share the same color in Figure 7. It is obvious that there is a weak correlation between HRV features and ANS-Load, most correlation coefficients are around 0, and the maximum is 0.273 for uLF in 60-s segment. On the contrary, the VG features show a stronger correlation with ANS-Load, especially the correlation coefficient of aND reaches 0.526 in 60-s. In addition, the absolute values of the correlation coefficients are all above 0.13. Besides, only few HRV features show significant correlation between features and ANS-Load (i.e., NNmean in 10-s, LFHF in 20-s, vLF in 40-s). However, almost all the correlations between VG features and ANS-Load are extremely significant, indicating the potential of VG features for ANS-Load quantification.

FIGURE 7

FIGURE 7. The Kendall rank correlation coefficients between features and segment length under different burst load. “*” stands for significant correlation between features and ANS-Load (p < 0.05) and “**” stands for extremely significant correlation (p < 0.01).

4 Discussion

A VG and SKNA based autonomic nervous activity analysis method was proposed in this paper. SKNA overcomes the sinoatrial level limitations of traditional HRV analysis, and a number of SKNA-derived metrics had been proposed for SNA quantification. However, the noise-susceptibility of these metrics required new methods for quantitatively ANS assessment, especially for short-time segment application. As a nonlinear analysis method, VG offered a new insight into ANS assessment in short-term segments and ANS-Load quantification.

The unique contribution of this paper was the first application of VG on eSKNA for ANS assessment. HRV was the most widely used ANS assessment method, and its indices from time-domain, frequency-domain and nonlinear-domain had been proved to be biomarkers for cardiac arrhythmias, brain injury and emotion (Fyfe-Johnson et al., 2016; Constantinescu et al., 2018; Galea et al., 2018). However, HRV from UST data could not fully show the nerve activity. The newly proposed noninvasive cardiac SNA assessment method (SKNA) and its derived metrices has been widely used in clinical events, such as acute myocardial infarction (He et al., 2020), neurologic recovery patients (Liu et al., 2021a), and sleep apnea (Kutkut et al., 2021). Nevertheless, these SKNA-derived metrics could only qualitatively analyze ANS and could not quantitatively reflect ANS intensity or ANS-Load. Thanks to the nonlinear dynamics analysis method–VG complex networks, we could evaluate the ANS from network aspect. Although many previous studies had investigated and compared VG and HRV in meditation analysis (Bhaduri and Ghosh, 2016), sleep assessment (Hou et al., 2016) and congestive heart failure (Gao et al., 2017), this paper was the first to employ VG on eSKNA for ANS assessment. The comparison of HRV and VG features between CH and CO (Figure 2) showed that VG features are superior to HRV features in the ANS analysis. There were no significant distribution differences between CH and CO for most HRV features, while almost all VG features were clearly distinguishable between the two groups.

The stability of HRV and VG features in UST segments were compared by quantifying their distribution differences against different data lengths. Most HRV features remained stable for both activated and inactivated segments under all data lengths. In addition, the time-domain features (NNmean, SDNN, RMSSD) and most frequency-domain features (LF, HF, LFHF) manifested conformity in these segments, indicating that most UST HRV features could be used as surrogates for short-term HRV features. These results were consistent with Castald’s (Castaldo et al., 2019) investigation that NNmean and HF displayed consistency across all of the excerpt lengths (30 s, 1 min, 2 min, 3 min, and 5 min) for mental stress assessment. However, Jin Woong et al. (Kim et al., 2021) studied UST HRV in non-static conditions by comparing UST HRV features (10, 30, 60, 120, 180, and 240-s) with those from 5-min HRV, the results showed that UST HRV variables derived from the static condition could not applied to the non-static conditions of daily life. Similarly, the CC, Trans, Dia and GE of VG features remained stable in both activated and inactivated segments across all data lengths, implying that these VG features could reveal the dynamical changes caused by the adjustment of autonomous neural system from UST segments. Likewise, Jiang et al. (2013) applied VG to heartbeat interval time series for meditation investigation, and they also tested the stability of VG features on different length data, the results showed that the data length had no prominent effect on the VG analysis. The reason may be that the degree distribution persisted the same form for different length of data in any activated and inactivated segments.

HRV had been used as a biomarker for SNA measurement, but seldom been used for quantifying ANS-Load. The correlations between HRV features and ANS-Load were studied at different data lengths, and the comparison results showed that time-domain and frequency-domain features (SDNN, LFHF, RMSSD, vLF, LF and HF) had the potential to quantify ANS in UST segments. From a multimodal perspective, Debnath et al. (2021) designed a template matching algorithm to calculate scaled and stretched HRV features, associated with other features, for sympathetic and parasympathetic parameters determination. However, the acquisition of these employed features was complicated, and it still required more other biomarkers or calculated features to improve the quantification accuracy for practical clinical applications. The SKNA had been applied to evaluate the ANS as a non-invasive method in many clinical applications (Zhang et al., 2022), and many SKNA-derived metrices (e.g., aSKNA, bSKNA) had been validated and used for ANS qualification. Nevertheless, these parameters were susceptible to noise. In this paper, the VG features on eSKNA were extracted and compared across different ANS-Load, the results showed that almost all the VG features were correlated to ANS-Load. The link-related features (CC, aPL, Dia and LD) increased as ANS-Load grow, while the degree-related features (aCC, Trans, aND, aDC and sM) presented a decreasing trend across increasing ANS-Load. The increase in SNA intensity was manifested as the rising number and duration of bursts in eSKNA, which leads to a growth in the possibility of links between any two nodes in the VG. However, these links only concentrated on certain nodes (peak points of bursts), it meant that the node degrees of the entire VG would be aggregated into these nodes, resulting the increase of community numbers and the decrease of average degree. In addition, the total number of links grow exponentially with the total number of nodes in the network, while the degree distribution did not change with the node numbers (Tessone et al., 2011).

One limitation of our study is the small number of participants, further studies with larger cohorts are needed to confirm and strengthen these results. Another limitation is the VG features are only compared with HRV features, its validity still needs comparison with demographic information and laboratory tests for practical clinical applications. In addition, the robustness against noise of this method needs more efforts.

5 Conclusion

In summary, a VG on eSKNA based autonomic nervous activity analysis method was proposed in this paper. The comparison results of the HRV and VG features on CH and CO segments showed the superiority of VG features in ANS analysis. Furthermore, the ANS analysis performance of VG features on eSKNA signals with different data lengths demonstrated the stability of VG features (aND, aDC, LD and sM) in discriminating activated and inactivated segments at different data lengths. In addition, the capability of HRV and VG features to quantify SNA intensity was also evaluated, and the results showed that VG features had the potential to determine ANS-Load.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by The Bioethics Committee of the First Affliated Hospital of Nanjing Medical University. The patients/participants provided their written informed consent to participate in this study.

Author contributions

ZC conducted this study and wrote the manuscript. HC collected and analyzed the data. YX produced the devices and processed the signals. FC performed the statistical analysis. YZ guided the signal acquisition and quality control. ZC and CC contributed to the guidance of the whole study. All the authors made contributions to this work and approved the submitted manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 62001105.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Reference

Balakrishnan R. (2004). The Energy of a Graph. Linear Algebra Appl. 387, 287–295. doi:10.1016/j.laa.2004.02.038

CrossRef Full Text | Google Scholar

Bhaduri A., Ghosh D. (2016). Quantitative Assessment of Heart Rate Dynamics during Meditation: An ECG Based Study with Multi-Fractality and Visibility Graph. Front. Physiol. 7, 44. doi:10.3389/fphys.2016.00044

PubMed Abstract | CrossRef Full Text | Google Scholar

Bodapati R. K., Kizer J. R., Kop W. J., Kamel H., Stein P. K. (2017). Addition of 24‐hour Heart Rate Variability Parameters to the Cardiovascular Health Study Stroke Risk Score and Prediction of Incident Stroke: The Cardiovascular Health Study. J. Am. Heart Assoc. 6, e004305. doi:10.1161/JAHA.116.004305

PubMed Abstract | CrossRef Full Text | Google Scholar

Castaldo R., Montesinos L., Melillo P., James C., Pecchia L. (2019). Ultra-short Term HRV Features as Surrogates of Short Term HRV: A Case Study on Mental Stress Detection in Real Life. BMC Med. Inf. Decis. Mak. 19, 12–13. doi:10.1186/s12911-019-0742-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chun-jing H., Shan O., Guo-dong L., Hao-xiong N., Yi-ran L., Ya-ping F. (2013). Effect of Cervical Sympathetic Block on Cerebral Vasospasm after Subarachnoid Hemorrhage in Rabbits. Acta Cir. Bras. 28, 89–93. doi:10.1590/S0102-86502013000200001

PubMed Abstract | CrossRef Full Text | Google Scholar

Constantinescu V., Matei D., Costache V., Cuciureanu D., Arsenescu-Georgescu C. (2018). Linear and Nonlinear Parameters of Heart Rate Variability in Ischemic Stroke Patients. Neurol. Neurochir. Pol. 52, 194–206. doi:10.1016/j.pjnns.2017.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Debnath S., Levy T. J., Bellehsen M., Schwartz R. M., Barnaby D. P., Zanos S., et al. (2021). A Method to Quantify Autonomic Nervous System Function in Healthy, Able-Bodied Individuals. Bioelectron. Med. 7, 13–17. doi:10.1186/s42234-021-00075-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Doytchinova A., Hassel J. L., Yuan Y., Lin H., Yin D., Adams D., et al. (2017). Simultaneous Noninvasive Recording of Skin Sympathetic Nerve Activity and Electrocardiogram. Heart rhythm. 14, 25–33. doi:10.1016/j.hrthm.2016.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Eroglu D., Marwan N., Stebich M., Kurths J. (2018). Multiplex Recurrence Networks. Phys. Rev. E 97, 012312. doi:10.1103/PhysRevE.97.012312

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang K., Sivakumar B., Woldemeskel F. M. (2017). Complex Networks, Community Structure, and Catchment Classification in a Large-Scale River Basin. J. Hydrol. X. 545, 478–493. doi:10.1016/j.jhydrol.2016.11.056

CrossRef Full Text | Google Scholar

Forte G., Favieri F., Casagrande M. (2019). Heart Rate Variability and Cognitive Function: A Systematic Review. Front. Neurosci. 13, 710. doi:10.3389/fnins.2019.00710

PubMed Abstract | CrossRef Full Text | Google Scholar

Fyfe-Johnson A. L., Muller C. J., Alonso A., Folsom A. R., Gottesman R. F., Rosamond W. D., et al. (2016). Heart Rate Variability and Incident Stroke: the Atherosclerosis Risk in Communities Study. Stroke 47, 1452–1458. doi:10.1161/STROKEAHA.116.012662

PubMed Abstract | CrossRef Full Text | Google Scholar

Galea O. A., Cottrell M. A., Treleaven J. M., O Leary S. P. (2018). Sensorimotor and Physiological Indicators of Impairment in Mild Traumatic Brain Injury: a Meta-Analysis. Neurorehabil. Neural Repair 32, 115–128. doi:10.1177/1545968318760728

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao Z., Cai Q., Yang Y., Dang W. (2017). Time-dependent Limited Penetrable Visibility Graph Analysis of Nonstationary Time Series. Phys. A Stat. Mech. its Appl. 476, 43–48. doi:10.1016/j.physa.2017.02.038

CrossRef Full Text | Google Scholar

Gao Z., Dang W., Wang X., Hong X., Hou L., Ma K., et al. (2021). Complex Networks and Deep Learning for EEG Signal Analysis. Cogn. Neurodyn. 15, 369–388. doi:10.1007/s11571-020-09626-1

PubMed Abstract | CrossRef Full Text | Google Scholar

He W., Tang Y., Meng G., Wang D., Wong J., Mitscher G. A., et al. (2020). Skin Sympathetic Nerve Activity in Patients with Obstructive Sleep Apnea. Heart rhythm. 17, 1936–1943. doi:10.1016/j.hrthm.2020.06.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou F. Z., Li F. W., Wang J., Yan F. R. (2016). Visibility Graph Analysis of Very Short-Term Heart Rate Variability during Sleep. Phys. A Stat. Mech. its Appl. 458, 140–145. doi:10.1016/j.physa.2016.03.086

CrossRef Full Text | Google Scholar

Jiang S., Bian C., Ning X., Ma Q. D. (2013). Visibility Graph Analysis on Heartbeat Dynamics of Meditation Training. Appl. Phys. Lett. 102, 253702. doi:10.1063/1.4812645

CrossRef Full Text | Google Scholar

Kim H., Cheon E., Bai D., Lee Y. H., Koo B. (2018). Stress and Heart Rate Variability: A Meta-Analysis and Review of the Literature. Psychiatry Investig. 15, 235–245. doi:10.30773/pi.2017.08.17

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim J. W., Seok H. S., Shin H. (2021). Is Ultra-short-term Heart Rate Variability Valid in Non-static Conditions? Front. Physiol. 12, 596060. doi:10.3389/fphys.2021.596060

PubMed Abstract | CrossRef Full Text | Google Scholar

Kusayama T., Wong J., Liu X., He W., Doytchinova A., Robinson E. A., et al. (2020). Simultaneous Noninvasive Recording of Electrocardiogram and Skin Sympathetic Nerve Activity (neuECG). Nat. Protoc. 15, 1853–1877. doi:10.1038/s41596-020-0316-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutkut I., Uceda D., Kumar A., Wong J., Li X., Wright K. C., et al. (2021). Skin Sympathetic Nerve Activity as a Biomarker for Neurologic Recovery during Therapeutic Hypothermia for Cardiac Arrest. Heart rhythm. 18, 1162–1170. doi:10.1016/j.hrthm.2021.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

León C., Carrault G., Pladys P., Beuchée A. (2020). Early Detection of Late Onset Sepsis in Premature Infants Using Visibility Graph Analysis of Heart Rate Variability. IEEE J. Biomed. Health Inf. 25, 1006–1017. doi:10.1109/JBHI.2020.3021662

CrossRef Full Text | Google Scholar

Li L., Alderson D., Doyle J. C., Willinger W. (2005). Towards a Theory of Scale-free Graphs: Definition, Properties, and Implications. Internet Math. 2, 431–523. doi:10.1080/15427951.2005.10129111

CrossRef Full Text | Google Scholar

Liu C., Lee C., Lin S., Tsai W. (2021). Temporal Clustering of Skin Sympathetic Nerve Activity Bursts in Acute Myocardial Infarction Patients. Front. Neurosci. 15, 720827. doi:10.3389/fnins.2021.720827

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu X., Kumar A., O'Neil J., Wong J., Saadoon O., Kadire S., et al. (2021). Skin Sympathetic Nerve Activity as a Biomarker of Fitness. Heart rhythm. 18, 2169–2176. doi:10.1016/j.hrthm.2021.08.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu Z., Dong W., Fu Y. (2015). Local Degree Blocking Model for Link Prediction in Complex Networks. Chaos 25, 013115. doi:10.1063/1.4906371

PubMed Abstract | CrossRef Full Text | Google Scholar

Mo J., Huang L., Peng J., Ocak U., Zhang J., Zhang J. H. (2019). Autonomic Disturbances in Acute Cerebrovascular Disease. Neurosci. Bull. 35, 133–144. doi:10.1007/s12264-018-0299-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Naredi S., Lambert G., Eden E., Zall S., Runnerstam M., Rydenhag B., et al. (2000). Increased Sympathetic Nervous Activity in Patients with Nontraumatic Subarachnoid Hemorrhage. Stroke 31, 901–906. doi:10.1161/01.STR.31.4.901

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan J., Tompkins W. J. (1985). A Real-Time QRS Detection Algorithm. IEEE Trans. Biomed. Eng. 32, 230–236. doi:10.1109/TBME.1985.325532

PubMed Abstract | CrossRef Full Text | Google Scholar

Santos L. D., Corrêa D. C., Walker D. M., de Godoy M. F., Macau E. E., et al. (2022). Characterisation of Neonatal Cardiac Dynamics Using Ordinal Partition Network. Med. Biol. Eng. Comput. 60, 829–842. doi:10.1007/s11517-021-02481-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider U., Bode F., Schmidt A., Nowack S., Rudolph A., Doelcker E., et al. (2018). Developmental Milestones of the Autonomic Nervous System Revealed via Longitudinal Monitoring of Fetal Heart Rate Variability. PLoS One 13, e0200799. doi:10.1371/journal.pone.0200799

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma R. R., Kumar A., Pachori R. B., Acharya U. R. (2019). Accurate Automated Detection of Congestive Heart Failure Using Eigenvalue Decomposition Based Features Extracted from HRV Signals. Biocybern. Biomed. Eng. 39, 312–327. doi:10.1016/j.bbe.2018.10.001

CrossRef Full Text | Google Scholar

Sternberg Z., Schaller B. (2020). Central Noradrenergic Agonists in the Treatment of Ischemic Stroke—An Overview. Transl. Stroke Res. 11, 165–184. doi:10.1007/s12975-019-00718-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Tessone C. J., Geipel M. M., Schweitzer F. (2011). Sustainable Growth in Complex Networks. EPL Europhys. Lett. 96, 58005. doi:10.1209/0295-5075/96/58005

CrossRef Full Text | Google Scholar

Thomas B. L., Claassen N., Becker P., Viljoen M. (2019). Validity of Commonly Used Heart Rate Variability Markers of Autonomic Nervous System Function. Neuropsychobiology 78, 14–26. doi:10.1159/000495519

PubMed Abstract | CrossRef Full Text | Google Scholar

Vest A. N., Da Poian G., Li Q., Liu C., Nemati S., Shah A. J., et al. (2018). An Open Source Benchmarked Toolbox for Cardiovascular Waveform and Interval Analysis. Physiol. Meas. 39, 105004. doi:10.1088/1361-6579/aae021

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing Y., Zhang Y., Xiao Z., Yang C., Li J., Cui C., et al. (2022). An Artifact-Resistant Feature SKNAER for Quantifying the Burst of Skin Sympathetic Nerve Activity Signal. Biosensors 12, 355. doi:10.3390/bios12050355

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing Y., Zhang Y., Yang C., Li J., Li Y., Cui C., et al. (2022). Design and Evaluation of an Autonomic Nerve Monitoring System Based on Skin Sympathetic Nerve Activity. Biomed. Signal Process. Control 76, 103681. doi:10.1016/j.bspc.2022.103681

CrossRef Full Text | Google Scholar

Xu P., Zhang R., Deng Y. (2018). A Novel Visibility Graph Transformation of Time Series into Weighted Networks. Chaos, Solit. Fractals 117, 201–208. doi:10.1016/j.chaos.2018.07.039

CrossRef Full Text | Google Scholar

Yperzeele L., Van Hooff R., Nagels G., De Smedt A., De Keyser J., Brouns R. (2015). Heart Rate Variability and Baroreceptor Sensitivity in Acute Stroke: a Systematic Review. Int. J. Stroke 10, 796–800. doi:10.1111/ijs.12573

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang J., Luo Y. (2017). “Degree Centrality, Betweenness Centrality, and Closeness Centrality in Social Network,” in Proceedings of the 2017 2nd International Conference on Modelling, Simulation and Applied Mathematics, 300–303.

CrossRef Full Text | Google Scholar

Zhang P., Liang J., Cai C., Tian Y., Dai M., Wong J., et al. (2019). Characterization of Skin Sympathetic Nerve Activity in Patients with Cardiomyopathy and Ventricular Arrhythmia. Heart rhythm. 16, 1669–1675. doi:10.1016/j.hrthm.2019.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang Y., Wang J., Xing Y., Cui C., Cheng H., Chen Z., et al. (2022). Dynamics of Cardiac Autonomic Responses during Hemodialysis Measured by Heart Rate Variability and Skin Sympathetic Nerve Activity: The Impact of Interdialytic Weight Gain. Front. Physiol. 13, 890536. doi:10.3389/fphys.2022.890536

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao M., Guan L., Wang Y. (2020). The Association of Autonomic Nervous System Function with Ischemic Stroke, and Treatment Strategies. Front. Neurol. 10, 1411. doi:10.3389/fneur.2019.01411

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao Z., Chen D., Wang L., Han C. (2018). Credit Risk Diffusion in Supply Chain Finance: a Complex Networks Perspective. Sustainability 10, 4608. doi:10.3390/su10124608

CrossRef Full Text | Google Scholar