The geometry of synchronization: quantifying the coupling direction of physiological signals of stress between individuals using inter-system recurrence networks

In the study of synchronization dynamics between interacting systems, several techniques are available to estimate coupling strength and coupling direction. Currently, there is no general ‘best’ method that will perform well in most contexts. Inter-system recurrence networks (IRN) combine auto-recurrence and cross-recurrence matrices to create a graph that represents interacting networks. The method is appealing because it is based on cross-recurrence quantification analysis, a well-developed method for studying synchronization between 2 systems, which can be expanded in the IRN framework to include N > 2 interacting networks. In this study we examine whether IRN can be used to analyze coupling dynamics between physiological variables (acceleration, blood volume pressure, electrodermal activity, heart rate and skin temperature) observed in a client in residential care with severe to profound intellectual disabilities (SPID) and their professional caregiver. Based on the cross-clustering coefficients of the IRN conclusions about the coupling direction (client or caregiver drives the interaction) can be drawn, however, deciding between bi-directional coupling or no coupling remains a challenge. Constructing the full IRN, based on the multivariate time series of five coupled processes, reveals the existence of potential feedback loops. Further study is needed to be able to determine dynamics of coupling between the different layers.

1 Introduction

In studies of the synchronization dynamics of interacting systems, many techniques are available to determine whether there is a (statistical) dependence between observed time series of the variables of interest, such as the cross-correlation function or the mutual information measure. A more complicated research question concerns causal inference from observed time series data, in which the goal is to quantify and qualify the nature of the coupling dynamics between the interacting systems in terms of strength and direction (unidirectional, symmetric, asymmetric bidirectional, etc.). Many different techniques are available like linear and nonlinear Granger causality (Granger, 1969; Marinazzo et al., 2011), various information theoretic measures such as transfer entropy (Schreiber, 2000), conditional mutual information (Paluš et al., 2001) and synergistic information (Quax et al., 2017), Cross-RQA (Zbilut et al., 1998; Shockley et al., 2002; Richardson et al., 2007) and Cross-Mapping techniques (Sugihara et al., 2012; Ma et al., 2014; Ma et al., 2017; Leng et al., 2020). The accurate identification of the coupling direction, however, remains a challenge (Feldhoff et al., 2012; Runge et al., 2019). Determining which technique is ‘best’ may depend on the specific context, because different techniques come with different sets of assumptions about the data generating process (e.g., linear, stationary, etc.) and different requirements for the observed data (e.g., time series length, noise levels). In addition, some measures of coupling direction have been shown to give inaccurate results based on the coupling strength and parameters such as frequency and phase differences between the signals (Feldhoff et al., 2012).

In the cognitive, behavioral and neurosciences, Cross-Recurrence Quantification analysis (CRQA) has frequently been used to study coupling dynamics of human behavior and physiology, for example, in postural sway, heart rate variability, interpersonal behavior and emotion dynamics (Shockley, 2005; Konvalinka et al., 2011; Cox et al., 2016; Main et al., 2016; Vink et al., 2018). In the present paper we demonstrate how measures obtained from inter-system recurrence networks (Feldhoff et al., 2012) can be used to analyze the synchronization between physiological signals that were simultaneously observed in two interacting individuals. Inter system recurrence networks are an extension of multilayer recurrence networks, in which the cross-recurrence matrix is used to quantify the interdependence between different network layers (Zou et al., 2019). The goal is not to present an optimal method for detecting coupling direction, but to explore the applicability of inter-system recurrence networks to analyze empirical data observed in a healthcare setting. One reason to choose the (multiplex) recurrence network framework is the potential for application in clinical practice (Hasselman, 2022). Due to affordable wearable sensors and electronic devices that allow experience sampling and momentary ecological assessment (self-reports of physiological, psychological, or emotional states), many healthcare institutes are looking for ways to integrate such process monitoring data into diagnostic and intervention (Fartacek et al., 2016; Olthof et al., 2019; Hasselman, 2022), and explore opportunities to personalize such procedures (Schiepek et al., 2016; Olthof et al., 2023).

The data analyzed in the present paper concern a client in residential care with severe to profound intellectual disabilities (SPID) and their professional caregiver, both equipped with wearable sensors which recorded a number of physiological signals (acceleration, blood volume pressure, electrodermal activity, heart rate and skin temperature) during regular daily shifts (Simons et al., 2021). A substantial percentage of people with SPID show severely challenging behavior, which most frequently includes aggression and self-injury (Janssen et al., 2002). In residential care, incidents in which SPID clients display challenging behavior are obviously a stressor for the client, but also for their environment, in particular the professional caregiver tasked with resolving as well as preventing such incidents. The frequency of occurrence and the severity of such incidents is known to affect caregivers’ wellbeing and consequently the quality of the care they are able to provide (Devereux et al., 2009). The focus of the analysis will be to assess the coupling strength and direction of physiological signals in the period before an incident with challenging behavior occurred (see Figure 4 for an example of the data). The goal is to gain insight in the role of physiological coupling dynamics between individuals leading up to an incident which may inform methods to reduce or prevent the frequency of occurrence and severity of such incidents. First, we formally introduce the concept of inter-system networks (IRN) and the measures that can be used to quantify coupling dynamics, second, we use simulated data to demonstrate the behavior of IRN measures in different coupling scenarios, finally we present an analysis of the empirical client-caregiver data using the same measures.

2 Inter-system recurrence networks (IRN)

Recurrence networks represent a graph-theoretical approach to time series analysis that is based on recurrence analysis, which evaluates the proximity of phase space trajectories. The structural properties of the recurrence network are associated to the geometry of the underlying attractor (Zou et al., 2019). In this section we describe how to construct a so-called inter-system recurrence network which is a representation of 2 (or more) interacting recurrence networks and examines to what extent the different systems can be understood to share the same state space or show evidence of coupled dynamics.

2.1 Constructing a recurrence network

An $M$ dimensional multivariate timeseries ${\left\{y_i\right\}}_{i=1}^{N_y}$ with $y_i=y\left(t_i\right)$ that captures the relevant degrees of freedom of a complex system $Y$ can be interpreted to represent its state evolution in an $M$ dimensional state space. Each state is represented by an m-tuple of coordinates, the state vector $\overrightarrow{y_i}$ describes a trajectory through the $M$ dimensional state space. To study the dynamics of system $Y$, first, a distance matrix is created that encodes the closeness of each observed state relative to all other states. Second, the selection of a distance threshold $\epsilon$ determines which states should be considered close, or recurring states. Applying the threshold to the distance matrix results in a binary auto-recurrence matrix

${\mathit{R}}_{ij}\left(\epsilon \right)={\mathit{R}}^Y\left(y_i,y_j|\epsilon \right)=\mathrm{\Theta }\left(\epsilon -‖\overrightarrow{y_i}-\overrightarrow{y_j}‖\right),i,j=1,\dots ,N(1)$

in which $‖\bullet ‖$ is a distance norm (e.g., Euclidean, Chebyshev) and $\mathrm{\Theta }\left(\bullet \right)$ the Heaviside function, returning a 1 if a distance value is below the threshold $\epsilon$, and 0 otherwise. The threshold $\epsilon$ determines how many recurrent states are found and is therefore monotonically related to the Recurrence Rate (RR), the number of $\epsilon$-recurrent points on the maximal possible number of points. Often, a predetermined recurrence rate is used to select a threshold $\epsilon$ (e.g., 1%–5%). The recurrence matrix is used in Recurrence Quantification Analysis (RQA) to obtain estimates of dynamical invariants and complexity measures from the (distribution of) $\epsilon -$ recurrent points and the vertical and diagonal line structures they form (Marwan et al., 2007).

The recurrence matrix can be used to generate a graph, a so-called recurrence network (RN) (Zou et al., 2019). In a recurrence network, each vertex (node) represents a time point and each edge (a connection between 2 vertices) indicates that the values observed at these time points are recurring values, given the threshold $\epsilon$. In general, to create an $\epsilon -RN$ of system $X$, the recurrence matrix $R_{ij}\left(\epsilon \right)=R^X\left({\epsilon }^X\right)$, with its main diagonal set to 0, is considered to be the adjacency matrix $A_{ij}\left(\epsilon \right)=A^X\left({\epsilon }^X\right)$ of a complex network (Zou et al., 2019). This is an unweighted, undirected graph, in which each vertex $V_i$ represents a state of $X$ observed at time point $t_i$. For example, the graph $G^X$ in Figure 1C, represents the $\epsilon$ recurring values of time series X in Figure 1A. That is, the values $‖\overrightarrow{X_i}-\overrightarrow{X_j}‖$ smaller than ${\epsilon }^X$ (the black circle in Figure 1A) are considered recurring values of X (Eq. 1). The top left part of the matrix in Figure 1B (blue cells) is the recurrence matrix $R^X\left({\epsilon }^X\right)$, turned into the adjacency matrix $A^X\left({\epsilon }^X\right)$, the rows and columns of the matrix represent the vertices of the network, a cell containing a 1 indicates the vertices are connected by an edge. In Figure 1C, vertices v2 and v3 are connected by an edge $E_{i\to j}$ with $i=2$ and $j=2$, indicating that the state at $V_{i=2}$ recurs at $V_{j=3}$. The construction of graph $G^Y$ in Figure 1C follows the same steps but is based on timeseries Y in Figure 1A.

FIGURE 1

FIGURE 1. (A) shows two time series X (blue) and Y (green) and thresholds ε^X = ε^Y (black circle) and ε^XY = ε^YX (grey circle). (B) shows the composite adjacency matrix built from the recurrence matrices of X (top left, blue cells) and Y (bottom right, green cells) and their cross-recurrence matrices (bottom left and top right, grey cells). (C) shows two recurrence networks G^X and G^Y, joined together through the cross-recurrence matrix as an interconnected recurrence network. To calculate triangle-based measures like the clustering coefficient, intersystem as well as intra system edges must be evaluated, for example, the triangle starting in v1: E(G^XY) = {(v1, u1), (u1, u2), (u2, v1)}, see the bold, underlined 1s in (B). Here, the cross-edge density 𝜌^XY = 𝜌^YX = 0.33, cross-clustering C^XY = 0.33 ≠ C^YX = 0.

2.2 Constructing an inter-system recurrence network

An extension of RQA developed to study synchronization and coupling dynamics is called Cross-RQA (Marwan and Kurths, 2002; Shockley et al., 2002) and involves evaluating whether the state evolution of two systems $X$ and $Y$ can be considered close, or recurring trajectories in a shared state space. The procedure to construct a cross-recurrence matrix is similar to the auto-recurrence matrix, but instead of calculating the distances between the states of $Y$ at different time points, the distances between the states of $X$ relative to states of $Y$ are evaluated as

${\mathit{C}\mathit{R}}_{ij}\left({\epsilon }^{XY}\right)={\mathit{C}\mathit{R}}^{XY}\left(x_i,y_j|{\epsilon }^{XY}\right)=\mathrm{\Theta }\left({\epsilon }^{XY}-‖\overrightarrow{x_i}-\overrightarrow{y_j}‖\right),i,j=1,\dots ,N(2)$

where $‖\bullet ‖$ is again a distance norm (e.g., Euclidean, Chebyshev) and $\mathrm{\Theta }\left(\bullet \right)$ the Heaviside function. The cross-recurrence matrix is generally not symmetrical around its diagonal, in fact, the asymmetry provides information about the coupling dynamics and can be visualized using Diagonal Cross-Recurrence Profiles (DCRP) (Marwan and Kurths, 2002; Wallot and Leonardi, 2018), also known as the $\tau$-recurrence rate (Marwan et al., 2007). Figure 2 shows Cross-Recurrence Plots (a visualization of the CR) and DCRPs for 2 coupled oscillators in 4 different coupling scenarios (see Section 3 for more details). There have been attempts to use statistical methods to study the DCRP patterns at a group level (Main et al., 2016), however, often DCRPs require visual inspection in order to make a reliable interpretation of the coupling strength and direction about individual cases. Several measures calculated from so-called inter-system recurrence networks (IRN) have been suggested that do not require such visual inspection (Donges et al., 2011; Feldhoff et al., 2012). Inter-system recurrence network can be constructed from the recurrence and cross-recurrence matrices of system $X$ and $Y$. The idea is to use Cross-RQA to determine whether there should be cross-network edges, such as the grey edges between $G^X$ and $G^Y$ in Figure 1C. More formally, to study the coupling dynamics between $G^X$ and $G^Y$ we can use the cross-recurrence matrices ${CR}^{XY}\left({\epsilon }^{XY}\right)$ and ${CR}^{YX}={\left[{CR}^{XY}\left({\epsilon }^{XY}\right)\right]}^T$, to build the inter-system recurrence matrix

$\mathbf{I}\mathbf{R}\left(\epsilon \right)=\left(\begin{array}{cc}{\mathbf{R}}^X\left({\epsilon }^X\right)& {\mathbf{CR}}^{XY}\left({\epsilon }^{XY}\right)\\ {\left[{\mathbf{CR}}^{XY}\left({\epsilon }^{XY}\right)\right]}^T& {\mathbf{R}}^Y\left({\epsilon }^Y\right)\end{array}\right)(3)$

with $\epsilon$ representing three threshold values ${{\epsilon }^X,{\epsilon }^Y,\epsilon }^{XY}$ that can be dissimilar. It is recommended to have less cross-system recurrences than intra-system recurrences in the network, which can be controlled by choosing a threshold based on the recurrence rate (RR) it yields $\left.{RR}^{XY}\le {\left(RR\right.}^X\approx {RR}^Y\right)$ (Feldhoff et al., 2012). $\text{IR}\left(\epsilon \right)$ is a symmetric matrix, which can be turned into the adjacency matrix

$\mathrm{A}\left(\epsilon \right)=IR\left(\epsilon \right)-{\mathbb{I}}_N(4)$

where ${\mathbb{I}}_N$ is the identity matrix with $N=N_X+N_Y$. This composite matrix is shown in Figure 1 Panel B. The inter-system network (IRN) is an undirected, unweighted, simple graph; however, it consists of two unipartite graphs (the recurrence networks of system $X$ and $Y$) and two bipartite graphs (the cross-recurrence networks for system $X\to Y$ and $Y\to X$).

FIGURE 2

FIGURE 2. (A) shows simulation of 201 time steps of 2 damped oscillators X and Y (standardized) in anti-phase for 4 different types of coupling scenarios: Uncoupled, X driving Y, Y driving X, and Competitive bi-directional. The panels in row (B) show the Cross Recurrence Plots (CRPs) based on the joint 2D state space of X (X, dX) and Y (Y, dY). Row (C) shows the Diagonal Cross-Recurrence Profiles (DCRP), representing the distribution of recurrent points around the main diagonal (0 on the x-axis of the DCRP). The grey lines are based on 39 surrogate DCRP analyses in which the temporal order of Y was randomized. If the (black) line falls above or below the grey line the observed value is above/below what may be expected based on chance (p < .05).

2.3 IRN measures

The first measure of interest for an IRN is the local vertex cross-degree:

$k_v^{XY}={\displaystyle \sum _{u\in V_Y}}{A}_{vu}(5)$

with $v\in V_X$, representing the number of edges that connect a vertex in $G^X$ to any vertex in $G^Y$. This is comparable to distinguishing between the in- and out-degree in a directed network, but in the present case the purpose is to distinguish between within network and between-network connections. In Figure 1C, the cross-degrees for vertices $\left\{v_1,v_2,v_3\right\}$ of $G^X$ are 2,1,0, respectively, which is the same as for vertices $\left\{u_1,u_2,u_3\right\}$ of $G^Y$.

2.3.1 The cross-edge density

${\rho }^{XY}=\frac{1}{N_X{N}_Y}{\displaystyle \sum _{v\in V_X,u\in V_Y}}{A}_{vu}(6)$

is a proportion representing the total number of cross-edges over the maximum possible number of edges which is equal to the cross-recurrence rate ${RR}^{XY}$. In Figure 1C the cross-edge density is $3/9=0.33$.

The local cross-clustering coefficient is an estimate of the probability that two neighbors of vertex $v\in V_X$ randomly drawn from graph $G^Y$ are also neighbors

$C_v^{XY}=\frac{1}{k_v^{YX}\left(k_v^{YX}-1\right)}{\displaystyle \sum _{u_i,u_j\in V_Y}}{A}_{v{u}_i}{A}_{u_i{u}_j}{A}_{u_{j}v},i,j=1,\dots ,N(7)$

with $i\ne j$ and for vertices with a cross-degree of 0 or 1, $C_v^{XY}=0$. This clearly distinguishes the cross-clustering coefficients from Cross-RQA measures, because their calculation requires the cross-recurrence matrix as well as the auto-recurrence matrices. In Figure 1C the local clustering coefficients for $\left\{v_1,v_2,v_3\right\}$ are 1,0,0, but they are all 0 for $\left\{u_1,u_2,u_3\right\}$. The global cross-clustering coefficient is simply the average over all vertices

$C^{XY}={〈C_v^{XY}〉}_{v\in V_X}(8)$

which gives $C^{XY}=0.33\ne C^{YX}=0$. If the distribution of edges within $G^Y$ and between $G^X$ and $G^Y$ is uncorrelated, the expected cross-clustering coefficient would be $C^{XY}\approx {\rho }^Y\approx {RR}^Y$, that is the cross-edges would be randomly and independently distributed between the graphs (Feldhoff et al., 2012). Coupling constants $C^{XY}\ll {\rho }^Y$ and $C^{XY}\gg {\rho }^Y$ indicate substantial (anti-)persistent correlations exist in the connectivity structure between the networks.

The cross-clustering coefficient measures provide information about the way the two systems are coupled and can be used to infer a coupling direction. In Section 3 we demonstrate the behavior of these measures for different coupling scenarios.

3 Simulation study: coupled oscillators

3.1 Method

The cross-clustering coefficients that can be computed from an IRN have been shown to take on different values for different coupling scenarios (Donges et al., 2011; Feldhoff et al., 2012). Predictions can be made about the expected relative magnitudes of the measures, here we demonstrate the behavior of these measures using two coupled damped oscillators, because our goal is to analyze physiological signals, which are usually periodic and may vary in amplitude and frequency over time. However, the expectations presented here do not differ from previous studies (Donges et al., 2011; Feldhoff et al., 2012).

If the oscillators are uncoupled, any cross-recurrences found will be spurious, however, they are expected because of the similarity between the oscillators, the only difference will be a phase lag. Therefore, the expectation is that ${\mathit{C}}^{\mathit{X}\mathit{Y}}\approx {\mathit{C}}^{\mathit{Y}\mathit{X}}$. If there would be a unidirectional coupling $\mathit{X}\to \mathit{Y}$ of sufficient strength, it is likely that for two vertices ${\mathit{v}}_{\mathit{i}},{\mathit{v}}_{\mathit{j}}\in {\mathit{V}}_{\mathit{X}}$ there exists a vertex ${\mathit{u}}_{\mathit{k}}\in {\mathit{V}}_{\mathit{Y}}$ that is connected to both ${\mathit{v}}_{\mathit{i}}$ and ${\mathit{v}}_{\mathit{j}}$. This would form a closed triplet from the perspective of ${\mathit{u}}_{\mathit{k}}$, which is a vertex of the driven system $\mathit{Y}$. Therefore, the expectation would be to see ${\mathit{C}}^{\mathit{Y}\mathit{X}}>{\mathit{C}}^{\mathit{X}\mathit{Y}}$ for a unidirectional coupling $\mathit{X}\to \mathit{Y}$. If the direction would be reversed and $\mathit{X}$ would be the driven system $\mathit{Y}\to \mathit{X}$, the opposite relation is expected ${\mathit{C}}^{\mathit{X}\mathit{Y}}>{\mathit{C}}^{\mathit{Y}\mathit{X}}$.

The analysis scripts are available at https://osf.io/5n9kv/.

3.1.1 Procedure

To study the behavior of inter-system recurrence network measures, two linear damped oscillators X and Y were simulated:

$\frac{d^{\mathbf{2}}X}{d^{\mathbf{2}}t}={-\mathrm{\eta }}_x\mathrm{X}+{\mathrm{\zeta }}_x\frac{dX}{dt}+{\mathrm{\gamma }}_{yx}\mathrm{Y}(9)$

$\frac{d^{\mathbf{2}}Y}{d^{\mathbf{2}}t}={-\mathrm{\eta }}_y\mathrm{Y}+{\mathrm{\zeta }}_y\frac{dY}{dt}+{\mathrm{\gamma }}_{xy}\mathrm{X}$

Parameters ${\mathrm{\eta }}_x$, ${\mathrm{\eta }}_y$ represent the frequency and ${\mathrm{\zeta }}_x$, ${\mathrm{\zeta }}_y$ the damping factor of X and Y respectively. The coupling dynamics are governed by the parameters ${\mathrm{\gamma }}_{xy}$ and ${\mathrm{\gamma }}_{yx}$, representing the effect of X on $d^2\mathrm{Y}$ and Y on $d^2\mathrm{X}$ respectively. For all simulations the initial values were $X_0=3$ and $Y_0=-3$, resulting in an antiphase pattern. The frequency parameters were ${\mathrm{\eta }}_x={\mathrm{\eta }}_y=0.3$ and the damping was set to ${\mathrm{\zeta }}_x={\mathrm{\zeta }}_y=-0.05$. The 4 coupling scenarios Uncoupled, X driving Y, Y driving X, and Bi-directional (competitive) were realized by changing the coupling parameters: ${\mathrm{\gamma }}_{xy}={\mathrm{\gamma }}_{yx}$ = 0 for the Uncoupled scenario; ${\mathrm{\gamma }}_{xy}=+0.05,{\mathrm{\gamma }}_{yx}=0.00$ for the X drives Y scenario; ${\mathrm{\gamma }}_{xy}=0.00,{\mathrm{\gamma }}_{yx}=-0.05$ for the Y drives X scenario and ${\mathrm{\gamma }}_{xy}=+0.05,{\mathrm{\gamma }}_{yx}=-0.05$ for the Bi-directional scenario. Timeseries of length 201 were simulated using package deSolve (Soetaert et al., 2010) for the R computing environment (R Core Team, 2022), with default ODE solver lsoda (Petzold, 1983).

To analyze the data, auto-recurrence matrices for $R^X\left({\epsilon }^X\right)$ and $R^Y\left({\epsilon }^Y\right)$ and cross-recurrence matrix ${CR}^{XY}\left({\epsilon }^{XY}\right)$ were generated based on the (standardized) state vectors $\left[X,dX\right]$ and $\left[Y,dY\right]$ which describe the 2D state space of the damped oscillators (i.e., dimension = 2, but no delay embedding procedure was required to reconstruct phase space). Threshold values ${\epsilon }^X$ and ${\epsilon }^Y$ were chosen to yield ${RR}^X={RR}^Y=.05$, and ${\epsilon }^{XY}$ to get ${RR}^{XY}=.03$. In addition, DCRPs were created to study the distribution of recurrent points around the main diagonal and evaluate how the interaction patterns are associated to the IRN measures.

3.2 Results

Figure 2 shows the cross-recurrence plots (CRP) in row B and the DCRPs in row C. Recurrent points in the diagonals of the upper triangle of the CRP are due to values occurring first in X (negative numbers on the x-axis of the DCRP), recurrent points in the lower triangle (positive diagonal numbers) are due to values occurring first in Y. If a peak is centered on the main diagonal this would suggest the systems are fully synchronized and display the same states at approximately the same point in time. If a peak is centered around diagonals in the lower or upper triangle of the matrix the inference can be made either system $X$ or $Y$ is leading the dynamics because states that occurred at an earlier point in time in $X$ recur at a later point in time in $Y$, and vice versa. This can be seen in Figure 2, row C for the scenarios $X\to Y$ and $Y\to X$. If two or more off-diagonal peaks are observed, to the left and right of the main diagonal, this can point to a bi-directional coupling. In Figure 2 such a pattern can be seen for both the uncoupled and bi-directional scenarios. For the uncoupled scenario the symmetrical pattern emerges because the oscillators have the same frequency and damping parameters, there is however no peak around the main diagonal. For the bi-directional coupling there is a peak at the lag which represents the phase difference, and the pattern is symmetrical around that diagonal.

The interpretation by visual inspection of the DCRPs is facilitated by our knowledge of the system under study. With an unknown system, the unidirectional coupling would not pose a problem in this case, but distinguishing between the uncoupled scenario and the bi-directional coupling would likely be uncertain. Figure 3C shows the IRN for the scenario “Y drives X” as two coupled recurrence networks (Figures 3A, B) in a spiral layout (Hasselman and Bosman, 2020).

FIGURE 3

FIGURE 3. (A,B) represent the recurrence networks of oscillator X and Y for the coupling scenario “Y drives X.” In the spiral graph layout, the temporal order is retained, with t = 1 at the lower end of the spiral and t = 201 at the upper end. The vertices are colored to distinguish 4 epochs of equal length. If edges connect vertices from the same epoch, the edge color is the same as the vertex color. The size of the vertices scales with the vertex degree of each node. (C) shows the inter-system recurrence network. The grey lines connecting the spirals represent the cross-edges.

Table 1 shows the results for the cross-clustering coefficients, as well as the different thresholds used to generate the (cross-)recurrence matrices. The magnitude differences are all in the expected direction, that is, the difference between the uncoupled and bi-directional scenarios are $\approx 0$, when $X\to Y$, $C^{YX}$ and are larger than $C^{XY}$, which is reversed for $Y\to X$. These results correspond to the interpretation based on the visual inspection of the DCRPs in Figure 2. In this stylized example, the uncoupled scenario yields exactly the same cross-clustering for the different directions. This can be used to distinguish it from the bidirectional scenario, which gives a difference close to zero. This will unlikely be possible with noisy real-world empirical data. The bi-directional clustering coefficients are similar to the largest coefficients obtained for the driven scenarios, which all represent the same coupling strength magnitude in the simulation (${\mathrm{\gamma }}_{xy}=+0.05,{\mathrm{\gamma }}_{yx}=-0.05$).

TABLE 1

TABLE 1. Inter-system recurrence network threshold parameters and clustering measures.

To study the relationship between coupling strength and the magnitude of the cross-clustering coefficient, further simulations were conducted in which ${\mathrm{\gamma }}_{xy}$ and ${\mathrm{\gamma }}_{yx}$ were varied from 0 to 0.1 in 40 steps. For each combination of ${\mathrm{\gamma }}_{xy}$ and ${\mathrm{\gamma }}_{yx}$ the difference ${∆C=C}^{XY}-C^{YX}$ was calculated, resulting in a 41 by 41 matrix of differences.

Figure 4 is a graphical representation of the 41 by 41 matrix. It shows that if both coupling strengths become very large, the two signals become very similar and as expected, the cross-clustering coefficients do not distinguish anymore between the true coupling differences (i.e., $C^{XY}-C^{YX}\approx 0$). Also, there are specific combinations in which the coupling strength difference appears to indicate the opposite direction of what is expected, this occurs around the transition where the coupling strengths become stronger than the damping factor and the oscillations exponentially increase instead of decrease. The inset panels labelled A, B and C in Figure 4 show this in some more detail. Panel A and C represent the expected relationship (Y leads X) based on the coupling strengths, but the value of ∆C indicates a reversal (X leads Y) in panel B. This occurs due to the initial phase difference, when the parameter settings produce oscillator dynamics with approximately equal amplitudes that decay very slowly.

FIGURE 4

FIGURE 4. Simulation results for the 41 × 41 matrix of different coupling strength values. Color coding represents the difference ${∆\mathit{C}=\mathit{C}}^{\mathit{X}\mathit{Y}}-{\mathit{C}}^{\mathit{Y}\mathit{X}}$. Red represents ${\mathit{C}}^{\mathit{X}\mathit{Y}}<{\mathit{C}}^{\mathit{Y}\mathit{X}}$, blue ${\mathit{C}}^{\mathit{X}\mathit{Y}}>{\mathit{C}}^{\mathit{Y}\mathit{X}}$ and white ${\mathit{C}}^{\mathit{X}\mathit{Y}}\approx {\mathit{C}}^{\mathit{Y}\mathit{X}}$. Three examples of oscillator dynamics are shown in the panels on the right of the matrix. The expected relationship (Y leads X) can be seen in (A,C), but there is a reversal in (B). See text for details.

4 Coupling between physiological processes

4.1 Method

Data were collected in the context of a larger study of the relation between physiological stress and challenging behavior in clients in residential care with severe to profound intellectual disabilities (SPID) (Simons et al., 2021). The goal of the study was to examine whether information about the dynamics of physiological stress in clients with SPID and their professional caregivers can be used to understand the emergence of challenging behavior and potentially prevent incidents from occurring. In the present study we examine whether a coupling direction can be detected in physiological signals between client and caregiver in the period before an incident occurred.

4.1.1 Participants, materials and procedure

The study was approved by the Faculty Ethics Committee of the Radboud University (NL-number: NL71683.091.19) 12 February 2020. Informed consent was obtained from caregivers and parents or legal representatives of clients. We analyzed 33 incidents experienced by one client-caregiver dyad over the course of 12 days between 3 August 2022 and 27 September 2022 during regular shifts of the caregiver. Both participants wore an Empatica E4 wristband (Milstein and Gordon, 2020; Schuurmans et al., 2020) which records acceleration (ACC), blood volume pressure (BVP), electrodermal activity (EDA), heart rate variability (HR) and temperature (TEMP). Whenever an incident with challenging behavior occurred the time was logged by pressing the event marker button on the wristband. After the incident, the time was also logged manually by the caregiver, providing more detailed information about the incident.

4.1.2 Data analysis

The raw data from the Empatica E4 wristband were imported using R package wearables (de Looff et al., 2022). Figure 5 is an example of the data centered on the time at which an incident occurred; the analyses are based on the 35 min before the incident. There were 4 incidents for which 2 or more timeseries were shorter than 25 min. The most prominent cause was that the incident occurred within 25 min of putting on the wristband, other cases were due to sensor failure. These timeseries were removed, leaving 29 incidents with 5 timeseries. All timeseries were resampled to the same frequency (1 Hz), yielding time series of length 2,100. To construct the IRN, two auto-recurrence matrices were created (one for the client and one for the caregiver) and one cross-recurrence matrix for each incident and each variable. The matrices are based on state space vectors constructed from delay embedded time series. For each incident and variable separately, the embedding lag and embedding dimension were determined for client and caregiver using the mutual information criterion to determine the optimal embedding lag and false nearest neighbor analysis to determine the embedding dimension (Marwan et al., 2007). One set of parameters was chosen for all embeddings (embedding lag = 100 and embedding dimension = 5). The recurrence threshold for the auto-recurrence matrices was chosen to get 5% recurrent points and the cross-recurrence matrix to get 4% recurrent points. Subsequently, the local and global cross-clustering coefficients were calculated.

FIGURE 5

FIGURE 5. Raw time series data from caregiver and patient 30 min before and after an incident with challenging behavior occurred (dotted line at 17:55). Shown are acceleration (ACC), blood volume pressure (BVP), electrodermal activity (EDA), heart rate (HR) and temperature (TEMP).

4.2 Results

Figure 6 displays the global cross-clustering coefficients for the 29 incidents, for the client (X) and the caregiver (Y). The location of the dots above or below the line $C^{XY}=C^{YX}$ is indicative of the coupling direction $X\to Y\left(C^{XY}<C^{YX}\right)$, and $Y\to X\left(C^{XY}>C^{YX}\right)$. This means the upper triangle shows cases in which the client “drives” the caregiver’s physiology, whereas the lower triangle shows cases in which the caregiver “drives” the client’s physiology. Points close to, or on the diagonal can represent uncoupled processes or bi-directional coupling dynamics.

FIGURE 6

FIGURE 6. Global cross-clustering coefficients for 29 incidents for 5 different inter-system networks based on different physiological time series observed in a client (X) and a caregiver (Y). The inset shows cross-clustering coefficients in range [0,0.05].

The points in the magnified area represent cases in which $C^{XY}\ll {\rho }^Y$ and $C^{YX}\ll {\rho }^X$, of which it is expected we should find signals with anti-persistent correlations. The results agree with these expectations as variables ACC and BVP predominantly yield values smaller than .05, whereas the coefficients for HR, EDA and TEMP are mostly larger than .05.

Table 2 shows the frequencies with which client and caregiver were driving the interaction, or whether the interaction was bi-directional or uncoupled, for each variable. If the difference $\left|C^{XY}-C^{YX}\right|\le .01$, this was counted as “Bi-directional/Uncoupled,” $C^{XY}-C^{YX}>.01$ as “Caregiver leading” and $C^{XY}-C^{YX}<-.01$ was counted as “Client leading.”

TABLE 2

TABLE 2. Frequency of coupling directions per variable across all 29 incidents (percentage).

Across all 29 incidents the three types of coupling dynamics occur with equal frequency (see the row Total in Table 2), however, there are differences between the variables. Most notable: Compared to the caregiver the client leads twice as often for the variable ACC, for BVP the bi-directional/uncoupled type was clearly dominant, and for EDA the caregiver most often was leading the interaction.

Finally, due to the choice for one sampling frequency and one set of embedding parameters for all variables, it is possible to calculate a cross-recurrence matrix for all variable pairs, so a full IRN representing all the observed variables can be constructed. The adjacency matrix has the auto-recurrence matrices of the 5 variables on its main diagonal, all off-diagonal cells represent cross-recurrence matrices between different variable pairs. Figure 7 is a representation of such a multiplex inter-system recurrence network for 1 incident. The edge weights represent the absolute clustering coefficient difference and have been colored to reflect the implied coupling direction. The size of a vertex represents its weighted degree, the vertex strength (Barrat et al., 2004). For this incident an interesting path can be seen that connects the vertices with the greatest strength: The edges connecting EDA–HR–TEMP, are characterized by the client driving the dynamics. A path connecting the second largest set of 3 vertices in terms of strength, is ACC–TEMP–HR. In this loop, the dynamics between ACC–TEMP and ACC–HR are driven by the caregiver, as mentioned before, TEMP–HR is driven by the client.

FIGURE 7

FIGURE 7. The complete multiplex inter-system recurrence network for 1 incident. Edges represent the differences between the cross-clustering coefficients for each variable. The colors are labelled according to the implied coupling direction. The vertex size reflects their strength (weighted degree).

5 General discussion and conclusion

The purpose of the present paper was to examine whether inter-system recurrence networks can be used to determine the coupling direction between physiological signals associated with stress measured simultaneously in a client-caregiver dyad, in which challenging behavior by the client can lead to incidents that are stressful for everyone involved. The question is whether such analyses can inform potential intervention or prevention efforts intended to reduce the number of incidents.

The simulation study in Section 3 provides results that are similar to previous studies (Feldhoff et al., 2012), showing that the magnitudes of the cross-coupling coefficients are associated to coupling strength and can be used to determine a coupling direction. In addition, we show that the qualification of the patterns of DCRPs agree with the inferred coupling direction. However, analogous to previous studies (Feldhoff et al., 2012), we also find parameter settings that may lead to uncertain or incorrect inferences about the coupling direction, given the true coupling strengths. In this particular case, the question is whether this should be considered an error of the method. If we were presented with the time series that yielded this result (e.g., the inset in Figure 4), absent any knowledge of the data generating process, we would conclude the estimated coupling direction to be accurate. It is indeed the case that the oscillator indicated by the cross-clustering coefficients is driving the interaction. The conclusion should be that these analyses do not recover the “true” coupling parameters, if such a thing is even possible, but instead estimate the observed, or implied coupling direction. An open question to be resolved is how to deal with coefficient differences close to 0. These could imply bidirectional, as well as completely uncoupled dynamics. One could argue that completely uncoupled, independent dynamics are unlikely in a measurement context in which individuals are in fact interacting, however, we cannot be completely certain based exclusively on the estimate.

The analysis of the physiological data related to 29 incidents with challenging behavior (Section 4) yielded interpretable results, as well as some corroboration of theoretical constraints on the measures described in Section 2. To start with the latter, Figure 6 clearly shows the expected relationship between (anti-)persistent correlations and edge density. The cross-clustering coefficients for the anti-persistent signal BVP are all equal or smaller than the recurrence rate of .05. For variable ACC there are 8 points outside the magnified region, but in those cases one of the two cross-clustering coefficients is always smaller than .05.

Based on Table 2 one can conclude that this dyad does not have a clear driver of the (physiological) interaction dynamics in the 35 min before an incident occurs. There are some differences between the individual variables, but Figure 7 shows that such differences at the aggregate level may be less informative to understand the interaction dynamics preceding a particular incident. The figure reveals that when all possible cross-edges between variables are considered, feedback loops can be identified that span two or more physiological signals. In Figure 7 there is clearly a strong coupling across heart rate, electrodermal activity and skin temperature that is completely driven by the client. A second motif connects heart rate and temperature to acceleration, which are both driven by the caregiver. What cannot be determined based on the current analyses is a temporal precedence for the emergence of these feedback loops. That is, we cannot claim the caregiver’s movements as measured by the acceleration variable drive the heart rate coupling dynamics which causes a reinforcing feedback loop in the client involving heartrate, electrodermal activity and temperature, that raise the stress levels of both client and caregiver to a level beyond which an incident is inevitable. The opposite scenario could also be true, in which the heart rate of the client causes the caregiver to drive the acceleration interaction dynamics, which in turn reinforces the driving of the heartrate dynamics by the client via temperature.

Future studies should explore whether imposing temporal constraints on the recurrence networks, i.e., turning them into directed networks (Hasselman, 2022) can provide more detailed information about temporal precedence in coupling dynamics. In regular directed networks, cross-clustering coefficients can be calculated to represent different types of connected triples, such as cycles, fanning in, fanning out and middleman patterns (Fagiolo, 2007). It is also possible to create weighted (and/or directed) recurrence networks (Hasselman and Bosman, 2020), for example, by keeping the distance values of the points that fall below the recurrence threshold $\epsilon$ in the matrix. Another option is to use recurrence times to give more weight to triangles representing states that will recur sooner rather than later in the calculation of the cross-clustering coefficient.

To summarize, the inter-system recurrence network approach to the analysis of coupled, multivariate systems produces estimates of coupling direction that are interpretable relative to the measurement context and provide insights into complex dependencies between different variables. These initial results suggest that IRN measures can at least inform post-hoc evaluations of interaction dynamics, which, in the present context, could occur as soon as the Empatica4 data become available. More research is required to determine whether the method can eventually be used to inform decisions related to intervention/prevention of behavior during real-time process monitoring of physiological and/or psychological variables.

Data availability statement

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

Ethics statement

The studies involving humans were approved by Faculty Ethics Committee of the Radboud University (NL-number: NL71683.091.19) 12 February 2020. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

FH: Conceptualization, Methodology, Visualization, Writing–original draft, Writing–review and editing. LU: Data curation, Project administration, Writing–review and editing. RK: Data curation, Project administration, Writing–review and editing. PL: Data curation, Methodology, Project administration, Writing–review and editing. RO: Data curation, Project administration, Writing–review and editing.

Funding

The authors declare that no financial support was received for the research, authorship, and/or publication of this article.

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.

References

Barrat, A., Barthélemy, M., Pastor-Satorras, R., and Vespignani, A. (2004). The architecture of complex weighted networks. Proceedings of the National Academy of Sciences of the United States of America.

Google Scholar

Cox, R. F., van der Steen, S., Guevara, M., de Jonge-Hoekstra, L., and van Dijk, M. (2016). “Chromatic and anisotropic cross-recurrence quantification analysis of interpersonal behavior,” in Recurrence plots and their quantifications: expanding horizons (Springer), 209–225.

CrossRef Full Text | Google Scholar

de Looff, P., Duursma, R., Noordzij, M., Taylor, S., Jaques, N., Scheepers, F., et al. (2022). Wearables: an R package with accompanying shiny application for signal analysis of a wearable device targeted at clinicians and researchers. Front. Behav. Neurosci. 16, 856544. doi:10.3389/fnbeh.2022.856544

PubMed Abstract | CrossRef Full Text | Google Scholar

Devereux, J., Hastings, R., and Noone, S. (2009). Staff stress and burnout in intellectual disability services: work stress theory and its application. J. Appl. Res. Intellect. Disabil. 22 (6), 561–573. doi:10.1111/j.1468-3148.2009.00509.x

CrossRef Full Text | Google Scholar

Donges, J. F., Schultz, H. C. H., Marwan, N., Zou, Y., and Kurths, J. (2011). Investigating the topology of interacting networks—theory and application to coupled climate subnetworks. Eur. Phys. J. B 84 (4), 635–651. Article 4. doi:10.1140/epjb/e2011-10795-8

CrossRef Full Text | Google Scholar

Fagiolo, G. (2007). Clustering in complex directed networks. Phys. Rev. E 76 (2), 026107. doi:10.1103/PhysRevE.76.026107

PubMed Abstract | CrossRef Full Text | Google Scholar

Fartacek, C., Schiepek, G., Kunrath, S., Fartacek, R., and Plöderl, M. (2016). Real-time monitoring of non-linear suicidal dynamics: methodology and a demonstrative case report. Front. Psychol. 7, 130. doi:10.3389/fpsyg.2016.00130

PubMed Abstract | CrossRef Full Text | Google Scholar

Feldhoff, J. H., Donner, R. V., Donges, J. F., Marwan, N., and Kurths, J. (2012). Geometric detection of coupling directions by means of inter-system recurrence networks. Phys. Lett. A 376 (46), 3504–3513. doi:10.1016/j.physleta.2012.10.008

CrossRef Full Text | Google Scholar

Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37 (3), 424–438. doi:10.2307/1912791

CrossRef Full Text | Google Scholar

Hasselman, F. (2022). Early warning signals in phase space: geometric resilience loss indicators from multiplex cumulative recurrence networks. Front. Physiology 13, 859127. doi:10.3389/fphys.2022.859127

CrossRef Full Text | Google Scholar

Hasselman, F., and Bosman, A. M. T. (2020). Studying complex adaptive systems with internal states: a recurrence network approach to the analysis of multivariate time-series data representing self-reports of human experience. Front. Appl. Math. Statistics 6. doi:10.3389/fams.2020.00009

CrossRef Full Text | Google Scholar

Janssen, C. G. C., Schuengel, C., and Stolk, J. (2002). Understanding challenging behaviour in people with severe and profound intellectual disability: a stress-attachment model. J. Intellect. Disabil. Res. 46 (6), 445–453. doi:10.1046/j.1365-2788.2002.00430.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Konvalinka, I., Xygalatas, D., Bulbulia, J., Schjødt, U., Jegindø, E.-M., Wallot, S., et al. (2011). Synchronized arousal between performers and related spectators in a fire-walking ritual. Proc. Natl. Acad. Sci. U. S. A. 108, 8514–8519. doi:10.1073/pnas.1016955108

PubMed Abstract | CrossRef Full Text | Google Scholar

Leng, S., Ma, H., Kurths, J., Lai, Y.-C., Lin, W., Aihara, K., et al. (2020). Partial cross mapping eliminates indirect causal influences. Nat. Commun. 11 (1), 2632. doi:10.1038/s41467-020-16238-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, H., Aihara, K., and Chen, L. (2014). Detecting causality from nonlinear dynamics with short-term time series. Sci. Rep. 4 (1), 7464. Article 1. doi:10.1038/srep07464

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, H., Leng, S., Tao, C., Ying, X., Kurths, J., Lai, Y.-C., et al. (2017). Detection of time delays and directional interactions based on time series from complex dynamical systems. Phys. Rev. E 96 (1), 012221. doi:10.1103/PhysRevE.96.012221

PubMed Abstract | CrossRef Full Text | Google Scholar

Main, A., Paxton, A., and Dale, R. (2016). An exploratory analysis of emotion dynamics between mothers and adolescents during conflict discussions. Emotion 16 (6), 913–928. doi:10.1037/emo0000180

PubMed Abstract | CrossRef Full Text | Google Scholar

Marinazzo, D., Liao, W., Chen, H., and Stramaglia, S. (2011). Nonlinear connectivity by Granger causality. NeuroImage 58 (2), 330–338. doi:10.1016/j.neuroimage.2010.01.099

PubMed Abstract | CrossRef Full Text | Google Scholar

Marwan, N., and Kurths, J. (2002). Nonlinear analysis of bivariate data with cross recurrence plots. Phys. Lett. A 302 (5–6), 299–307. doi:10.1016/s0375-9601(02)01170-2

CrossRef Full Text | Google Scholar

Marwan, N., Romano, M. C., Thiel, M., and Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Phys. Reports-Review Sect. Phys. Lett. 438 (5–6), 237–329. doi:10.1016/j.physrep.2006.11.001

CrossRef Full Text | Google Scholar

Milstein, N., and Gordon, I. (2020). Validating measures of electrodermal activity and heart rate variability derived from the Empatica E4 utilized in research settings that involve interactive dyadic states. Front. Behav. Neurosci. 14, 148. doi:10.3389/fnbeh.2020.00148

PubMed Abstract | CrossRef Full Text | Google Scholar

Olthof, M., Hasselman, F., Aas, B., Lamoth, D., Scholz, S., Daniels-Wredenhagen, N., et al. (2023). The best of both worlds? General principles of psychopathology in personalized assessmen. J. Psychopathol. Clin. Sci. 132, 808–819. doi:10.1037/abn0000858

PubMed Abstract | CrossRef Full Text | Google Scholar

Olthof, M., Hasselman, F., Strunk, G., van Rooij, M., Aas, B., Helmich, M. A., et al. (2019). Critical fluctuations as an early-warning signal for sudden gains and losses in patients receiving psychotherapy for mood disorders. Clin. Psychol. Sci. 8 (1), 25–35. doi:10.1177/2167702619865969

CrossRef Full Text | Google Scholar

Paluš, M., Komárek, V., Hrnčíř, Z., and Štěrbová, K. (2001). Synchronization as adjustment of information rates: detection from bivariate time series. Phys. Rev. E 63 (4), 046211. doi:10.1103/PhysRevE.63.046211

CrossRef Full Text | Google Scholar

Petzold, L. (1983). Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM J. Sci. Stat. Comput. 4 (1), 136–148. doi:10.1137/0904010

CrossRef Full Text | Google Scholar

Quax, R., Har-Shemesh, O., and Sloot, P. M. A. (2017). Quantifying synergistic information using intermediate stochastic variables. Entropy 19 (2), 85. Article 2. doi:10.3390/e19020085

CrossRef Full Text | Google Scholar

R Core Team (2022). R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.

Google Scholar

Richardson, D. C., Dale, R., and Kirkham, N. Z. (2007). The art of conversation is coordination: common ground and the coupling of eye movements during dialogue. Psychol. Sci. 18 (5), 407–413. doi:10.1111/j.1467-9280.2007.01914.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Runge, J., Bathiany, S., Bollt, E., Camps-Valls, G., Coumou, D., Deyle, E., et al. (2019). Inferring causation from time series in Earth system sciences. Nat. Commun. 10 (1), 2553. Article 1. doi:10.1038/s41467-019-10105-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiepek, G., Stöger-Schmidinger, B., Aichhorn, W., Schöller, H., and Aas, B. (2016). Systemic case formulation, individualized process monitoring, and state dynamics in a case of dissociative identity disorder. Front. Psychol. 7, 1545. doi:10.3389/fpsyg.2016.01545

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Schuurmans, A. A. T., de Looff, P., Nijhof, K. S., Rosada, C., Scholte, R. H. J., Popma, A., et al. (2020). Validity of the Empatica E4 wristband to measure heart rate variability (hrv) parameters: a comparison to electrocardiography (ecg). J. Med. Syst. 44 (11), 190. doi:10.1007/s10916-020-01648-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Shockley, K. (2005). “Cross recurrence quantification of interpersonal postural activity,” in Tutorials in contemporary nonlinear methods for the behavioral sciences, 142–177. Available at: http://nsf.gov/sbe/bcs/pac/nmbs/chap4.pdf.

Google Scholar

Shockley, K., Butwill, M., Zbilut, J. P., and Webber, C. L. (2002). Cross recurrence quantification of coupled oscillators. Phys. Lett. A 305 (1–2), 59–69. doi:10.1016/s0375-9601(02)01411-1

CrossRef Full Text | Google Scholar

Simons, R., Koordeman, R., de Looff, P., and Otten, R. (2021). Physiological measurements of stress preceding incidents of challenging behavior in people with severe to profound intellectual disabilities: longitudinal study protocol of single-case studies. JMIR Res. Protoc. 10, e24911. doi:10.2196/24911

PubMed Abstract | CrossRef Full Text | Google Scholar

Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010). Solving differential equations in R: package deSolve. J. Stat. Softw. 33, 1–25. doi:10.18637/jss.v033.i09

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugihara, G., May, R., Ye, H., Hsieh, C., Deyle, E., Fogarty, M., et al. (2012). Detecting causality in complex ecosystems. Sci. (New York, N.Y.) 338 (6106), 496–500. doi:10.1126/science.1227079

PubMed Abstract | CrossRef Full Text | Google Scholar

Vink, R., Hasselman, F., Cillessen, A. H., Wijnants, M. L., and Bosman, A. M. (2018). Does competence determine who leads in a dyadic cooperative task? A study of children with and without a neurodevelopmental disorder. Complexity.

Google Scholar

Wallot, S., and Leonardi, G. (2018). Analyzing multivariate dynamics using cross-recurrence quantification analysis (CRQA), diagonal-cross-recurrence Profiles (DCRP), and multidimensional recurrence quantification analysis (MdRQA) – a tutorial in R. Front. Psychol. 9, 2232. doi:10.3389/fpsyg.2018.02232

PubMed Abstract | CrossRef Full Text | Google Scholar

Zbilut, J. P., Giuliani, A., and Webber, C. L. (1998). Detecting deterministic signals in exceptionally noisy environments using cross-recurrence quantification. Phys. Lett. A 246 (1–2), 122–128. doi:10.1016/s0375-9601(98)00457-5

CrossRef Full Text | Google Scholar

Zou, Y., Donner, R. V., Marwan, N., Donges, J. F., Kurths, J., and Yao, Q. (2019). MiR-146a attenuates liver fibrosis by inhibiting transforming growth factor-β1 mediated epithelial-mesenchymal transition in hepatocytes. Phys. Rep. 787, 1–8. doi:10.1016/j.cellsig.2019.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar