| Title: | Network Physiology and Time-Delay-Stability Organ-Interaction Networks |
|---|---|
| Description: | Tools for Network Physiology: the study of dynamic interactions among diverse physiological systems (brain, heart, respiration, muscle). Implements Time Delay Stability (TDS) analysis to identify stable dynamic couplings between signals, builds organ-interaction networks whose links are significance-tested against phase-randomised surrogates, and quantifies how those networks reconfigure across physiological states (e.g. sleep stages). Designed to operate on the multimodal PhysioExperiment data model, but the core engine works on plain multivariate time series. Methods follow Bashan et al. (2012) <doi:10.1038/ncomms1705> and Bartsch et al. (2015) <doi:10.1371/journal.pone.0142143>. |
| Authors: | Yusuke Matsui [aut, cre] |
| Maintainer: | Yusuke Matsui <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.1 |
| Built: | 2026-10-03 13:42:22 UTC |
| Source: | https://github.com/x-biosignal/PhysioNetPhysiology |
Aggregate a multiplex into a single weighted network
aggregateNetwork(x, method = c("sum", "mean"))aggregateNetwork(x, method = c("sum", "mean"))
x |
A |
method |
|
A single weighted adjacency matrix.
L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) aggregateNetwork(multiplexNetwork(list(alpha = L1, beta = L2)))L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) aggregateNetwork(multiplexNetwork(list(alpha = L1, beta = L2)))
Builds a continuous cardiac phase that advances linearly from 0 to
between consecutive R-peaks (wrapped per beat). Samples outside the
first/last beat are NA.
cardiacPhaseFromPeaks(rpeaks, n)cardiacPhaseFromPeaks(rpeaks, n)
rpeaks |
Integer sample indices of R-peaks (1-based). |
n |
Length of the output phase vector (number of samples). |
Numeric vector of wrapped cardiac phase in [0, 2*pi), length n.
rpeaks <- seq(10, 1000, by = 25) ph <- cardiacPhaseFromPeaks(rpeaks, 1000)rpeaks <- seq(10, 1000, by = 25) ph <- cardiacPhaseFromPeaks(rpeaks, 1000)
High-level convenience wrapper: extracts respiratory phase (band-passed Hilbert) and cardiac phase (from R-peaks), scans n:m phase locking, and returns the dominant heartbeats-per-breath ratio, its locking index, and the synchrogram.
cardiorespiratoryCoupling( resp, rpeaks, sampling_rate, resp_band = c(0.1, 0.5), max_ratio = 8L )cardiorespiratoryCoupling( resp, rpeaks, sampling_rate, resp_band = c(0.1, 0.5), max_ratio = 8L )
resp |
Respiration signal. |
rpeaks |
Integer R-peak sample indices (same time base as |
sampling_rate |
Sampling rate in Hz. |
resp_band |
Respiratory passband in Hz (default |
max_ratio |
Largest heartbeats-per-breath ratio to consider (default 8). |
An object of class "cardioresp": a list with ratio (a
"heart:resp" string), n, m, lambda (locking index), sync_index
(1:1 phase-locking value), synchrogram, and the table of scanned ratios.
t <- seq(0, 120, by = 0.01); fs <- 100 resp <- sin(2 * pi * 0.25 * t) card <- sin(2 * pi * 1.0 * t) # 4 beats per breath rpeaks <- which(diff(sign(card)) > 0) cr <- cardiorespiratoryCoupling(resp, rpeaks, fs) cr$ratiot <- seq(0, 120, by = 0.01); fs <- 100 resp <- sin(2 * pi * 0.25 * t) card <- sin(2 * pi * 1.0 * t) # 4 beats per breath rpeaks <- which(diff(sign(card)) > 0) cr <- cardiorespiratoryCoupling(resp, rpeaks, fs) cr$ratio
Partitions a network into communities (modules) and reports the modularity Q.
The membership can be passed to moduleRoles() / participationCoefficient(),
closing the loop those role metrics need. Delegates to igraph.
communityDetection( x, method = c("louvain", "fast_greedy", "walktrap", "leading_eigen", "infomap"), seed = NULL )communityDetection( x, method = c("louvain", "fast_greedy", "walktrap", "leading_eigen", "infomap"), seed = NULL )
x |
A |
method |
Community-detection algorithm: |
seed |
Optional integer seed (some algorithms are stochastic). |
A community_partition object: membership (named integer vector),
modularity (Q), n_communities, sizes, and method.
networkMetrics(), moduleRoles()
set.seed(1) A <- matrix(0, 8, 8) A[1:4, 1:4] <- 0.9; A[5:8, 5:8] <- 0.9; A[1, 5] <- A[5, 1] <- 0.2 diag(A) <- 0 communityDetection(A)set.seed(1) A <- matrix(0, 8, 8) A[1:4, 1:4] <- 0.9; A[5:8, 5:8] <- 0.9; A[1, 5] <- A[5, 1] <- 0.2 diag(A) <- 0 communityDetection(A)
Epochs a brain signal around each R-peak, baseline-corrects, and averages to the heartbeat-evoked potential.
heartbeatEvokedPotential( eeg, rpeaks, sampling_rate, window = c(-0.1, 0.6), baseline = c(-0.1, 0) )heartbeatEvokedPotential( eeg, rpeaks, sampling_rate, window = c(-0.1, 0.6), baseline = c(-0.1, 0) )
eeg |
Numeric brain signal (one channel). |
rpeaks |
Integer R-peak sample indices (same time base as |
sampling_rate |
Sampling rate in Hz. |
window |
Epoch window |
baseline |
Baseline window |
An object of class "hep": a list with time (s, relative to the
R-peak), hep (averaged waveform), epochs (epochs x time matrix), n
(number of epochs), and sampling_rate.
set.seed(1); fs <- 250; n <- 30 * fs rpeaks <- seq(fs, n - fs, by = round(0.9 * fs)) eeg <- rnorm(n) hep <- heartbeatEvokedPotential(eeg, rpeaks, fs)set.seed(1); fs <- 250; n <- 30 * fs rpeaks <- seq(fs, n - fs, by = round(0.9 * fs)) eeg <- rnorm(n) hep <- heartbeatEvokedPotential(eeg, rpeaks, fs)
Amplitude of the HEP in a component window
hepAmplitude(hep, component = c(0.2, 0.4), fun = c("mean", "peak"))hepAmplitude(hep, component = c(0.2, 0.4), fun = c("mean", "peak"))
hep |
A |
component |
|
fun |
Summary over the window: |
A single amplitude value.
set.seed(1); fs <- 250; n <- 30 * fs rpeaks <- seq(fs, n - fs, by = round(0.9 * fs)) hep <- heartbeatEvokedPotential(rnorm(n), rpeaks, fs) hepAmplitude(hep, component = c(0.2, 0.4))set.seed(1); fs <- 250; n <- 30 * fs rpeaks <- seq(fs, n - fs, by = round(0.9 * fs)) hep <- heartbeatEvokedPotential(rnorm(n), rpeaks, fs) hepAmplitude(hep, component = c(0.2, 0.4))
Tests whether the HEP component amplitude is larger than expected if the brain
signal were not time-locked to the heartbeat, by recomputing the component
amplitude for n_perm datasets with circularly shifted R-peaks.
hepSignificance( eeg, rpeaks, sampling_rate, window = c(-0.1, 0.6), baseline = c(-0.1, 0), component = c(0.2, 0.4), n_perm = 200 )hepSignificance( eeg, rpeaks, sampling_rate, window = c(-0.1, 0.6), baseline = c(-0.1, 0), component = c(0.2, 0.4), n_perm = 200 )
eeg, rpeaks, sampling_rate, window, baseline
|
As in
|
component |
Component window in seconds (see |
n_perm |
Number of permutations (default 200). |
A list with observed (|component amplitude|), null (permutation
amplitudes), and p_value (one-sided).
set.seed(1); fs <- 250; n <- 20 * fs rpeaks <- seq(fs, n - fs, by = round(0.9 * fs)) hepSignificance(rnorm(n), rpeaks, fs, n_perm = 50)$p_valueset.seed(1); fs <- 250; n <- 20 * fs rpeaks <- seq(fs, n - fs, by = round(0.9 * fs)) hepSignificance(rnorm(n), rpeaks, fs, n_perm = 50)$p_value
Analytic-signal (Hilbert) instantaneous phase, optionally after isolating a
frequency band. Physiological oscillations (respiration, a cardiac rhythm)
should be reasonably narrow-band for the phase to be meaningful; pass band
to bandpass first.
instantaneousPhase(x, sampling_rate = 1, band = NULL)instantaneousPhase(x, sampling_rate = 1, band = NULL)
x |
Numeric signal. |
sampling_rate |
Sampling rate in Hz (default 1). |
band |
Optional |
Numeric vector of wrapped phases in (-pi, pi].
t <- seq(0, 10, by = 0.01) ph <- instantaneousPhase(sin(2 * pi * 1 * t), sampling_rate = 100)t <- seq(0, 10, by = 0.01) ph <- instantaneousPhase(sin(2 * pi * 1 * t), sampling_rate = 100)
Cosine similarity between the (upper-triangular) edge-weight vectors of each pair of layers: 1 = identical coupling structure, 0 = disjoint.
layerSimilarity(x)layerSimilarity(x)
x |
A |
A n_layers x n_layers similarity matrix.
L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) layerSimilarity(multiplexNetwork(list(alpha = L1, beta = L2)))L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) layerSimilarity(multiplexNetwork(list(alpha = L1, beta = L2)))
Combines the within-module degree z-score and participation coefficient into the seven-role classification of Guimera & Amaral (2005).
moduleRoles(adj, communities)moduleRoles(adj, communities)
adj |
A square adjacency matrix. |
communities |
A vector of module labels, one per node. |
A data frame with node, z (within-module degree z-score), P
(participation coefficient), and role.
adj <- matrix(0, 6, 6) adj[1, 2:3] <- adj[2:3, 1] <- 1 adj[4, 5] <- adj[5, 4] <- adj[5, 6] <- adj[6, 5] <- 1 adj[3, 4] <- adj[4, 3] <- 1 moduleRoles(adj, c(1, 1, 1, 2, 2, 2))adj <- matrix(0, 6, 6) adj[1, 2:3] <- adj[2:3, 1] <- 1 adj[4, 5] <- adj[5, 4] <- adj[5, 6] <- adj[6, 5] <- 1 adj[3, 4] <- adj[4, 3] <- 1 moduleRoles(adj, c(1, 1, 1, 2, 2, 2))
Per-layer and overlapping (multiplex) degree
multiplexDegree(x)multiplexDegree(x)
x |
A |
A list with per_layer (nodes x layers degree matrix) and
overlapping (summed degree across layers per node).
L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) mp <- multiplexNetwork(list(alpha = L1, beta = L2)) multiplexDegree(mp)L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) mp <- multiplexNetwork(list(alpha = L1, beta = L2)) multiplexDegree(mp)
Construct a multiplex network
multiplexNetwork(layers, node_names = NULL)multiplexNetwork(layers, node_names = NULL)
layers |
A named list of adjacency matrices, all the same size and node order (one layer per coupling channel or frequency band). |
node_names |
Optional node names (default: the row names of the first
layer, or |
An object of class "multiplex": a list with layers, nodes,
n_layers, layer_names.
L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) mp <- multiplexNetwork(list(alpha = L1, beta = L2))L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) mp <- multiplexNetwork(list(alpha = L1, beta = L2))
How evenly a node's connections spread across layers. 0 = all links confined to one layer; 1 = links distributed uniformly across all layers (Battiston et al. 2014).
multiplexParticipation(x)multiplexParticipation(x)
x |
A |
A named numeric vector per node, in [0, 1].
L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) multiplexParticipation(multiplexNetwork(list(alpha = L1, beta = L2)))L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) multiplexParticipation(multiplexNetwork(list(alpha = L1, beta = L2)))
Computes the standard network-science metrics on a weighted physiological
network: node-level (degree, strength, weighted clustering, betweenness,
closeness, eigenvector centrality, local efficiency) and global (density,
transitivity, characteristic path length, global/local efficiency, degree
assortativity, and modularity from a community partition). Edge weights are
treated as coupling strengths; path-based metrics use distances 1/weight
(stronger coupling = shorter path). Delegates to igraph.
networkMetrics(x, directed = FALSE, community_method = "louvain")networkMetrics(x, directed = FALSE, community_method = "louvain")
x |
A |
directed |
Treat the network as directed (default |
community_method |
Community algorithm for the modularity summary
(see |
A network_metrics object: global (named list), node (data frame
of per-node metrics), membership, and directed.
communityDetection(), smallWorldness(), moduleRoles()
set.seed(1) A <- matrix(runif(64, 0, 0.3), 8, 8); A <- (A + t(A)) / 2; diag(A) <- 0 nm <- networkMetrics(A) nm$global$global_efficiencyset.seed(1) A <- matrix(runif(64, 0, 0.3), 8, 8); A <- (A + t(A)) / 2; diag(A) <- 0 nm <- networkMetrics(A) nm$global$global_efficiency
Computes for every coprime (n, m) in the requested range
and returns the dominant ratio. The frequency ratio f1 : f2 implied by the
winning (n, m) is m : n.
nmPhaseSynchronization(phi1, phi2, orders_n = 1:8, orders_m = 1:8)nmPhaseSynchronization(phi1, phi2, orders_n = 1:8, orders_m = 1:8)
phi1, phi2
|
Instantaneous phases. |
orders_n, orders_m
|
Integer vectors of orders to scan. |
A list with table (data frame of n, m, lambda), and the
dominant n, m, lambda, and ratio (a "f1:f2" string = m:n).
t <- seq(0, 120, by = 0.01) resp <- sin(2 * pi * 0.25 * t) card <- sin(2 * pi * 1.0 * t) # 4x faster: locked 4:1 pr <- instantaneousPhase(resp, 100) pc <- instantaneousPhase(card, 100) nmPhaseSynchronization(pc, pr)$ratiot <- seq(0, 120, by = 0.01) resp <- sin(2 * pi * 0.25 * t) card <- sin(2 * pi * 1.0 * t) # 4x faster: locked 4:1 pr <- instantaneousPhase(resp, 100) pc <- instantaneousPhase(card, 100) nmPhaseSynchronization(pc, pr)$ratio
, where is node
i's summed connection weight to module s and its total strength.
0 = all links inside one module; approaches 1 as links spread evenly over
many modules.
participationCoefficient(adj, communities)participationCoefficient(adj, communities)
adj |
A square (weighted or binary) adjacency matrix. |
communities |
A vector of module labels, one per node. |
A named numeric vector of participation coefficients in [0, 1).
adj <- matrix(0, 4, 4); adj[1,2] <- adj[2,1] <- 1; adj[3,4] <- adj[4,3] <- 1 participationCoefficient(adj, c(1, 1, 2, 2))adj <- matrix(0, 4, 4); adj[1,2] <- adj[2,1] <- 1; adj[3,4] <- adj[4,3] <- 1 participationCoefficient(adj, c(1, 1, 2, 2))
The n:m phase-synchronization index
, in [0, 1].
For two rhythms whose frequency ratio is f1 : f2 = m : n, the term
n*phi1 - m*phi2 is (near-)constant and approaches 1.
phaseSynchronizationIndex(phi1, phi2, n = 1L, m = 1L)phaseSynchronizationIndex(phi1, phi2, n = 1L, m = 1L)
phi1, phi2
|
Instantaneous phases (e.g. from |
n, m
|
Integer orders (default 1:1). |
A single locking index in [0, 1].
t <- seq(0, 60, by = 0.01) p1 <- instantaneousPhase(sin(2 * pi * 1 * t), 100) p2 <- instantaneousPhase(sin(2 * pi * 1 * t), 100) phaseSynchronizationIndex(p1, p2)t <- seq(0, 60, by = 0.01) p1 <- instantaneousPhase(sin(2 * pi * 1 * t), 100) p2 <- instantaneousPhase(sin(2 * pi * 1 * t), 100) phaseSynchronizationIndex(p1, p2)
Builds the time-by-nodes matrix that timeDelayStability() consumes, where
each node is one physiological variable (an organ system or a derived feature
such as an EEG band-power, heart rate, respiration amplitude, or an EMG
envelope). Inputs are a named list of equal-length numeric vectors, a numeric
matrix, or a data frame.
physioNodeMatrix( x, nodes = NULL, standardize = FALSE, assay = NULL, channels = NULL, feature = c("raw", "envelope"), target_rate = NULL )physioNodeMatrix( x, nodes = NULL, standardize = FALSE, assay = NULL, channels = NULL, feature = c("raw", "envelope"), target_rate = NULL )
x |
A named list of equal-length numeric vectors, a numeric matrix
(time x nodes), a data frame of numeric columns, or a |
nodes |
Optional character vector of node names (defaults to the names of
|
standardize |
If |
assay |
For a |
channels |
For a |
feature |
Per-channel feature: |
target_rate |
Optional rate (Hz) to down-sample each channel to, by
block-averaging. Required for a |
In Network Physiology, coupling is assessed between the dynamics of these
variables, so nodes are typically instantaneous features sampled on a common,
relatively slow grid (e.g. 1 Hz), not the raw high-rate waveforms. Assemble
such features with the relevant Physio* package and pass them here.
A numeric matrix, time (rows) by nodes (columns), with column names.
nodes <- list( hr = cumsum(rnorm(500)), resp = sin(seq(0, 20 * pi, length.out = 500)), emg = abs(rnorm(500)) ) M <- physioNodeMatrix(nodes) dim(M)nodes <- list( hr = cumsum(rnorm(500)), resp = sin(seq(0, 20 * pi, length.out = 500)), emg = abs(rnorm(500)) ) M <- physioNodeMatrix(nodes) dim(M)
Draws the classic Time-Delay-Stability trace: the per-window optimal lag with windows inside stable periods highlighted. Long flat stretches are the stable couplings TDS rewards.
plotDelayStability(x, i, j, ...)plotDelayStability(x, i, j, ...)
x |
A |
i, j
|
The two nodes (names or indices) whose coupling to plot. |
... |
Passed to |
invisible(NULL); called for the plot.
Plot a heartbeat-evoked potential
plotHEP(x, main = "Heartbeat-evoked potential", ...)plotHEP(x, main = "Heartbeat-evoked potential", ...)
x |
A |
main |
Title. |
... |
Passed to |
invisible(NULL).
Plot a cardiorespiratory synchrogram
plotSynchrogram( x, sampling_rate = NULL, main = "Cardiorespiratory synchrogram", ... )plotSynchrogram( x, sampling_rate = NULL, main = "Cardiorespiratory synchrogram", ... )
x |
A |
sampling_rate |
Sampling rate, to label the x-axis in seconds (optional). |
main |
Plot title. |
... |
Passed to |
invisible(NULL).
Circular node layout; edge width and colour encode TDS strength.
plotTDSnetwork(x, main = "TDS organ-interaction network", ...)plotTDSnetwork(x, main = "TDS organ-interaction network", ...)
x |
A |
main |
Plot title. |
... |
Ignored. |
invisible(NULL).
Heatmap of per-state node degree (states x nodes): how each system's number of stable couplings changes across physiological states.
plotTDSreconfiguration(x, main = "TDS network reconfiguration", ...)plotTDSreconfiguration(x, main = "TDS network reconfiguration", ...)
x |
A |
main |
Plot title. |
... |
Ignored. |
invisible(NULL).
Quantifies small-world organization on the binarized network: the
Humphries-Gurney sigma = (C/C_rand)/(L/L_rand) (small-world when > 1, comparing
clustering and path length to degree-preserving random graphs) and the Telford
omega = L_rand/L - C/C_latt (near 0 = small-world; -> 1 random; -> -1 lattice).
Delegates the randomization and lattice references to igraph.
smallWorldness(x, n_rand = 50L, seed = NULL)smallWorldness(x, n_rand = 50L, seed = NULL)
x |
A |
n_rand |
Number of degree-preserving random graphs (default 50). |
seed |
Optional integer seed for the randomization. |
A small_world object: sigma, omega, and the components C, L,
C_rand, L_rand, C_latt, n_rand.
Humphries & Gurney 2008; Telford et al. 2011.
g <- igraph::sample_smallworld(1, 40, 4, 0.05) # a small-world graph sw <- smallWorldness(as.matrix(igraph::as_adjacency_matrix(g)), n_rand = 20, seed = 1) sw$sigmag <- igraph::sample_smallworld(1, 40, 4, 0.05) # a small-world graph sw <- smallWorldness(as.matrix(igraph::as_adjacency_matrix(g)), n_rand = 20, seed = 1) sw$sigma
Builds the (N*L) x (N*L) block matrix with the layer adjacencies on the
diagonal blocks and interlayer coupling on the identity off-diagonal blocks
(categorical coupling of a node to its own replicas across layers).
supraAdjacency(x, interlayer = 1)supraAdjacency(x, interlayer = 1)
x |
A |
interlayer |
Interlayer coupling weight (default 1). |
A (N*L) x (N*L) numeric matrix.
L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) supraAdjacency(multiplexNetwork(list(alpha = L1, beta = L2)))L1 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3) L2 <- matrix(c(0,0,1, 0,0,0, 1,0,0), 3) supraAdjacency(multiplexNetwork(list(alpha = L1, beta = L2)))
The classic Schafer-Rosenblum synchrogram: the respiratory phase observed at
each heartbeat, wrapped into m respiratory cycles. Under n:m locking (n
heartbeats per m breaths) the points collapse onto n horizontal bands.
synchrogram(resp_phase, rpeaks, m = 1L)synchrogram(resp_phase, rpeaks, m = 1L)
resp_phase |
Instantaneous respiratory phase (length = n samples), e.g.
from |
rpeaks |
Integer R-peak sample indices. |
m |
Number of respiratory cycles to wrap into (default 1). |
A data frame with beat (R-peak sample index) and psi (respiratory
phase at that beat, in [0, m)), of class "synchrogram".
t <- seq(0, 120, by = 0.01); fs <- 100 resp_phase <- instantaneousPhase(sin(2 * pi * 0.25 * t), fs) rpeaks <- which(diff(sign(sin(2 * pi * 1.0 * t))) > 0) # ~1 Hz beats sg <- synchrogram(resp_phase, rpeaks, m = 1)t <- seq(0, 120, by = 0.01); fs <- 100 resp_phase <- instantaneousPhase(sin(2 * pi * 0.25 * t), fs) rpeaks <- which(diff(sign(sin(2 * pi * 1.0 * t))) > 0) # ~1 Hz beats sg <- synchrogram(resp_phase, rpeaks, m = 1)
Thresholds a TDS matrix into a weighted, undirected network: nodes are
physiological systems and edges are stable dynamic couplings whose weight is
the TDS percentage. Edges below threshold are removed.
tdsNetwork( x, threshold = NULL, surrogate = FALSE, X = NULL, n_surrogates = 50, quantile = 0.95 )tdsNetwork( x, threshold = NULL, surrogate = FALSE, X = NULL, n_surrogates = 50, quantile = 0.95 )
x |
A |
threshold |
Minimum TDS percentage for an edge to be kept. If |
surrogate |
Logical; if |
X |
Original node matrix, required when |
n_surrogates, quantile
|
Passed to |
An object of class "tds_network": a list with adjacency (weighted,
thresholded), threshold, nodes, degree (number of links per node),
strength (summed edge weight per node), and edges (a data frame of the
surviving links).
set.seed(1) n <- 3000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0 y <- c(rep(0, 5), x[seq_len(n - 5)]) + rnorm(n, sd = 0.3) z <- rnorm(n) X <- cbind(brain = x, heart = y, muscle = z) res <- timeDelayStability(X, window = 100, step = 50, max_lag = 20, min_stable = 3) net <- tdsNetwork(res, threshold = 40) net$edgesset.seed(1) n <- 3000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0 y <- c(rep(0, 5), x[seq_len(n - 5)]) + rnorm(n, sd = 0.3) z <- rnorm(n) X <- cbind(brain = x, heart = y, muscle = z) res <- timeDelayStability(X, window = 100, step = 50, max_lag = 20, min_stable = 3) net <- tdsNetwork(res, threshold = 40) net$edges
Builds one TDS network per physiological state and summarises how the network reconfigures between states. TDS is computed continuously over the whole recording; each window is assigned to a state, and a pair's per-state TDS is the percentage of that state's windows that fall inside a stable coupling period.
tdsNetworkByState( X, states, sampling_rate = 1, window = NULL, step = NULL, max_lag = NULL, window_sec = NULL, step_sec = NULL, max_lag_sec = NULL, stability_tol = 1L, min_stable = 4L, threshold = NULL, n_surrogates = 50, min_windows = 3L )tdsNetworkByState( X, states, sampling_rate = 1, window = NULL, step = NULL, max_lag = NULL, window_sec = NULL, step_sec = NULL, max_lag_sec = NULL, stability_tol = 1L, min_stable = 4L, threshold = NULL, n_surrogates = 50, min_windows = 3L )
X |
Node matrix (or any input accepted by |
states |
A per-sample state label vector of length |
sampling_rate, window, step, max_lag, window_sec, step_sec, max_lag_sec, stability_tol, min_stable
|
Windowing / stability arguments passed to |
threshold |
Link threshold (TDS percentage). If |
n_surrogates |
Surrogates for the automatic threshold (default 50). |
min_windows |
Minimum number of windows a state must contain to be analysed (default 3); smaller states are dropped with a warning. |
An object of class "tds_state": a list with states (analysed
levels), tds (named list of per-state TDS matrices), networks (named
list of tdsNetwork() results), threshold, window_state, and
reconfiguration (a data frame of per-state link count, density, and mean
degree, plus a state x node degree matrix in degree).
set.seed(1) n <- 4000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0 half <- n / 2 # heart couples to brain only in state "B" y <- rnorm(n, sd = 0.3) y[(half + 6):n] <- x[half:(n - 6)] + rnorm(half - 5, sd = 0.3) X <- cbind(brain = x, heart = y) st <- rep(c("A", "B"), each = half) res <- tdsNetworkByState(X, st, window = 100, step = 50, max_lag = 20, min_stable = 3, threshold = 30) res$reconfigurationset.seed(1) n <- 4000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0 half <- n / 2 # heart couples to brain only in state "B" y <- rnorm(n, sd = 0.3) y[(half + 6):n] <- x[half:(n - 6)] + rnorm(half - 5, sd = 0.3) X <- cbind(brain = x, heart = y) st <- rep(c("A", "B"), each = half) res <- tdsNetworkByState(X, st, window = 100, step = 50, max_lag = 20, min_stable = 3, threshold = 30) res$reconfiguration
Returns, for every node pair, the TDS in each of two states and their difference — the edges that appear, vanish, or strengthen between states.
tdsReconfiguration(x, from, to)tdsReconfiguration(x, from, to)
x |
A |
from, to
|
State labels to compare. |
A data frame of from, to (node pair), the TDS in each state, and
delta, ordered by absolute change.
Generates n_surrogates phase-randomised surrogate datasets (each column
randomised independently, breaking cross-column timing while preserving each
column's power spectrum), recomputes the TDS matrix for each, and returns a
high quantile of the pooled off-diagonal surrogate TDS values. Links in the
real network exceeding this threshold are unlikely to arise from matched
spectra alone.
tdsSurrogateThreshold(X, n_surrogates = 50, quantile = 0.95, ...)tdsSurrogateThreshold(X, n_surrogates = 50, quantile = 0.95, ...)
X |
Node matrix (or any input accepted by |
n_surrogates |
Number of surrogate datasets (default 50). |
quantile |
Quantile of the surrogate TDS distribution to use as the threshold (default 0.95). |
... |
Windowing / stability arguments passed to |
A single numeric TDS percentage threshold, with attribute
"surrogate_tds" holding the pooled surrogate values.
set.seed(1) n <- 2000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0 y <- c(rep(0, 5), x[seq_len(n - 5)]) + rnorm(n, sd = 0.3) X <- cbind(x = x, y = y) thr <- tdsSurrogateThreshold(X, n_surrogates = 20, window = 100, step = 50, max_lag = 20, min_stable = 3) thrset.seed(1) n <- 2000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)); x[is.na(x)] <- 0 y <- c(rep(0, 5), x[seq_len(n - 5)]) + rnorm(n, sd = 0.3) X <- cbind(x = x, y = y) thr <- tdsSurrogateThreshold(X, n_surrogates = 20, window = 100, step = 50, max_lag = 20, min_stable = 3) thr
Computes pairwise Time Delay Stability (TDS) for every pair of columns in
X. Each column is treated as one physiological node (e.g. an EEG band-power,
heart rate, respiration, an EMG envelope). The result is a symmetric matrix
of TDS strengths (percentage of time each pair is stably coupled).
timeDelayStability( X, sampling_rate = 1, window = NULL, step = NULL, max_lag = NULL, window_sec = NULL, step_sec = NULL, max_lag_sec = NULL, stability_tol = 1L, min_stable = 4L, nodes = NULL )timeDelayStability( X, sampling_rate = 1, window = NULL, step = NULL, max_lag = NULL, window_sec = NULL, step_sec = NULL, max_lag_sec = NULL, stability_tol = 1L, min_stable = 4L, nodes = NULL )
X |
A numeric matrix or data frame, time (rows) by nodes (columns), or an
object accepted by |
sampling_rate |
Sampling rate in Hz of the rows of |
window, step, max_lag
|
Window length, step, and maximum lag in samples
(override the |
window_sec, step_sec, max_lag_sec
|
The same quantities in seconds.
Defaults: |
stability_tol |
Maximum change in optimal lag (in samples) between consecutive windows for the pair to be considered stable (default 1). |
min_stable |
Minimum number of consecutive windows forming a stable period (default 4). |
nodes |
Optional character vector of node names (defaults to the column
names of |
Windowing may be given in samples (window, step, max_lag) or in seconds
(window_sec, step_sec, max_lag_sec, using sampling_rate). Sample
arguments take precedence. Sensible defaults derive step and max_lag from
the window when unspecified.
An object of class "tds": a list with tds (nodes x nodes TDS
percentage matrix, diagonal 0), stable_lag (nodes x nodes dominant stable
lag in samples), pairs (per-pair detail incl. tau0 and in_stable),
nodes, n_windows, and params.
Bashan et al. (2012) doi:10.1038/ncomms1705.
set.seed(1) n <- 3000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)) x[is.na(x)] <- 0 y <- c(rep(0, 5), x[seq_len(n - 5)]) + rnorm(n, sd = 0.3) # y lags x by 5 z <- rnorm(n) # unrelated X <- cbind(x = x, y = y, z = z) res <- timeDelayStability(X, window = 100, step = 50, max_lag = 20, min_stable = 3) res$tdsset.seed(1) n <- 3000 x <- as.numeric(stats::filter(rnorm(n), rep(1/5, 5), sides = 2)) x[is.na(x)] <- 0 y <- c(rep(0, 5), x[seq_len(n - 5)]) + rnorm(n, sd = 0.3) # y lags x by 5 z <- rnorm(n) # unrelated X <- cbind(x = x, y = y, z = z) res <- timeDelayStability(X, window = 100, step = 50, max_lag = 20, min_stable = 3) res$tds
Standardises each node's within-module strength relative to the other nodes in its module. High values mark within-module hubs.
withinModuleDegreeZ(adj, communities)withinModuleDegreeZ(adj, communities)
adj |
A square adjacency matrix. |
communities |
A vector of module labels, one per node. |
A named numeric vector of z-scores.
adj <- matrix(0, 6, 6) adj[1, 2:3] <- adj[2:3, 1] <- 1 # module 1: node 1 hub adj[4, 5] <- adj[5, 4] <- adj[5, 6] <- adj[6, 5] <- 1 # module 2: 4-5-6 path adj[3, 4] <- adj[4, 3] <- 1 # bridge withinModuleDegreeZ(adj, c(1, 1, 1, 2, 2, 2))adj <- matrix(0, 6, 6) adj[1, 2:3] <- adj[2:3, 1] <- 1 # module 1: node 1 hub adj[4, 5] <- adj[5, 4] <- adj[5, 6] <- adj[6, 5] <- 1 # module 2: 4-5-6 path adj[3, 4] <- adj[4, 3] <- 1 # bridge withinModuleDegreeZ(adj, c(1, 1, 1, 2, 2, 2))