A new directionality index based on high-resolution joint symbolic dynamics to assess information transfer in multivariate networks

The detection and quantification of coupling strength and direction are important aspects for achieving a deeper understanding of physiological regulatory processes in the field of network physiology. Due to the limitations of established approaches, we developed directionality indices based on simple mathematical symbolization principles and simple computational procedures that allow a quick and comprehensive understanding of the underlying couplings. We introduced a new directionality index (D_HRJSD) derived from the pattern family density matrix of the High-Resolution Joint Symbolic Dynamics (HRJSD) approach and its multivariate version (mHRJSD) to determine coupling direction and driver-response relationships. The mHRJSD approach contains the multivariate directionality index D_mHRJSD (D_mHRJSD(x,y|z), D_mHRJSD(x,z|y), and D_mHRJSD(y,z|x)), allowing us to determine the primary driver **D_mHRJSD, the secondary driver *D_mHRJSD, and the dominant responder ⁻D_mHRJSD in multivariate systems that are at least weakly coupled. Different linear and non-linear bi- and multivariate coupled systems (Gaussian autoregressive models) with different mutual influences were generated to validate these indices. The simulation results showed that D_HRJSD was able to correctly detect the dominant coupling direction in linear bivariate coupled systems but was partly able to detect the dominant coupling direction in non-linear bivariate coupled systems. The proposed directionality index D_mHRJSD derived from the mHRJSD approach is able to correctly detect the driver-responder relationships in linear coupled systems. The main advantages of the newly introduced directionality indices include their insensitivity to non-stationary time series, their ability to capture couplings through a simple, fast, and easy-to-implement symbolization procedure, and their scale invariance. Additionally, they are independent of time series length, model order selection, and the procedure for determining their significance level.

1 Introduction

The complex interplay between the central nervous system (CNS) and autonomic nervous system (ANS) with their large number of subsystems (parasympathetic and sympathetic activity) is also known as the central autonomic network (CAN) (Bashan et al., 2012; Bartsch et al., 2015; Ivanov et al., 2016). It has been shown that the output of CAN is directly linked to ANS (heart rate) as well and that sensory information from peripheral end organs provides feedback to the CAN (i.e., baroreceptor reflex). Information transfer between the CNS and ANS operates as a feedback-feedforward system, dynamically responding to the body’s significant demands. These brain-heart interactions are involved in multiple bodily processes, including sensing, integration, and activity regulation, to maintain homeostasis (Craig, 2002). Communication between the brain and the heart is bidirectional and occurs through different neural mechanisms, such as the vagal and spinal pathways (Chen et al., 2021). Moreover, the relationships between cortical network segregations and cardiac dynamics as cardiac sympathetic–vagal oscillations may provide valuable insights into the affective state of healthy and diseased-related alterations in network physiology (Candia-Rivera et al., 2024).

The role of the cerebral cortex in autonomic control of the cardiovascular system is gaining increased attention in medicine. Different cardiovascular control centers in the brainstem deal with different reflex mechanisms of cardiovascular adjustment (i.e., the cardiopulmonary reflex, the chemoreflex, and the baroreflex) (Dampney, 1994). Here, neurons in the caudal and rostral ventrolateral medulla affect efferent sympathetic reflexes and contribute to the maintenance of heart rate and blood pressure via the intermediolateral cell column of the spinal cord. The two medullary areas, the ambiguous nucleus and the vague nerve’s dorsal motor nucleus, are preganglionic parasympathetic neurons mediating the efferent parasympathetic reflex mechanism (McAllen, 1976; Taylor et al., 2001).

The parasympathetic nervous system is responsible for the “rest and digest” function while sitting, resting, and relaxing. It constricts the pupils, slows the heart rate and contractility, contracts the bronchial musculature, stimulates bronchial secretions, and enhances gut motility for digestion. The preganglionic neurons synapse with postganglionic neurons in the parasympathetic ganglion located next to or in the effector end organs. The sympathetic nervous system dominates during “fight-or-flight” reactions and during exercise and thus prepares the body for stressful physical activity. Sympathetic nervous activity increases the flow of well-oxygenated blood and is rich in nutrients to the tissues that need it, particularly the working skeletal muscles. The preganglionic sympathetic neurons arise from the thoracic and lumbar regions of the spinal cord (segments T1 through L2) and are located about halfway between the CNS and the effector tissue (McCorry, 2007). The preganglionic neurons of both the sympathetic and parasympathetic divisions release the neurotransmitter acetylcholine. The postganglionic neurons of the parasympathetic system also release acetylcholine, whereas the postganglionic sympathetic neurons release norepinephrine (Rea, 2016). The cardiac or respiration-related activity (parasympathetic) is connected to preganglionic neurons. It has been shown that brain regions like the insula, thalamus, hypothalamus, amygdala, and medial prefrontal cortex are involved in autonomic regulation at rest and during cognitive or emotional stress conditions, proven by functional brain imaging (Ziegler et al., 2009; Shoemaker et al., 2015). Beissner et al. (2013) showed that largely divergent brain networks were associated with sympathetic and parasympathetic activity. The ventromedial prefrontal cortex (VMPFC), the perigenual anterior cingulate cortex (pACC), the dorsal anterior cingulate cortex (dACC), the posterior cingulate cortex (PCC), the insular cortices, and the amygdala seem to be the main cortical and subcortical areas involved in ANS regulation processes that are created by a network of interactions related to the task and autonomic division.

Investigating the coupling between these ANS subsystems, their variability, and brain activity may lead to a better understanding of pathophysiological regulatory processes within the central autonomic network. For the quantitative analysis of the brain-heart (CNS-ANS) network coupling pathways and its integrated interacting subsystems, such as the cardiovascular and cardiorespiratory systems, several linear/non-linear univariate and multivariate approaches are available.

These approaches focus on characterizing the multivariate information transfer. These concepts (Schulz et al., 2013a; Bartsch et al., 2015; Faes et al., 2015; Ivanov et al., 2016) are applicable in the following domains: entropy, Granger causality, non-linear prediction, phase synchronization, symbolization, recurrence quantification analysis (RQA), and functional connectivity analysis techniques (Marwan et al., 2013; Aguirre et al., 2018; Lombardi et al., 2019). Studies have demonstrated that information transfer between the cardiovascular and cardiorespiratory systems exhibits strong non-linearity (Novak et al., 1993), and therefore, linear approaches alone cannot fully quantify physiological and pathophysiological regulatory processes. There is no generally superior approach capable of considering all aspects of coupling analysis (linearity, non-linearity, causality, multivariate analysis, directionality, coupling strength) and its quantitative evaluation. Some of these approaches include one or more of these aspects, but usually not to a sufficient extent, so the time series with their mutual interactions and couplings can only be interpreted and analyzed incompletely and in parts. Furthermore, many of these approaches are not standardized and not user-friendly (degrees of freedom, preconditions, model selection and model order estimation, scale dependency). They are based on purely mathematical concepts, making it difficult to select the “right” approach to apply them to quantify physiological and pathophysiological regulatory processes.

Computerized quantitative models are essential to integrating and evaluating the information generated by these methods. For example, the Physiome Project provides a framework for qualitatively evaluating such information. The term “Physiome” is derived from “physio,” meaning “life,” and “ome,” meaning “as a whole,” (Leem, 2016). The project developed a multiscale modeling framework to understand physiological functions, enabling models to be combined and linked hierarchically.

In this study, we employed the High-Resolution Joint Symbolic Dynamics (HRJSD) approach, which represents an enhanced version of the classical Joint Symbolic Dynamics (JSD). The HRJSD approach was initially developed for the analysis of non-linear cardiovascular and cardiorespiratory couplings in acute schizophrenia (Schulz et al., 2013b; Schulz et al., 2013c; Schulz et al., 2017a; Schulz et al., 2017b). This was based on analyzing dynamic processes through symbols (Voss et al., 1996). The HRJSD is founded upon a symbolization procedure that enables a coarse-grain quantitative evaluation of the classification and characterization of short-term regulatory bivariate coupling patterns that are predominant in the interaction generated by the ANS. The HRJSD employs a redundancy reduction strategy that groups single-word types into eight pattern families, thereby enabling a comprehensive quantification of bivariate short-term autonomic coupling patterns. Based on the concept of frequent deterministic pattern classification, the bivariate redundancy reduction strategy overcomes the limitations of classical univariate symbolization strategies. It facilitates a transition between univariate and bivariate symbolic analyses, significantly advancing over the standard JSD and Symbolic Coupling Traces (SCT) (Wessel et al., 2009, 2011). The HRJSD method overcomes the issues encountered by JSD in differentiating between decreases and steady states, as well as between small and large changes in autonomic regulation due to threshold level l = 0 and the alphabet A = {0,1} for symbol transformation. It is impossible to differentiate between noise, artifacts (for example, those generated by undersampling or ectopic events), and fluctuations arising from (patho) physiological regulatory processes when using JSD. However, both approaches have the main advantages of being insensitive to non-stationary time series and capable of capturing non-linear couplings through a straightforward procedure.

Therefore, the main objective of this study was to develop new directionality indices derived from coupling approaches that are based on simple mathematical principles, such as symbolization, and simple calculation procedures, enabling a comprehensive understanding of the underlying couplings in a fast and easy way, and do not have the limitations of already established approaches. Here, we developed new directionality indices enabling the assessment of coupling directions in bivariate and multivariate systems based on the high-resolution joint symbolic dynamics approach (HRJSD).

2 Materials and methods

2.1 Basics of high-resolution joint symbolic dynamics – HRJSD

Baumert et al. (2002) developed the joint symbolic dynamics (JSD) method to analyze nonlinear couplings between systolic blood pressure (SP) and heart rate (BBI) time series, relying on the analysis of dynamic processes using symbols (Kurths et al., 1995).

Schulz et al. (2013c) introduced an enhanced version of the classical JSD to overcome their restrictions: high-resolution joint symbolic dynamics (HRJSD), which is based on three symbols and a symbol-transformation threshold, which can be used to quantify short-term non-linear coupling by means of symbols.

Therefore, the direct analysis of successive signal amplitudes is based on discrete states (symbols). In short, HRJSD works by transforming the two investigated time series (here: BBI and SP) into symbol sequences based on their signal amplitudes using a given alphabet A = {0, 1, 2}. The bivariate sample vector X (Equation 1) of the two-time series with x_BBI and x_SP is transformed into a bivariate symbol vector S (Equation 2), where n are the nth beat-to-beat values of BBI and SP, respectively.

$X={\left\{{\left[x_n^{BBI},x_n^{SP}\right]}^T\right\}}_{n=0,1,\dots }x\in R\left(\begin{array}{l}R\mathrm{is\; the\; subset\; of}\\ \mathrm{real\; positive\; numbers}\end{array}\right)\text{ (1)}$

$\begin{array}{ll}S={\left\{{\left[s_n^{BBI},s_n^{SP}\right]}^T\right\}}_{n=0,1,\dots }s\in 0,1,2& (2)\end{array}$

The bivariate symbol vector S is defined using the following definitions (Equations 3, 4):

$S_n^{BBI\text{ (3)}}=\{\begin{array}{l}0:\left(x_{n+1}^{BBI}-x_n^{BBI}\right)<-l^{BBI}\hfill \\ 1:-l^{BBI}\le \left(x_{n+1}^{BBI}-x_n^{BBI}\right)\le l^{BBI}\hfill \\ 2:\left(x_{n+1}^{BBI}-x_n^{BBI}\right)>l^{BBI}\hfill \end{array}$

$S_n^{SP\text{ (4)}}=\{\begin{array}{l}0:\left(x_{n+1}^{SP}-x_n^{SP}\right)<-l^{SP}\hfill \\ 1:-l^{SP}\le \left(x_{n+1}^{SP}-x_n^{SP}\right)\le l^{SP}\hfill \\ 2:\left(x_{n+1}^{SP}-x_n^{SP}\right)>l^{SP}\hfill \end{array}$

Increasing values were coded as “2,” decreasing values as “0,” and unchanging (no or little variability) values as “1,” respectively. Afterward, S was subdivided into short words (sequences of symbols, bins) w_k of length k = 3. In this study, an adapted threshold for the individual physiological dynamic variability was equal to 25% of the standard deviation of the time series. The derived different word types from the BBI (w_BBI) and SP ( $w_{SP}$ ) time series (word types ranging from: 000, 001,…, 221, 222) were organized into a normalized 27 × 27 vector matrix W_n ranging from word type (000,000)^T to (222,222)^T. These single-word types $w_{\mathrm{B}\mathrm{B}\mathrm{I}},w_{SP}$ (total number of all word type combinations 27 × 27 = 729) were afterward grouped into eight pattern families’ wf, whereby the sum of probabilities of all single-word family occurrences p(wf) was normalized to 1. The eight pattern families (E0, E1, E2, LU1, LD1, LA1, P, V) describe different aspects of autonomic modulation of the BBI- and SP time series and were sorted into an 8 × 8 pattern family density matrix Wf, resulting in 64 coupling patterns. The pattern famines are defined as follows:

• E0, E1, and E2: Words consisting of three equal symbols (no variation of symbols) of type ‘0,’ ‘1’, and ‘2,’ respectively.

• LU1 and LD1: Words consisting of two different symbols with low increasing behavior (LU1) and low decreasing behavior (LD1).

• LA1: Words consisting of two different alternating symbols of type ‘0’ and ‘2’ with an increasing-decreasing behavior.

• P and V: Words consisting of three different symbols with peak-like behavior (P) and valley-like behavior (V) (Schulz et al., 2013c)

From the word distribution density matrix Wf, the normalized joint probability of the occurrence of each word was estimated (Schulz et al., 2013b; Schulz et al., 2015).

From the matrix Wf, the sum of each (n = 8) column $c{f}_{SP}$ (cfE0, cfE1, cfE2, cfLU1, cfLD1, cfLA1, cfP, cfV) and the sum of each (n = 8) row rf_BBI (rfE0, rfE1, rfE2, rfLU1, rfLD1, rfLA1, rfP, rfV) were calculated (Figure 1).

Figure 1

Figure 1. Basic principle of HRJSD. (a) Transformation of the bivariate sample vector X (BBI = beat-to-beat intervals [msec]; SP = systolic blood pressure [mmHg]) into the bivariate symbol vector S (0: decreasing values, 1: equal values, 2: increasing values) and word distribution density matrix W_n (27 × 27). (b) Word pattern family distribution density matrix Wf (8 × 8) with eight pattern families wf created from 27 single-word types w_{BBI, SP}. Rows represent pattern families of BBI interval changes, and column pattern families of SP changes, rf_BBI (row): sum of the specific word family, cf_SP (column): sum of the specific word family.

2.2 Directionality index – bivariate system

To evaluate the physiological states of highly complex biological systems, it is important and necessary to determine synchronization processes within coupled complex systems and the predominant direction of their coupling. Causality can be defined by using the directionality of time to establish a causal ordering of dependent time series, encompassing both direct and indirect influences from one process to another. This definition can be applied to bivariate (two-time series, Figure 2A) and multivariate (more than two-time series, Figures 2B, 3) analysis. In the case of multivariate analysis, a distinction can be made between direct coupling (from one-time series to another) and indirect coupling (effects mediated through one or more other time series) (Schulz et al., 2013a).

Figure 2

Figure 2. Examples of directional dependencies for direct and indirect couplings. Interdependence structure for (a) a bivariate and (b) a multivariate case. (a) Direct coupling exists for x₁↔x₂; (b) direct coupling exists for x₁ → x₂ and x₂↔x₃ and indirect coupling between x₁ → x₃ mediated by x₂ (direction of coupling: →,← unidirectional, ↔ bidirectional). Adopted to Schulz et al. (2013a).

Figure 3

Figure 3. Simulated multivariate systems with their mutual influence between the time series x₁, x₂, and x₃. Arrows indicating the causal coupling direction from one system to another (e.g. x₁ → x₂ means a unidirectional driving (a, b) from system 1 (x₁) to system 2 (x₂), and x₂ → x₃ means a bidirectional driving (c) between system 2 (x₂) to system 3 (x₃)).

The coupling direction can be determined, e.g., from the amplitudes of the system (properties of the system state) by calculating their mutual predictability (Schiff et al., 1996), from mutual nearest neighbors (Arnhold et al., 1999; Quiroga et al., 2000) in the reconstructed state space, or by applying information theoretical approaches (Schreiber, 2000; Palus et al., 2001; Palus and Stefanovska, 2003).

In the field of symbolization, so far, no approaches have been available to determine the coupling direction for bivariate or multivariate systems. The first attempts at this were integrated into Symbolic Coupling Traces (SCT) (Wessel et al., 2011). SCTs can detect delayed couplings (time lags) but cannot assess the coupling direction or the driver-response relationships. To address this gap, we introduced a directionality index (D_HRJSD) derived from the 8 × 8 pattern family density matrix Wf from the HRJSD approach. This index is able to determine the dominant coupling direction and assess the driver-response relationships in bivariate (n = 2) (Figure 2A) and multivariate (n = 3) systems (Figure 2B).

For a bivariate system (x,y), the columns cf_x (n = 8) and the rows rf_y (n = 8) from the matrix Wf (Figure 1) were used to calculate D_HRJSD(x,y) (Equation 5):

$\begin{array}{ll}{D}_{\mathrm{HRJSD}}\left(x,y\right)=-\left({\displaystyle \sum }_{i=1}^n\frac{c{f}_{\mathrm{x}}\left(i\right)-r{f}_{\mathrm{y}}\left(i\right)}{c{f}_{\mathrm{x}}\left(i\right)+r{f}_{\mathrm{y}}\left(i\right)}\right)/n& (5)\end{array}$

If D_HRJSD(x,y) is positive, driving (→) from system 1 (x) to system 2 (y) predominates (Equation 6) and becomes negative for the opposite case (Equation 7).

$\begin{array}{ll}{D}_{\mathrm{HRJSD}}\left(x,y\right)>0;x\to y& (6)\end{array}$

$\begin{array}{ll}{D}_{\mathrm{HRJSD}}\left(x,y\right)<0;y\to x& (7)\end{array}$

2.2.1 Simulated coupled linear and non-linear systems to validate D_HRJSD

Simulated data were used to validate D_HRJSD. Therefore, two different multivariate models were applied (Baccala and Sameshima, 2001; Montalto et al., 2014), each with 100 simulated time series:

• Linear time series with a normal distribution of the variables, generated by a linear Gaussian AR model and

• Non-linear time series generated by a non-linear Gaussian AR model.

Three distinct multivariate coupled systems were generated for both linear and nonlinear models, incorporating varying mutual influences (unidirectional and bidirectional) between the time series (Figure 3).

The following equations were used for the three linear Gaussian autoregressive models (Baccala and Sameshima, 2001; Montalto et al., 2014) (Equations 8–10):

Linear system 1, LS1 (Figure 3A):

$\begin{array}{c}{x}_1\left(n\right)=0.95\sqrt{2}{x}_1\left(n-1\right)-0.9025x_1\left(n-2\right)+w_1\left(n\right)\\ x_2\left(n\right)=-0.5x_1\left(n-1\right)+w_2\left(n\right)\\ x_3\left(n\right)=0.4x_1\left(n-2\right)+w_3\left(n\right)\end{array}\text{ (8)}$

Linear system 2, LS2 (Figure 3B):

$\begin{array}{c}{x}_1\left(n\right)=0.95\sqrt{2}{x}_1\left(n-1\right)-0.9025x_1\left(n-2\right)+w_1\left(n\right)\\ x_2\left(n\right)=0.5x_1\left(n-2\right)+w_2\left(n\right)\\ x_3\left(n\right)=-0.4x_1\left(n-3\right)-0.2x_2\left(n-2\right)+w_3\left(n\right)\end{array}\text{ (9)}$

Linear system 3, LS3 (Figure 3C):

$\begin{array}{c}{x}_1\left(n\right)=0.95\sqrt{2}{x}_1\left(n-1\right)-0.9025x_1\left(n-2\right)+w_1\left(n\right)\\ x_2\left(n\right)=0.5x_1\left(n-2\right)+0.4x_3\left(n-1\right)+w_2\left(n\right)\\ x_3\left(n\right)=-0.4x_1\left(n-3\right)-0.2x_2\left(n-2\right)+w_3\left(n\right)\end{array}\text{ (10)}$

Where w₁(n), w₂(n), and w₃(n) were drawn from Gaussian noise with zero mean and unit variance. For linear system 3, a closed loop from x₃(n) back to x₂(n) via a direct connection was integrated, with x₃ as the predominant driver.

For the non-linear models (Montalto et al., 2014) (Equations 11–13), x₂(n) was modified by a quadratic term of $x_1^2$ . Thus, the three linear model equations changed to:

Non-linear system 1, NLS1 (Figure 3A):

$\begin{array}{c}{x}_1\left(n\right)=0.95\sqrt{2}{x}_1\left(n-1\right)-0.9025x_1\left(n-2\right)+w_1\left(n\right)\\ x_2\left(n\right)=-0.5x_1^2\left(n-1\right)+w_2\left(n\right)\\ x_3\left(n\right)=0.4x_1\left(n-2\right)+w_3\left(n\right)\end{array}\text{ (11)}$

Non-linear system 2, NLS2 (Figure 3B):

$\begin{array}{c}{x}_1\left(n\right)=0.95\sqrt{2}{x}_1\left(n-1\right)-0.9025x_1\left(n-2\right)+w_1\left(n\right)\\ x_2\left(n\right)=0.5x_1^2\left(n-2\right)+w_2\left(n\right)\\ x_3\left(n\right)=-0.4x_1\left(n-3\right)-0.2x_2\left(n-2\right)+w_3\left(n\right)\end{array}\text{ (12)}$

Non-linear system 3, NLS3 (Figure 3C):

$\begin{array}{c}{x}_1\left(n\right)=0.95\sqrt{2}{x}_1\left(n-1\right)-0.9025x_1\left(n-2\right)+w_1\left(n\right)\\ x_2\left(n\right)=0.5x_1^2\left(n-2\right)+0.5x_3\left(n-1\right)+w_2\left(n\right)\\ x_3\left(n\right)=-0.4x_1\left(n-3\right)-0.2x_2\left(n-2\right)+w_3\left(n\right)\end{array}\text{ (13)}$

Where w₁(n), w₂(n), and w₃(n) were drawn from Gaussian noise with zero mean and unit variance. For non-linear system 3, a closed loop from x₂(n) back to x₃(n) via a direct connection was integrated, with x₂ as the predominant driver.

2.3 Directionality index – multivariate system

No methodical approach based on symbolization allows the detection of the coupling direction in multivariate coupled systems. Rather, no coupling approach is available that can determine the primary driver, the secondary driver, and the dominant responder in multivariate weakly coupled systems. To overcome the limitation of analyzing bivariate couplings, only the HRJSD approach was adopted in a further step for the quantification of multivariate couplings – the multivariate High-Resolution Joint Symbolic Dynamics (mHRJSD) (Schulz et al., 2018). In short, mHRJSD works in that way, that the set of three investigated time series (e.g., X, Y, and Z) was transformed into symbol sequences based on their signal amplitudes using a given alphabet A = {0, 1, and 2}. The trivariate sample vector X of these time series, x^X, x^Y and x^Z, were then transformed into a trivariate symbol vector S, where n were the nth beat-to-beat values of X, Y, and Z, respectively. Here, all single-word types w_{X, Y, and Z} were grouped into eight pattern families’ w_f whereby the probabilities of all single-word family occurrences p(w_f) were also normalized to 1. These eight pattern families were sorted into an 8 × 8 × 8 pattern family density matrix Wf. Furthermore, from the matrix Wf, the sum of each (n = 8) x-, y-, and z-plane (pf_X, pf_Y, pf_Z) as pfE0, pfE1, pfE2, pfLU1, pfLD1, pfLA1, pfP, and pfV were calculated describing how one family pattern in one time series is coupled with all other eight pattern families of the other two time series.

For the mHRJSD approach, the proposed directionality index was further extended to determine the dominant coupling direction and assess the driver-response relationships in multivariate (n = 3) systems.

For a multivariate system (x,y,z), the single-word family occurrences p(w_f) from the x-, y- and z-plane (pf_X, pf_Y, pf_Z) from the 8 × 8 × 8 pattern family density matrix Wf were used to calculate D_mHRJSD (Equations 14–16). Thereby, for each coupling pathway, one directionality index was calculated (e.g., two interacting time series: x and y with z as the covariate |). Thus, for the coupling between the time series x and y with covariate z, the directionality index (Equation 14) is defined as:

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(x,y|z\right)=-\left({\displaystyle \sum }_{i=1}^n\frac{p{f}_{\mathrm{x}|\mathrm{z}}\left(i\right)-p{f}_{\mathrm{y}|\mathrm{z}}\left(i\right)}{p{f}_{\mathrm{x}|\mathrm{z}}\left(i\right)+p{f}_{\mathrm{y}|\mathrm{z}}\left(i\right)}\right)/n& (14)\end{array}$

For the coupling between the time series x and z with covariate y, the directionality index (Equation 15) is defined as:

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(x,z|y\right)=-\left({\displaystyle \sum }_{i=1}^n\frac{p{f}_{\mathrm{x}|\mathrm{y}}\left(i\right)-p{f}_{\mathrm{z}|\mathrm{y}}\left(i\right)}{p{f}_{\mathrm{x}|\mathrm{y}}\left(i\right)+p{f}_{\mathrm{z}|\mathrm{y}}\left(i\right)}\right)/n& (15)\end{array}$

For the coupling between the time series y and z with covariate x, the directionality index (Equation 16) is defined as:

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(y,z|x\right)=-\left({\displaystyle \sum }_{i=1}^n\frac{p{f}_{\mathrm{y}|\mathrm{x}}\left(i\right)-p{f}_{\mathrm{z}|\mathrm{x}}\left(i\right)}{p{f}_{\mathrm{y}|\mathrm{x}}\left(i\right)+p{f}_{\mathrm{z}|\mathrm{x}}\left(i\right)}\right)/n& (16)\end{array}$

If D_mHRJSD(x,y|z) is positive, driving (→) from system 1 (x) to system 2 (y) predominates (Equation 17) and becomes negative for the opposite case (Equation 18).

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(x,y|z\right)>0;x|z\to y|z& (17)\end{array}$

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(x,y|z\right)<0;y|z\to x|z& (18)\end{array}$

If D_mHRJSD(x,z|y) is positive, driving (→) from system 1 (x) to system 2 (z) predominates (Equation 19) and becomes negative for the opposite case (Equation 20).

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(x,z|y\right)>0;x|y\to z|y& (19)\end{array}$

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(x,z|y\right)<0;z|y\to x|y& (20)\end{array}$

If D_mHRJSD(y,z|x) is positive, driving (→) from system 1 (y) to system 2 (z) predominates (Equation 21) and becomes negative for the opposite case (Equation 22).

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(y,z|x\right)>0;y|x\to z|x& (21)\end{array}$

$\begin{array}{ll}{D}_{\mathrm{mHRJSD}}\left(y,z|x\right)<0;z|x\to y|x& (22)\end{array}$

Thus, three indices were derived, which are subsequently used to determine the strongest driver and the most dominant responder in the overall system. Therefore, all three indices were compared to whether they were greater (+) or less than 0 (−) (Table 1). This means that if the index is greater than zero for the multivariate interaction, “+” is set; if the index is less than zero, “−” is set. Afterward, the sum of the three comparisons of the three directionality indices was determined and counted. If a time series is present twice (++ or −− or +− or −+) as a driver, it dominates the overall system and is determined as the primary driver **D_mHRJSD, and if a time series is present only once (+ or −), then it is determined as the secondary driver *D_mHRJSD of the overall system, and the non-occurring time series is the dominant responder ⁻D_mHRJSD.

Table 1

Table 1. Determination of the primary driver (**D_mHRJSD), secondary driver (*D_mHRJSD), and the dominant responder (⁻D_mHRJSD) in a multivariate system derived from the directionality indices D_mHRJSD(x,y|z), D_mHRJSD(x,z|y), and D_mHRJSD(y,z|x).

For exapmple: D_mHRJSD(x,y|z) > 0, D_mHRJSD(x,z|y) > 0, and D_mHRJSD(y,z|x) > 0, resulting in (+++ equal to x-x-y) or D_mHRJSD(x,y|z) < 0, D_mHRJSD(x,z|y) > 0, and D_mHRJSD(y,z|x) < 0, resulting in (−+− equal to y-x-z). The indices **D_mHRJSD, *D_mHRJSD, and ⁻D_mHRJSD are determined by their absolute values in descending order of importance of their values. For the first example (+++), we have information flows from x → y, x → z, and y → z resulting in x as **D_mHRJSD, y as *D_mHRJSD, and z as ⁻D_mHRJSD. For the second example (−+−), we have information flows from y → x, x → z, and z → y, resulting in a closed-loop where **D_mHRJSD, *D_mHRJSD, and ⁻D_mHRJSD cannot be clearly determined (these options are represented by the last two lines in the Table 1).

The simulated linear and non-linear AR systems were validated with two additional approaches: the normalized short-time partial directed coherence (NSTPDC) (Adochiei et al., 2013) and the multivariate transfer entropy (MuTE) (Montalto et al., 2014). Both methods allow us to determine the coupling direction. NSTPDC mainly detects linear coupling, whereas MuTE mainly detects non-linear coupling. In short, NSTPDC is based on an m-dimensional AR model with the order p and allows determining linear Granger causality in the frequency domain. Mute is an information-theoretical approach that detects the information transfer between multivariate joint processes and discovers purely non-linear interactions with a range of interaction delays. We applied in-house software in the programming environment MatlabR2013b.

3 Results

3.1 Results of simulated systems to validate D_HRJSD

All three methods, HRJSD, NSTPDC, and MuTE, calculated a directionality index D (D_HRJSD, D_NSTPDC, D_MuTE), which was used for validation. These three indices have in common that if the index is positive, driving (→) from system 1 (x) to system 2 (y) predominates and becomes negative for the opposite case that system 2 (y) is driving system 1 (x) (Table 2). Table 2 presents the results of the simulated linear and non-linear AR models with the underlying simulated driver-response relationships in a multivariate context and the results of the applied directionality indices.

Table 2

Table 2. Results of simulated linear and non-linear autoregressive (AR) systems to validate the directionality index D_HRJSD (blue: driver variable).

For example, here, the first two rows are explained where the information flows were simulated with a linear AR model, linearly coupled variables, and simulated driver-response relationships 1 → 2 and 1 → 3 (blue represents the simulated driver variable). That means time series 1 transfers information to time series 2 and 3. The coupling directions are from 1 to 2 and 1 to 3.

Linear system 1:

1 → 2 and 1 → 3: D_HRJSD, D_NSTPDC, and D_MuTE are positive; correct classification of the predominating coupling directions (1 is the driver).

Linear system 2:

1 → 2, 1 → 3, 2 → 3: D_HRJSD, D_NSTPDC, and D_MuTE are positive; correct classification of the predominating coupling directions (1 and 2 are drivers).

Linear system 3:

1 → 2, 1 → 3: D_HRJSD, D_NSTPDC, and D_MuTE are positive; correct classification of the predominating coupling directions (1 is the driver).

2⇄3: D_HRJSD, D_NSTPDC, and D_MuTE are negative; correct classification of the predominating coupling direction (3 is the driver).

For the linear AR model with purely linear couplings among the three variables (1, 2, 3), all directionality indices (D_HRJSD, D_NSTPDC, D_MuTE) were able to correctly detect the predominating coupling directions and the related driver variable.

Non-linear system 1 (Table 2):

1 → 2: D_NSTPDC and D_MuTE are positive; correct classification of the predominating coupling direction (1 is the driver).

D_HRJSD is negative; incorrect classification of the predominating coupling direction. D_HRJSD detects variable 2 as the driver.

1 → 3: D_HRJSD, D_NSTPDC, and D_MuTE are positive; correct classification of the predominating coupling direction (1 is the driver).

Non-linear system 2:

1 → 2, 1 → 3, 2 → 3: D_NSTPDC and D_MuTE are positive; correct classification of the predominating coupling directions (1 and 2 are drivers).

D_HRJSD is negative; incorrect classification of the predominating coupling directions. D_HRJSD detects variables 2 and 3 as the drivers.

Non-linear system 3:

1 → 2, 1 → 3, 2 ⇄3: D_NSTPDC and D_MuTE are negative; correct classification of the predominating coupling directions (1 and 2 are drivers).

D_HRJSD is negative; incorrect classification of the predominating coupling directions. D_HRJSD detects variables 2 and 3 as the drivers.

For the non-linear AR model with purely non-linear couplings among the three variables (1, 2, 3), only NSTPDC and MuTE were able to correctly detect the predominating coupling directions and the related driver variable. D_HRJSD could partly detect the dominating coupling direction in non-linear systems (non-linear system 1).

The proposed directionality index D_HRJSD derived from the HRJSD approach is able to correctly detect the dominating coupling direction in linear bivariate coupled systems but is only partly able to detect the dominating coupling direction in non-linear bivariate coupled systems. Due to this limitation, in detailed investigations to determine the coupling direction, other methods should be used in addition to D_HRJSD (e.g., MuTE), which can also correctly determine the dominant driver-response relationships in pure non-linear systems.

3.2 Results of simulated systems to validate D_mHRJSD

Similarly to validating D_HRJSD for the bivariate system, simulated data was used to validate D_mHRJSD. Therefore, a multivariate linear Gaussian AR model was applied to generate a set of multivariate linear time series (n = 100) with a normal distribution of the variables (Baccala and Sameshima, 2001; Montalto et al., 2014). For the linear model, two different multivariate coupled systems were generated (Equations 8, 9) with different mutual influences (unidirectional, bidirectional) between the time series (Figures 3B,C). Non-linear AR models were not applied since D_HRJSD seems to only partly detect the correct driver-responder relationship between non-linear coupled time series.

The results of the two multivariate coupled linear AR systems showed that the determination of the multivariate directionality index D_mHRJSD works properly, as well as the determination of the primary driver **D_mHRJSD, the secondary driver *D_mHRJSD, and dominant responder ⁻D_mHRJSD in the multivariate systems (Table 3).

Table 3

Table 3. Determination of the primary driver (**D_mHRJSD), secondary driver (*D_mHRJSD), and the dominant responder (⁻D_mHRJSD) derived from the directionality indices D_mHRJSD(x,y|z), D_mHRJSD(x,z|y), and D_mHRJSD(y,z|x) for two simulated multivariate coupled systems (LS2, LS3).

Linear system 2, LS2:

x₁ → x₂: D_mHRJSD(x₁,x₂|x₃) is positive; correct classification of the dominating coupling direction (x₁ is driver).

x₁ → x₃: D_mHRJSD(x₁,x₃|x₂) is positive; correct classification of the dominating coupling direction (x₁ is driver).

x₂ → x₃: D_mHRJSD(x₂,x₃|x₁) is positive; correct classification of the dominating coupling direction (x₂ is driver).

From this results that:

D_mHRJSD(x₁,x₂|x₃) = x₁ and D_mHRJSD(x₁,x₃|x₂) = x₁ ⇒ **D_mHRJSD = x₁.

D_mHRJSD(x₂,x₃|x₁) = x₂ ⇒ *D_mHRJSD = x₂.

⇒ ⁻D_mHRJSD = x₃.

For the coupled multivariate linear AR model (LS2), the correct driver-responder relationships were classified with x₁ as the primary driver, x₂ as the secondary driver, and x₃ as the responder of the system, as it was simulated.

Linear system 3, LS3:

x₁ → x₂: D_mHRJSD(x₁,x₂|x₃) is positive; correct classification of the dominating coupling direction (x₁ is driver).

x₁ → x₃: D_mHRJSD(x₁,x₃|x₂) is positive; correct classification of the dominating coupling direction (x₁ is driver).

x₂⇄x₃: D_mHRJSD(x₂,x₃|x₁) is negative; correct classification of the dominating coupling direction (x₃ is driver).

From this results that:

D_mHRJSD(x₁,x₂|x₃) = x₁ and D_mHRJSD(x₁,x₃|x₂) = x₁ ⇒ **D_mHRJSD = x₁.

D_mHRJSD(x₂,x₃|x₁) = x₃ ⇒ *D_mHRJSD = x₃.

⇒ ⁻D_mHRJSD = x₂.

For the coupled multivariate linear AR model (LS3), the correct driver-responder relationships were classified with x₁ as the primary driver, x₃ as the secondary driver, and x₂ as the responder of the system, as it was simulated.

The mHRJSD approach contains multivariate directionality indices D_mHRJSD (D_mHRJSD(x,y|z), D_mHRJSD(x,z|y), and D_mHRJSD(y,z|x)), allowing us to determine the primary driver **D_mHRJSD, the secondary driver *D_mHRJSD, and the dominant responder ⁻D_mHRJSD in multivariate systems. Therefore, it has to be assumed that the time series to be analyzed are at least weakly coupled with each other. Limiting factors are that the proposed directionality index D_mHRJSD derived from the mHRJSD approach can only correctly detect the driver-responder relationships in linear coupled systems and cannot detect the driver-responder relationships in non-linear coupled systems. The mHRJSD approach can evaluate direct causal information transfer in multivariate systems. Despite this limitation of D_mHRJSD, the feature to assess the driver-response relationships in multivariate systems is not implemented in any of the existing symbolization approaches and thus clearly complements the existing coupling approaches.

4 Discussion

The HRJSD approach emphasizes a clear characterization of how the couplings are composed by regulatory aspects of the ANS; it is able to quantify the coupling direction (directionality index: D_HRJSD) in linear and non-linear coupled systems, which was not possible with existing symbolization approaches, neither for bivariate nor for multivariate systems, and assesses the driver-response relationships in bivariate (n = 2) and multivariate (n = 3) systems.

The newly developed directionality indices derived from (m)HRJSD are based on simple mathematical symbolization principles and simple calculation procedures, enabling a comprehensive understanding of the underlying couplings in a fast and easy way, and do not have the limitations of existing approaches. The main advantages of the (m) HRJSD-derived directionality indices are that they are insensitive to non-stationary time series; they are able to capture couplings through a simple, fast, and easy-to-implement symbolization procedure; they are scale-invariant; they are independent of time series length, model order selection, and significance level determination procedure. Moreover, clear advantages of using D_HRJSD over, e.g., MuTE are that with D_HRJSD, a multivariate interaction can be classified in the overall system in such a way that it can be determined which variable in the system generates the primary and secondary information flow and which variable only acts as a responder in the system. This assessment is currently only possible for linear systems. Further research with other models and time delays will provide further insight into D_HRJSD. Most of the already established coupling approaches, e.g., partial directed coherence (PDC) or directed transfer function, depend on the reliability of the fitted multivariate autoregressive model (MAR) (i.e., optimal model order, epoch length) and a significance level has to be used for both to avoid spurious interactions (Schulz et al., 2013a). In general, most of these coupling approaches have high degrees of freedom, are not standardized in their preconditions (e.g., preprocessing steps, parameter settings, time-series length, model order selection, significance level determination, scale-independent data, and stationarity), and have been validated using physiological and pathophysiological cases (Schulz et al., 2013a).

The HRJSD approach includes different threshold levels and a directionality index, D_HRJSD. The validation studies showed that the directionality index D_HRJSD is able to correctly detect the dominating coupling direction in linear coupled systems but is partly able to detect the dominating coupling direction in non-linear coupled systems. It is an intriguing observation that, within non-linear systems, where linear couplings may also prevail, the method can precisely detect these linear couplings. This suggests that the method can be employed to identify the part of the underlying coupling in a system, irrespective of its linear or non-linear nature and regardless of the type (linear or non-linear) of coupling between the variables (Table 2, NLS1, 1→3).

The HRJSD approach and other symbolization approaches are only able to analyze bivariate couplings, whereas the HRJSD approach is also able to determine the driver-responder relationship; facing this, the mHRJSD approach is able to quantify multivariate couplings and to determine dominant driver-responder relationship in multivariate coupled systems. These are outstanding new features for coupling analyses based on symbolizations. The mHRJSD approach facilitates multivariate analysis by incorporating a third time series, enabling the coarse-grained evaluation of time series dynamics. As outstanding and unique features of the mHRJSD approach are the implemented multivariate directionality indices D_mHRJSD (D_mHRJSD(x,y|z), D_mHRJSD(x,z|y), and D_mHRJSD(y,z|x)) allowing us to determine the primary driver **D_mHRJSD, the secondary driver *D_mHRJSD, and the dominant responder ⁻D_mHRJSD in multivariate systems (assumption: weakly coupled system). The simulation procedure revealed that the proposed directionality index D_mHRJSD derived from the mHRJSD approach is able to correctly detect the driver-responder relationships in linear coupled systems. Moreover, the mHRJSD approach is able to evaluate the direct causal information transfer in multivariate systems.

However, given the limitation of D_HRJSD (e.g., MuTE) in fully determining coupling direction in nonlinear coupled systems, additional methods should be employed. The method-specific characteristics of the different coupling approaches operate in different domains in the assessment of coupling (strength and direction) and causality, and there is currently no superior approach that can combine all the advantages in a single approach. While non-linear methods study complex signal interactions, linear methods favor the frequency domain representation of biological signals (characterizing the connectivity between specific oscillatory components) (Schulz et al., 2013a).

Due to this restriction so far, it would be of great interest to test other models and see how the newly proposed directional indices are able to evaluate driver-response relationships. Here, the physiome as a computerized quantitative model could be an option. The physiome is the definitive quantitative and integrated description of the functional behavior of the physiological state of an individual or species. The physiome describes the physiological dynamics of the normal, intact organism. It is built upon information and structure, namely the genome, proteome, and morpheme. The physiome must define relationships from genome to organism and from functional behavior to gene regulation. Quantitative models in physiome are unlike any other database-driven research area, including bioinformatics, network biology, or big data analysis. The physiomic model is a repository of previous data and a tool for testing and predicting results by varying factors. Thereby, models can be defined at various levels of abstraction: the conceptual level, the mathematical level, the formulation level, and the solution level (Leem, 2016). Physiome and network physiology are very similar as they pursue integrative concepts. The multidisciplinary field of network physiology studies how different physiological systems and subsystems interact to regulate their functions. It looks at how these systems interact across different spatial and temporal scales, from cellular to organism levels. It also looks at how they communicate and work together to generate different physiological states and behaviors in health and disease (Bartsch et al., 2015b; Ivanov et al., 2016a; Ivanov, 2021).

The cardiovascular, cardiorespiratory, and central networks are intricate physiological systems with both direct and indirect interactions. In the investigation of these networks, bivariate approaches are frequently employed. Nevertheless, it seems reasonable to posit that multivariate approaches will become more prevalent in place of bivariate ones, given that they enhance the characterization of causal or non-causal interrelationships between the networks. For instance, multivariate coupling analysis, for example, heart rate, systolic blood pressure, respiration, and central activity, may provide more advanced information about the complex autonomic network in physiological and pathophysiological conditions than uni- and bivariate approaches.

For instance, interactions within the cardiorespiratory network are primarily reflected in respiratory sinus arrhythmia (RSA), characterized by rhythmic fluctuation in cardiac cycle intervals in relation to respiration. Two principal mechanisms have been proposed to explain the RSA. The first is the central influence of respiration on vagal cardiac motor neurons, and the second is the impact of respiration on intrathoracic pressure. In the context of central networks, the concepts of functional connectivity and neurophysiological brain processes are significant. These include bottom-up and top-down processing, whereby information flows from higher brain regions to peripheral end organ systems or vice versa (Schulz et al., 2013a; Schulz et al., 2018).

In further steps, the new directionality indices will be applied and tested based on already analyzed and existing clinical data (cardiovascular system, cardiorespiratory system, central autonomic network) (Schulz et al., 2016; Schulz et al., 2018; Schulz et al., 2019) and compared with other methods in terms of their classification rate and accuracy. Moreover, further testing is necessary to ensure the robustness of the new directionality indexes and how time delays in the simulated models, as well as within the calculation procedure, influence the directionality results and the driver-response relationships.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

SS: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. AS: Writing – review & editing. K-JB: Writing – review & editing. JH: Writing – review & editing. GS: Writing – review & editing. AV: Writing – review & editing, Methodology, Supervision, Validation.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This study was supported by funding from the Software AG Stiftung (SE-P 14572).

Acknowledgments

I would also like to thank my former colleague and friend, Dipl.-Ing. Rico Schröder, who helped me with his expertise in solving problems that arose during the study. The successful, long-term collaboration with him has strongly influenced the content of this study.

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.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

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.

References

Adochiei, F., Schulz, S., Edu, I., Costin, H., and Voss, A. (2013). A new normalised short time PDC for dynamic coupling analyses. Biomed Tech (Berl) 58:1. doi: 10.1515/bmt-2013-4167

PubMed Abstract | Crossref Full Text | Google Scholar

Aguirre, R. R., Mustafa, M. Z., Dumenigo, A., Schulz, S., Voss, A., Goubran, B., et al. (2018). Influence of acute antipsychotic treatment on cardiorespiratory coupling and heart rate variability. Cureus 10:e2066. doi: 10.7759/cureus.2066

PubMed Abstract | Crossref Full Text | Google Scholar

Arnhold, J., Grassberger, P., Lehnertz, K., and Elger, C. E. (1999). A robust method for detecting interdependences: application to intracranially recorded EEG. Physica D: Nonlinear Phenomena 134, 419–430. doi: 10.1016/S0167-2789(99)00140-2

Crossref Full Text | Google Scholar

Baccala, L. A., and Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. Biol. Cybern. 84, 463–474. doi: 10.1007/PL00007990

PubMed Abstract | Crossref Full Text | Google Scholar

Bartsch, R. P., Liu, K. K., Bashan, A., and Ivanov, P. (2015). Network physiology: how organ systems dynamically interact. PLoS One 10:e0142143. doi: 10.1371/journal.pone.0142143

PubMed Abstract | Crossref Full Text | Google Scholar

Bashan, A., Bartsch, R. P., Kantelhardt, J. W., Havlin, S., and Ivanov, P. (2012). Network physiology reveals relations between network topology and physiological function. Nat. Commun. 3:702. doi: 10.1038/ncomms1705

PubMed Abstract | Crossref Full Text | Google Scholar

Baumert, M., Walther, T., Hopfe, J., Stepan, H., Faber, R., and Voss, A. (2002). Joint symbolic dynamic analysis of beat-to-beat interactions of heart rate and systolic blood pressure in normal pregnancy. Med. Biol. Eng. Comput. 40, 241–245. doi: 10.1007/BF02348131

PubMed Abstract | Crossref Full Text | Google Scholar

Beissner, F., Meissner, K., Bar, K. J., and Napadow, V. (2013). The autonomic brain: an activation likelihood estimation meta-analysis for central processing of autonomic function. J. Neurosci. 33, 10503–10511. doi: 10.1523/JNEUROSCI.1103-13.2013

PubMed Abstract | Crossref Full Text | Google Scholar

Candia-Rivera, D., Chavez, M., and De Vico Fallani, F. (2024). Measures of the coupling between fluctuating brain network organization and heartbeat dynamics. Netw Neurosci 8, 557–575. doi: 10.1162/netn_a_00369

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, W. G., Schloesser, D., Arensdorf, A. M., Simmons, J. M., Cui, C., Valentino, R., et al. (2021). The emerging science of Interoception: sensing, integrating, interpreting, and regulating signals within the self. Trends Neurosci. 44, 3–16. doi: 10.1016/j.tins.2020.10.007

PubMed Abstract | Crossref Full Text | Google Scholar

Craig, A. D. (2002). How do you feel? Interoception: the sense of the physiological condition of the body. Nat. Rev. Neurosci. 3, 655–666. doi: 10.1038/nrn894

PubMed Abstract | Crossref Full Text | Google Scholar

Dampney, R. A. (1994). Functional organization of central pathways regulating the cardiovascular system. Physiol. Rev. 74, 323–364. doi: 10.1152/physrev.1994.74.2.323

PubMed Abstract | Crossref Full Text | Google Scholar

Faes, L., Marinazzo, D., Jurysta, F., and Nollo, G. (2015). Linear and non-linear brain-heart and brain-brain interactions during sleep. Physiol. Meas. 36, 683–698. doi: 10.1088/0967-3334/36/4/683

PubMed Abstract | Crossref Full Text | Google Scholar

Ivanov, P. C. (2021). The new field of network physiology: building the human Physiolome. Front. Netw. Physiol. 1:711778. doi: 10.3389/fnetp.2021.711778

PubMed Abstract | Crossref Full Text | Google Scholar

Ivanov, P. C., Liu, K. K., and Bartsch, R. P. (2016). Focus on the emerging new fields of network physiology and network medicine. New J. Phys. 18:100201. doi: 10.1088/1367-2630/18/10/100201

PubMed Abstract | Crossref Full Text | Google Scholar

Kurths, J., Voss, A., Saparin, P., Witt, A., Kleiner, H. J., and Wessel, N. (1995). Quantitative analysis of heart rate variability. Chaos 5, 88–94. doi: 10.1063/1.166090

PubMed Abstract | Crossref Full Text | Google Scholar

Leem, C. H. (2016). Perspectives of physiome research. Integr. Med. Res. 5, 37–40. doi: 10.1016/j.imr.2015.12.004

PubMed Abstract | Crossref Full Text | Google Scholar

Lombardi, A., Guaragnella, C., Amoroso, N., Monaco, A., Fazio, L., Taurisano, P., et al. (2019). Modelling cognitive loads in schizophrenia by means of new functional dynamic indexes. NeuroImage 195, 150–164. doi: 10.1016/j.neuroimage.2019.03.055

PubMed Abstract | Crossref Full Text | Google Scholar

Marwan, N., Zou, Y., Wessel, N., Riedl, M., and Kurths, J. (2013). Estimating coupling directions in the cardiorespiratory system using recurrence properties. Philos. Trans. A Math. Phys. Eng. Sci. 371:20110624. doi: 10.1098/rsta.2011.0624

PubMed Abstract | Crossref Full Text | Google Scholar

McAllen, R. M. (1976). Proceedings: inhibition of the baroreceptor input to the medulla by stimulation of the hypothalamic defence area. J. Physiol. 257, 45P–46P

PubMed Abstract | Google Scholar

McCorry, L. K. (2007). Physiology of the autonomic nervous system. Am. J. Pharm. Educ. 71:78. doi: 10.5688/aj710478

PubMed Abstract | Crossref Full Text | Google Scholar

Montalto, A., Faes, L., and Marinazzo, D. (2014). MuTE: a MATLAB toolbox to compare established and novel estimators of the multivariate transfer entropy. PLoS One 9:e109462. doi: 10.1371/journal.pone.0109462

PubMed Abstract | Crossref Full Text | Google Scholar

Novak, V., Novak, P., De Champlain, J., Le Blanc, A. R., Martin, R., and Nadeau, R. (1993). Influence of respiration on heart rate and blood pressure fluctuations. J. Appl. Physiol. 74, 617–626. doi: 10.1152/jappl.1993.74.2.617

PubMed Abstract | Crossref Full Text | Google Scholar

Palus, M., Komarek, V., Hrncir, Z., and Sterbova, K. (2001). Synchronization as adjustment of information rates: detection from bivariate time series. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 63:046211. doi: 10.1103/PhysRevE.63.046211

PubMed Abstract | Crossref Full Text | Google Scholar

Palus, M., and Stefanovska, A. (2003). Direction of coupling from phases of interacting oscillators: an information-theoretic approach. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 67:055201. doi: 10.1103/PhysRevE.67.055201

PubMed Abstract | Crossref Full Text | Google Scholar

Quiroga, R. Q., Arnhold, J., and Grassberger, P. (2000). Learning driver-response relationships from synchronization patterns. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 61, 5142–5148. doi: 10.1103/PhysRevE.61.5142

PubMed Abstract | Crossref Full Text | Google Scholar

Rea, P. (2016). “Chapter 1 - overview of the nervous system” in Essential clinically applied anatomy of the peripheral nervous system in the head and neck. ed. P. Rea (London, United Kingdom: Academic Press), 1–20.

Google Scholar

Schiff, S. J., So, P., Chang, T., Burke, R. E., and Sauer, T. (1996). Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 54, 6708–6724. doi: 10.1103/PhysRevE.54.6708

PubMed Abstract | Crossref Full Text | Google Scholar

Schreiber, T. (2000). Measuring information transfer. Phys. Rev. Lett. 85, 461–464. doi: 10.1103/PhysRevLett.85.461

Crossref Full Text | Google Scholar

Schulz, S., Adochiei, F. C., Edu, I. R., Schroeder, R., Costin, H., Bar, K. J., et al. (2013a). Cardiovascular and cardiorespiratory coupling analyses: a review. Philos. Trans. A Math. Phys. Eng. Sci. 371:20120191. doi: 10.1098/rsta.2012.0191

Crossref Full Text | Google Scholar

Schulz, S., Bolz, M., Bär, K. J., and Voss, A. (2016). Central- and autonomic nervous system coupling in schizophrenia. Philos. Trans. A Math. Phys. Eng. Sci. 374:20150178. doi: 10.1098/rsta.2015.0178

PubMed Abstract | Crossref Full Text | Google Scholar

Schulz, S., Castro, M.R., Giraldo, B., Haueisen, J., and Voss, A. (2017a). "Multivariate high resolution joint symbolic dynamics (mHRJSD): a new tool to analyze couplings in physiological networks.", in: Biomedical engineering / Biomedizinische Technik).

Google Scholar

Schulz, S., Haueisen, J., Bar, K. J., and Andreas, V. (2015). High-resolution joint symbolic analysis to enhance classification of the cardiorespiratory system in patients with schizophrenia and their relatives. Philos trans a math Phys. Eng. Sci. 373:20140098. doi: 10.1098/rsta.2014.0098

PubMed Abstract | Crossref Full Text | Google Scholar

Schulz, S., Haueisen, J., Bär, K.J., and Voss, A. Quantification of cardiorespiratory coupling in acute schizophrenia applying high resolution joint symbolic dynamics. Computing in Cardiology Conference (CinC), (2013b), 101–104

Google Scholar

Schulz, S., Haueisen, J., Bar, K. J., and Voss, A. (2018). Multivariate assessment of the central-cardiorespiratory network structure in neuropathological disease. Physiol. Meas. 39:074004. doi: 10.1088/1361-6579/aace9b

PubMed Abstract | Crossref Full Text | Google Scholar

Schulz, S., Haueisen, J., Bär, K.-J., and Voss, A. (2019). Altered causal coupling pathways within the central-autonomic-network in patients suffering from schizophrenia. Entropy 21:733. doi: 10.3390/e21080733

PubMed Abstract | Crossref Full Text | Google Scholar

Schulz, S., Ricoy Castro, M., Giraldo, B., Haueisen, J., Bär, K.J., and Voss, A. (2017b). "Altered central cardiovascular network pattern in neuropathological disease - application of the three dimensional high resolution joint symbolic dynamics", in: Computing in Cardiology, 260–475.

Google Scholar

Schulz, S., Tupaika, N., Berger, S., Haueisen, J., Bär, K. J., and Voss, A. (2013c). Cardiovascular coupling analysis with high-resolution joint symbolic dynamics in patients suffering from acute schizophrenia. Physiol. Meas. 34, 883–901. doi: 10.1088/0967-3334/34/8/883

PubMed Abstract | Crossref Full Text | Google Scholar

Shoemaker, J. K., Norton, K. N., Baker, J., and Luchyshyn, T. (2015). Forebrain organization for autonomic cardiovascular control. Auton. Neurosci. 188, 5–9. doi: 10.1016/j.autneu.2014.10.022

PubMed Abstract | Crossref Full Text | Google Scholar

Taylor, E. W., Al-Ghamdi, M. S., Ihmied, I. H., Wang, T., and Abe, A. S. (2001). The neuranatomical basis of central control of cardiorespiratory interactions in vertebrates. Exp. Physiol. 86, 771–776. doi: 10.1111/j.1469-445X.2001.tb00043.x

PubMed Abstract | Crossref Full Text | Google Scholar

Voss, A., Kurths, J., Kleiner, H. J., Witt, A., Wessel, N., Saparin, P., et al. (1996). The application of methods of non-linear dynamics for the improved and predictive recognition of patients threatened by sudden cardiac death. Cardiovasc. Res. 31, 419–433. doi: 10.1016/S0008-6363(96)00008-9

PubMed Abstract | Crossref Full Text | Google Scholar

Wessel, N., Suhrbier, A., Riedl, M., Marwan, N., Malberg, H., Bretthauer, G., et al. (2009). Detection of time-delayed interactions in biosignals using symbolic coupling traces. EPL 87:10004. doi: 10.1209/0295-5075/87/10004

Crossref Full Text | Google Scholar

Wessel, N., Suhrbier, A., Riedl, M., Marwan, N., Malberg, H., Bretthauer, G., et al. (2011). Symbolic coupling traces for causality analysis of cardiovascular control. Annu Int Conf IEEE Eng Med Biol Soc 2011, 5935–5938. doi: 10.1109/IEMBS.2011.6091468

Crossref Full Text | Google Scholar

Ziegler, G., Dahnke, R., Yeragani, V. K., and Bar, K. J. (2009). The relation of ventromedial prefrontal cortex activity and heart rate fluctuations at rest. Eur. J. Neurosci. 30, 2205–2210. doi: 10.1111/j.1460-9568.2009.07008.x

PubMed Abstract | Crossref Full Text | Google Scholar