Multimodal Topological Characterization of Dynamic Causal Networks in Multivariate Time Series

1. Introduction

Internet-connected devices, social networks, and IoT sensors have generated massive volumes of multivariate time series. Their joint dynamics, capturing variables such as pressure and humidity together with temperature, provides a richer view of complex phenomena than univariate series alone [1,2].

Characterizing causal structures in multivariate time series has become a central challenge in neuroscience, economics, and environmental sciences. Classical methods such as Granger causality [3] provide estimators for the presence and direction of causal relations, yet the direct interpretation of the resulting p-value matrices is difficult in high-dimensional or non-stationary regimes [4,5], and these methods do not by themselves offer compact descriptors capable of detecting regime transitions or comparing causal structures over time [6,7].

The method proposed here does not aim to compete with consolidated causal discovery methods (Transfer Entropy [8], PCMCI [9], Convergent Cross Mapping [10], neural-network-based approaches [11,12]), but rather to complement them: while these estimators aim to identify which causal relations exist, our pipeline characterizes how these relations organize themselves topologically over time. We unify: (i) sliding-window estimation of Granger causality matrices, (ii) their transformation into complex networks via the Vector Visibility Graph (VVG) [13,14], and (iii) the analysis of these networks through three coordinated descriptors: betweenness centrality [15] captures local mediating structure; NetLSD [16] provides a global spectral view; and the Deterministic Tourist Walk (DTW) characterizes long-term dynamics through attractors and basins. We show that, although mathematically independent, these descriptors converge to a common critical threshold, providing an internal cross-validation that substitutes for comparison with external baselines in the specific task of topological characterization.

The remainder of the paper is organized as follows. The rest of Section 1 introduces the theoretical building blocks (Granger causality, VVG, betweenness centrality, NetLSD, and DTW). Section 2 reviews related work. Section 3 describes the methodology. Section 4 reports synthetic and real-data results, and Section 6 concludes.

1.1. Granger Causality

Granger causality (GC) formalizes the predictive influence between time series: $X_1$ is said to Granger-cause $X_2$ if past values of $X_1$ provide statistically significant information for predicting $X_2$ beyond what is already contained in the past of $X_2$ [3]. Considering two stochastic signals modeled as univariate autoregressions $X_i\left(t\right)={\sum }_{k=1}^p{a}_{i,k}{X}_i(t-k)+{ϵ}_i\left(t\right)$, for $i\in \{1,2\}$, with the corresponding bivariate models

$$\begin{array}{cc}\hfill X_1\left(t\right)& ={\displaystyle \sum _{i=1}^p}{a}_i{X}_1(t-i)+{\displaystyle \sum _{i=1}^p}{\beta }_i{X}_2(t-i)+{ϵ}_{1\mid 2}\left(t\right),\hfill \\ \hfill X_2\left(t\right)& ={\displaystyle \sum _{i=1}^p}{b}_i{X}_2(t-i)+{\displaystyle \sum _{i=1}^p}{\gamma }_i{X}_1(t-i)+{ϵ}_{2\mid 1}\left(t\right).\hfill \end{array}$$

(1)

$X_2$ Granger-causes $X_1$ if the cross-coefficients ${\beta }_i$ are jointly different from zero. The null hypothesis $H_0:{\beta }_1=\dots ={\beta }_p=0$ is tested through a nested F-test, $F=((RS{S}_r-RS{S}_u)/m)/(RS{S}_u/(T-k))$, where $RS{S}_r$ and $RS{S}_u$ are the residual sums of squares of the restricted and unrestricted models, $m=p$ is the number of restrictions, T the sample size, and k the number of parameters of the unrestricted model. We adopt $\alpha =0.05$ and select the lag order p via the Bayesian Information Criterion (BIC) [17].

1.2. Vector Visibility Graph

The Vector Visibility Graph extends the classical visibility graph to multivariate series, mapping them onto a directed complex network [13,14]. Given a series $X_t={\left\{x_{i,t}\right\}}_{i=1}^m$, $t=1,\dots ,N$, each vector ${\mathbf{X}}_t=[x_{1,t},\dots ,x_{m,t}]$ is a node. The projection of ${\mathbf{X}}_b$ onto ${\mathbf{X}}_a$ ($a<b$) is $\parallel {\mathbf{X}}_a^b\parallel ={\sum }_i{x}_{i,a}{x}_{i,b}/{\sum }_i{x}_{i,a}^2$. The vectors ${\mathbf{X}}_a$ and ${\mathbf{X}}_b$ have directed visibility if, for every c such that $a<c<b$,

$$∥{\mathbf{X}}_a^c∥<∥{\mathbf{X}}_a^b∥+\left(\parallel {\mathbf{X}}_a\parallel -∥{\mathbf{X}}_a^b∥\right){\displaystyle \frac{t_b-t_c}{t_b-t_a}}.$$

(2)

Geometrically, ${\mathbf{X}}_b$ is visible from ${\mathbf{X}}_a$ if no intermediate vector dominates the line of sight, where the obstruction is measured by the scalar projection of each vector onto the direction of ${\mathbf{X}}_a$. Time instants with similar and persistent causal structures tend to be mutually visible, forming hubs of high connectivity. All topological metrics below are computed on the underlying undirected version of the VVG, since the objective is to characterize the intensity and structure of causal coupling rather than its direction.

1.3. Betweenness Centrality

Betweenness centrality [15] quantifies the importance of a node v along the shortest paths of the network: $C_B\left(v\right)={\sum }_{s\ne v\ne t}{\sigma }_{st}\left(v\right)/{\sigma }_{st}$, with ${\sigma }_{st}$ the number of shortest paths between s and t and ${\sigma }_{st}\left(v\right)$ the number passing through v. For an undirected network with N nodes, the reference value adopted in this work is $C_B\left(\mathrm{hub}\right)=(N-1)(N-2)/2$, the centrality of the hub in a canonical star topology.

1.4. Network Laplacian Spectral Descriptor (NetLSD)

NetLSD [16] is a permutation-invariant, size-invariant, and scale-adaptive spectral descriptor. Given $G=(V,E)$ with $n=\left|V\right|$, the normalized Laplacian $\mathcal{L}=I-D^{-1/2}A{D}^{-1/2}$ admits the spectral decomposition $\mathcal{L}=\mathsf{\Phi }\mathsf{\Lambda }{\mathsf{\Phi }}^{\top }$, with eigenvalues ${\mu }_i\in [0,2]$. The heat equation ${\partial }_{t}u=-\mathcal{L}u$ admits $H_t=e^{-t\mathcal{L}}$, and NetLSD uses the trace

$$h_t=\mathrm{tr}\left(H_t\right)={\displaystyle \sum _{j=1}^n}{e}^{-t{\mu }_j},$$

(3)

sampled at logarithmically spaced $t\in [{10}^{-2},{10}^2]$. We normalize by the empty graph signature, $\tilde{h}\left(G\right)=h\left(G\right)/h\left({\overline{K}}_n\right)$ with $h_t\left({\overline{K}}_n\right)=1/n$, and use the $L_2$ distance $d_{\mathrm{NetLSD}}(G_1,G_2)={\parallel \tilde{h}\left(G_1\right)-\tilde{h}\left(G_2\right)\parallel }_2$ between normalized signatures. For each pair of coupling values $({\lambda }_i,{\lambda }_j)$ and autoregressive coefficient $\varphi$, we compute the mean over $R=100$ independent realizations,

$$\overline{d}({\lambda }_i,{\lambda }_j;\varphi )={\displaystyle \frac{1}{R}}{\displaystyle \sum _{r=1}^R}{∥\tilde{h}\left(G_{{\lambda }_i,\varphi }^{\left(r\right)}\right)-\tilde{h}\left(G_{{\lambda }_j,\varphi }^{\left(r\right)}\right)∥}_2.$$

(4)

Diagonal elements (${\lambda }_i={\lambda }_j$) quantify intra-regime variability (spectral reproducibility); off-diagonal elements (${\lambda }_i\ne {\lambda }_j$) quantify inter-regime discriminability. The heatmaps are normalized row-wise.

1.5. Deterministic Tourist Walk (DTW)

The DTW is a partially self-avoiding deterministic walk [18,19]. An agent departs from each node $v_0\in V$ and moves following the maximum-degree rule

$$v_{k+1}=arg\underset{u\in \mathcal{N}\left(v_k\right)\setminus {\mathcal{V}}_k}{max}deg\left(u\right),$$

(5)

with ${\mathcal{V}}_k=\{v_0,\dots ,v_{k-1}\}$. Ties are broken by node index; if all neighbors have already been visited, the tourist revisits the highest-degree neighbor, generating a limit cycle. Each trajectory splits into a transient phase and an attractor phase, characterized by the means and standard deviations of the transient $\overline{t},{\sigma }_t$ and period $\overline{p},{\sigma }_p$, the number of unique attractors $N_{\mathcal{A}}$, and the dominance of the largest basin

$$D_{max}\left(G\right)={\displaystyle \frac{{max}_{\mathcal{A}}\left|\mathrm{Basin}\left(\mathcal{A}\right)\right|}{N}}.$$

(6)

High $D_{max}$ indicates structural coherence; low values indicate dynamical fragmentation.

2. Related Work

Complex networks offer a powerful framework to analyze time series and detect regime changes, characterize nonlinear dependencies, and quantify information flow [20,21,22,23]. Methods for time-varying causal networks include kernel reweighted group lasso for dynamic VAR models [24], dynamic graph structure learning [25], and spectral co-clustering on evolving directed graphs. Visibility graphs [13,14,20,26] and their multivariate and vector generalizations [6,7,27] bridge time series and network analysis, while NetLSD currently represents the state of the art in compact, invariant graph representation. Despite this variety, there is a lack of compact, interpretable topological descriptors capable of characterizing time-varying causal networks in a multimodal fashion. The proposed framework fills this gap by integrating sliding-window Granger estimation with the VVG and validating it through three complementary descriptors. The DTW, conceptually simple and discriminative, is employed here as a long-term dynamical descriptor, complementing NetLSD (spectral) and betweenness centrality (local). Table 1 summarizes the differentiation: the central axis is not the causal estimator (treated as an interchangeable input) but the multimodal topological analysis of the resulting networks.

3. Materials and Methods

Figure 1 shows the pipeline, structured in five stages: generation of time series, extraction of causality matrices via sliding windows, vectorization of the matrices, construction of the VVG, and topological analysis through the three descriptors.

3.1. Synthetic Time Series

3.1.1. Scenario A (Coupled AR(1) Chain)

Four series $(A,B,C,D)$ of length $N=5,000$ are generated by an AR(1) process with linear lagged coupling:

$$\begin{array}{cccc}\hfill A\left(t\right)& =\varphi A(t-1)+{\eta }_1\left(t\right),\hfill & \hfill B\left(t\right)& =\varphi B(t-1)+\lambda A(t-5)+{\eta }_2\left(t\right),\hfill \\ \hfill C\left(t\right)& =\varphi C(t-1)+\lambda B(t-5)+{\eta }_3\left(t\right),\hfill & \hfill D\left(t\right)& =\varphi D(t-1)+\lambda C(t-5)+{\eta }_4\left(t\right),\hfill \end{array}$$

(7)

with ${\eta }_i\sim \mathcal{N}(0,1)$ i.i.d. We use $\varphi \in \{0.0,0.2,0.6,0.8\}$ and $\lambda \in \{0.0,0.1,\dots ,1.0,1.5,2.0,\dots ,5.5\}$ (20 values), yielding 80 combinations with $R=100$ independent simulations each (8,000 realizations). Initial values for $t\in \{0,\dots ,4\}$ are sampled from $\mathcal{N}(0,1)$. The coefficient $\varphi$ controls intrinsic memory ($\varphi =0$ produces white noise, large $\varphi$ introduces strong autocorrelation), and $\lambda$ regulates the coupling intensity along the chain. Figure 2 illustrates, for $\varphi =0.6$, the cases $\lambda =0$ and $\lambda =1$. We also consider an abrupt change at $t=2,500$: (i) $\lambda :0\to 1$ and (ii) $\lambda :1\to 0$ (Figure 2, right panels).

3.1.2. Scenario B (Coupled Chaotic Logistic Maps)

A bivariate experiment with coupled chaotic oscillators (the coupling constant is denoted $\kappa$ to avoid notational collision):

$$\begin{array}{cc}\hfill f_r\left(x\right)& =rx(1-x),\hfill \\ \hfill A_{t+1}& =(1-\kappa )f_r\left(A_t\right)+\kappa f_r\left(B_t\right)+{\xi }_t^A,B_{t+1}=(1-\kappa )f_r\left(B_t\right)+\kappa f_r\left(A_t\right)+{\xi }_t^B,\hfill \end{array}$$

(8)

with ${\xi }_t^A,{\xi }_t^B\sim \mathcal{N}(0,{\sigma }^2)$, $\sigma =0.1$, $A_0=0.123,B_0=0.654$, and $r=4$. For $r=4$, in the deterministic case, the Lyapunov exponent of the logistic map is ${\mathsf{\Lambda }}_L=ln2$. Transverse stability of the synchronized solution requires ${\mathsf{\Lambda }}_L+ln|1-2\kappa |<0$, so that $\kappa \in ({\kappa }_c,1-{\kappa }_c)$, with ${\kappa }_c=(1-e^{-{\mathsf{\Lambda }}_L})/2=1/4$, i.e., $\kappa \in (0.25,0.75)$ [28,29]. In the presence of independent noise, synchronization becomes approximate, characterized by the reduction of the transverse error within this coupling range.

3.2. Sliding-Window Causal Vectors and VVG Construction

Within each window $[t,t+L]$, the Granger causality p-value matrix is

$$C_{ij}\left(t\right)=\left\{\begin{array}{cc}p{\text{-}\mathrm{value}}_F(T_i\to T_j;p),\hfill & i\ne j,\hfill \\ 1,\hfill & i=j,\hfill \end{array}\right.$$

(9)

and is linearized into $\mathbf{v}\left(t\right)=\mathrm{vec}\left(C\left(t\right)\right)\in {\mathbb{R}}^{m^2}$. We tested $L\in \{50,100,200,400,800,1,600,3,200,4,000\}$ and report $L=50$ for Scenario A (maximum sensitivity to $\lambda$) and $L=400$ for Scenario B. The sequence ${\left\{\mathbf{v}\left(t\right)\right\}}_{t=L}^N$ defines the nodes of the VVG. A directed edge $i\to j$ is created when $\mathbf{v}\left(i\right),\mathbf{v}\left(j\right)$ satisfy Equation (2); the topological analysis uses the underlying undirected version.

3.3. Topological Descriptors and Validation

The aggregate measure of structural intensity is the normalized maximum betweenness centrality

$${\tilde{C}}_B^{max}(\lambda ;\varphi )={\displaystyle \frac{C_B^{max}(\lambda ;\varphi )}{(N_v-1)(N_v-2)/2}},$$

(10)

with $N_v=N-L$ nodes (for $L=50$, $N=5,000$: $N_v=4,950$, denominator $12,243,826$). For NetLSD, we compute Equation (4) for all triples $(\varphi ,{\lambda }_i,{\lambda }_j)$ with $R=100$. For the DTW, a walk is initiated from every node following Equation (5), and the aggregate metrics $\overline{t},{\sigma }_t,\overline{p},{\sigma }_p,N_{\mathcal{A}},D_{max}$ are computed and averaged over the 100 realizations. All results are reported as means with standard error.

3.3.1. Computational Complexity

Let $N_v$ be the number of VVG nodes and E the number of edges. The Granger stage costs $O(N_v·m^2·L{p}^2)$; VVG construction is $O\left(N_v^2\right)$ for sparse graphs (worst case $O\left(N_v^3\right)$); Brandes’s betweenness algorithm [15] is $O\left(N_{v}E\right)$; the NetLSD spectral decomposition is $O\left(N_v^3\right)$ in general, $O\left(kE\right)$ with Lanczos for sparse graphs; the DTW is $O\left(N_v^2\overline{d}\right)$, with $\overline{d}$ the mean degree. For Scenario A ($N=5,000$, $L=50$, $m=4$, $p\le 5$), one realization runs in 2–5 min on a single CPU core (Python); parallelization over the $R=100$ realizations scales linearly with the number of cores.

4. Results

4.1. Degree Distribution in Causal Networks

4.1.1. Scenario A

Figure 3 compares the VVG degree distributions with windows $L=50$ and $L=100$ for $\varphi =0$, contrasting $\lambda =0$ and $\lambda =1$. With $L=50$, the topological difference is immediate: the $\lambda =0$ distribution is concentrated, with degrees typically below 200, while $\lambda =1$ exhibits tails extending up to degrees of the order of the network size (∼5,000), indicating the emergence of giant hubs. In the VVG, the degree of a node t measures how many other instants $t^{\prime }$ can see the causal vector $\mathbf{v}\left(t\right)$ without geometric obstruction. For $\lambda =0$, vectors are dominated by noise and the visibility geometry is essentially random, consistent with VVGs from white-noise series [30]. For $\lambda =1$, the chain $A\to B\to C\to D$ imprints persistent signatures on $\mathbf{v}\left(t\right)$, and instants with pronounced causal structure become geometric reference points. When L is increased to 100, the two distributions become almost indistinguishable, signaling a loss of sensitivity. The window controls temporal resolution: small L preserves the variability of $\mathbf{v}\left(t\right)$, while large L acts as a temporal low-pass filter, homogenizing $\mathbf{v}\left(t\right)$ and destroying topological discrimination. The choice $L=50$ is the empirical balance for Scenario A.

The circular visualization in Figure 6 (top) complements this analysis. Communities were identified by the Girvan–Newman algorithm [31]; in all analyzed cases the algorithm detected exactly two communities. For $\lambda =0$, communities are diffuse and weakly delimited; for $\lambda =1$, densely connected communities emerge with visible hubs near the inter-community boundary, acting as temporal bridges between weakly and strongly coupled segments of the series, consistent with recurrence-network results [20,21].

Figure 4 reveals one of the most expressive findings of this section: under an abrupt $\lambda$ change at $t=2,500$, the degree distribution is clearly bimodal, with one mode at degrees 5–8 (uncoupled segment) and another at 15–20 (coupled segment). This bimodality is absent in series of constant coupling and emerges precisely because the VVG simultaneously integrates causal vectors from two distinct regimes. Comparison with a purely stochastic series of the same duration (Figure 4b) confirms that the second mode is a robust structural signature rather than a fluctuation artifact, and Figure 4c shows that the bimodality disappears for $L=3,000$ due to temporal dilution. Practical implication: simple inspection of the degree distribution of the global VVG can signal a structural change without partitioning the series or specifying the change point a priori.

4.1.2. Scenario B

The coupled logistic maps offer a qualitatively distinct robustness test: deterministic and chaotic dynamics, with a precise theoretical synchronization threshold ${\kappa }_c=1/4$ [28,29]. With $L=400$ (Figure 5), the contrast between the two regimes is even sharper than in Scenario A. For $\kappa =0.20$ (outside the synchronization window), the distribution follows a power law typical of VVGs from non-synchronized chaotic series [30]. For $\kappa =0.70$ (inside the synchronization window), the distribution becomes bimodal with a tail up to ∼2,000. When the maps synchronize, their states become essentially identical, the causal vectors $\mathbf{v}\left(t\right)$ converge to a stationary state of high coherence, and highly coherent vectors share nearly total mutual visibility, generating hubs. The bimodality is informative: the first mode reflects transient instants in which synchronization is partially broken by the noise $\sigma =0.1$, the second mode corresponds to instants of full synchronization, and the separation between them encodes the fraction of time spent in full synchronization. The circular visualizations (Figure 6, bottom) confirm this picture: for $\kappa =0.20$, a diffuse structure without hubs; for $\kappa =0.70$, high-connectivity hubs concentrated in a restricted region of the circle, forming a dense core of fully synchronized instants. Taken together, Scenarios A and B demonstrate that the degree distribution of the causal VVG is a robust topological discriminator sensitive both to the intensity of the coupling and to the underlying dynamical mechanism (linear vs. chaotic, stochastic vs. deterministic).

Figure 5. VVG degree distribution with $L=400$ for coupled logistic maps. Blue: $\kappa =0.70$ (synchronized); orange: $\kappa =0.20$ (non-synchronized). The bimodal distribution with a tail up to ∼2,000 in the synchronized regime contrasts with the power-law distribution of the non-synchronized regime.

Figure 5. VVG degree distribution with $L=400$ for coupled logistic maps. Blue: $\kappa =0.70$ (synchronized); orange: $\kappa =0.20$ (non-synchronized). The bimodal distribution with a tail up to ∼2,000 in the synchronized regime contrasts with the power-law distribution of the non-synchronized regime.

Figure 6. Circular visualization of the VVG. Two communities detected by the Girvan–Newman algorithm [31] in all cases. (a) Scenario A with $\lambda =0$. (b) Scenario A with $\lambda =1$. (c) Scenario B with $\kappa =0.20$. (d) Scenario B with $\kappa =0.70$. Coupled/synchronized regimes (right column) show densely connected communities with visible inter-community hubs, absent in the uncoupled/non-synchronized regimes (left column).

Figure 6. Circular visualization of the VVG. Two communities detected by the Girvan–Newman algorithm [31] in all cases. (a) Scenario A with $\lambda =0$. (b) Scenario A with $\lambda =1$. (c) Scenario B with $\kappa =0.20$. (d) Scenario B with $\kappa =0.70$. Coupled/synchronized regimes (right column) show densely connected communities with visible inter-community hubs, absent in the uncoupled/non-synchronized regimes (left column).

4.2. Structural Characterization via Betweenness Centrality

Figure 7 shows the normalized mean maximum betweenness centrality (Equation (10)) as a function of $\lambda$ for the four values of $\varphi$. For each $\varphi$, the 100 trials are fitted by a saturated exponential

$${\tilde{C}}_B^{max}(\lambda ;\varphi )=a\left(\varphi \right)-b\left(\varphi \right)e^{-c\left(\varphi \right)\lambda },$$

(11)

whose parameters and goodness-of-fit are reported in Table 2. The centrality grows from ${\tilde{C}}_B^{max}\approx 0.04$ at $\lambda =0$ to a saturation plateau $a\left(\varphi \right)$, reflecting the competition between stochastic disorder (dominant for small $\lambda$, random visibility geometry, no dominant mediating nodes) and the hierarchical organization induced by the chain $A\to B\to C\to D$ (dominant for large $\lambda$, instants of intense causal transmission become hubs).

All four coefficients yield $R^2>0.96$ and ${\chi }_{\nu }^2\in [0.88,1.70]$ with $p\left({\chi }^2\right)>0.09$, confirming that the fits are statistically compatible with the experimental errors. Comparison with a logistic sigmoid of the same number of parameters [32] shows that the exponential is preferable for $\varphi \in \{0.0,0.6\}$ ($\Delta \mathrm{AIC}=-4.38$ and $-4.77$; strong evidence) and statistically indistinguishable for $\varphi \in \{0.2,0.8\}$; in the latter cases, the saturated exponential is adopted for the interpretability of its parameters: $a\left(\varphi \right)$ is the topological saturation level, $b\left(\varphi \right)$ the total transition amplitude, and $c\left(\varphi \right)$ the rate of emergence of the causal signal over the noise. Both $a\left(\varphi \right)$ and $b\left(\varphi \right)$ decrease monotonically with $\varphi$ (autoregressive memory distributes the causal signal more uniformly, attenuating topological concentration), while $c\left(\varphi \right)$ grows from 0.643 ($\varphi =0$) to 1.840 ($\varphi =0.8$): autocorrelation increases the signal-to-noise ratio of the Granger estimator at fixed $\lambda$, accelerating the transition although the saturation level reached is lower. The standard error of c for $\varphi =0.8$ is substantially larger (${\sigma }_c/c\approx 27\%$) because saturation occurs rapidly and few observational points precede the plateau; the claim of monotonic growth of $c\left(\varphi \right)$ is therefore robust but with reduced certainty for the extreme value $\varphi =0.8$.

4.3. Spectral Validation: NetLSD

Figure 8 shows the heatmaps of the normalized mean $L_2$ distance between NetLSD signatures (Equation (4)) for $L=50$ and the four values of $\varphi$. For all $\varphi$, the heatmaps display a bipartite structure with a transition band that splits the parameter space into two qualitatively distinct blocks.

The red block corresponds to pairs $({\lambda }_i,{\lambda }_j)$ with significant coupling difference, such as $({\lambda }_0=0.0,{\lambda }_{20}=5.5)$: VVGs generated in such distinct regimes exhibit markedly different spectral signatures, resulting in high $L_2$ distances and reflecting the sensitivity of NetLSD to the coupling constant. The blue block, in turn, corresponds to pairs $({\lambda }_i,{\lambda }_j)$ whose coupling values are close to each other: since VVGs generated in similar regimes share the same dominant topological structure, their NetLSD signatures converge and the residual $L_2$ distances (∼0.09–0.22) are small, representing intrinsic variability between realizations rather than genuine structural difference. In particular, diagonal elements (${\lambda }_i={\lambda }_j$) quantify intra-regime spectral reproducibility. The transition band corresponds to the interval of $\lambda$ over which the spectrum of the normalized Laplacian undergoes its largest reorganization, captured by the heat trace at multiple scales.

The threshold ${\lambda }^{*}$ shifts monotonically: ${\lambda }^{*}\approx 0.9$–$1.0$ for $\varphi =0.0$ and $0.2$, broadening for $\varphi =0.6$ and shifting to ${\lambda }^{*}\approx 0.5$–$0.7$ for $\varphi =0.8$. The residual distances in the blue block grow from 0.09–0.13 (low $\varphi$) to 0.13–0.22 ($\varphi =0.8$). All heatmaps are approximately symmetric, consistent with the symmetry of the $L_2$ metric. The coincidence between this spectral threshold and the point at which betweenness centrality begins its systematic growth (Section 4.2) is a non-trivial coherence result and provides internal cross-validation of the method.

4.4. Dynamical Analysis via DTW

Figure 9a shows the mean transient $\overline{t}$ as a function of $\lambda$. The pattern is non-monotonic: for all $\varphi$, $\overline{t}$ decreases from $\lambda =0$, reaches a minimum in the range $\lambda \approx 0.5$–$0.7$, and recovers toward a plateau comparable to the baseline. For $\varphi =0$, $\overline{t}$ drops from ∼170 at $\lambda =0$ to ${\overline{t}}_{min}\approx 105$ at $\lambda \approx 0.6$, recovering to ∼150–175 for $\lambda \ge 1.5$. For $\varphi =0.8$ the same pattern holds with globally suppressed values (minimum ∼70, plateau ∼90). A Wilcoxon test (Holm correction [33]) confirms the valley as significant only for $\varphi =0.2$, $\lambda =0.5$ ($W=3,560$, $p_{\mathrm{corr}}=0.003$); for the other $\varphi$, the pattern is consistent in the means but indistinguishable from noise at the individual-realization level, reflecting limited statistical power rather than absence of the phenomenon. The mechanistic interpretation is robust: at $\lambda =0$, the absence of hubs forces long walks; at moderate $\lambda$, incipient hubs quickly capture the tourist but are not yet globally connected, yielding short cycles and minimum transients; at $\lambda \gtrsim 1$, giant hubs dominate the topology and the tourist circulates in cycles involving them, raising $\overline{t}$ back to the plateau.

Figure 9b shows the monotonic collapse of $N_{\mathcal{A}}$ with two phases: rapid decay for $\lambda \lesssim 1.5$ and a plateau for $\lambda \gtrsim 1.5$. From $N_{\mathcal{A}}\approx 860$–890 at $\lambda =0$ (high dynamical fragmentation), progressively more connected hubs attract trajectories into a smaller set of limit cycles. The Kruskal–Wallis test on plateau points ($\lambda \ge 3.0$) strongly rejects equality between groups ($H=522.4$, $p=6.7\times {10}^{-113}$), and pairwise Mann–Whitney comparisons with Holm correction (Table 3) confirm all six pairs as significantly distinct. This formally establishes that autoregressive memory preserves dynamical diversity even under strong causal coupling: each $\varphi$ converges to a structurally distinct plateau, around 240–290 attractors for $\varphi \in \{0.0;0.2\}$ and 300–400 for $\varphi \in \{0.6;0.8\}$—a result consistent with the role of autocorrelation in maintaining dynamically diverse trajectories even when the network topology is dominated by high-connectivity hubs.

Figure 9c shows the sigmoidal behavior of $D_{max}$ with three phases: low plateau ($\lambda <0.5$, no dominant attractor), abrupt growth around $\lambda \approx 0.9$–$1.0$, and a new plateau ($\lambda >2$) whose level depends strongly on $\varphi$. Fitting a logistic $L/[1+e^{-k(\lambda -{\lambda }^{*})}]$ yields ${\lambda }^{*}=1.14$ (95% CI $[0.98;1.30]$) for $\varphi =0.0$, $1.22$ ($[1.11;1.33]$) for $\varphi =0.2$, $0.87$ ($[0.73;1.00]$) for $\varphi =0.6$, and $0.06$ ($[-0.05;0.17]$) for $\varphi =0.8$. The CIs overlap and are compatible with ${\lambda }^{*}\approx 1$ for $\varphi \in \{0.0;0.2;0.6\}$, formally confirming the convergence already identified by betweenness centrality and NetLSD. For $\varphi =0.8$, the fit returns ${\lambda }^{*}=0.06$ with a CI including zero, indicating that the transition is already complete before the smallest $\lambda$ in the grid and the sigmoidal fit loses resolution in this regime. The sigmoidal character of $D_{max}\left(\lambda \right)$ is a classical phase-transition signature [34], with the order parameter remaining small in the disordered phase and saturating in the ordered phase; its abrupt growth near $\lambda \approx 1$ is the dynamical counterpart of the centrality saturation.

Figure 9d integrates $\overline{t}$ and $D_{max}$ in a single scatter plot, showing a quasi-orthogonal separation of $\varphi$ (cluster position) and $\lambda$ (displacement within cluster), useful for identifying the dynamical regime of an unknown real system from its position in this plane.

4.4.1. Internal Cross-Validation

The three descriptors converge into a coherent picture. First, the threshold ${\lambda }^{*}\approx 0.9$–$1.0$ identified by NetLSD coincides with (a) the onset of systematic centrality growth ($c\left(\varphi \right)>0$ for all $\varphi$, $R^2>0.96$) and (b) the sigmoidal inflection of $D_{max}$, whose 95% CIs overlap among $\varphi \in \{0.0;0.2;0.6\}$ and contain $\lambda =1$. Shortest paths, Laplacian spectrum, and attractor basins formally point to the same critical point. Second, the U-shaped minimum of $\overline{t}$ at $\lambda \approx 0.5$–$0.7$ lies slightly below ${\lambda }^{*}$, consistent with the interpretation that this region precedes topological consolidation. Third, the effect of $\varphi$ is quantitatively consistent and statistically confirmed across descriptors: monotonic reduction of $a\left(\varphi \right)$ and $b\left(\varphi \right)$; leftward shift of the NetLSD boundary; and significant distinction of all $N_{\mathcal{A}}$ plateau pairs. This triple convergence (local, spectral, dynamical) constitutes the internal cross-validation of the framework, the most rigorous form of robustness evaluation available in the absence of established external baselines for the specific task of topological characterization.

4.5. Case Study: Meteorological Variables and Rice Prices

To illustrate applicability to real data, we conducted a case study on meteorological time series and rice prices in Southern Brazil during 2015–2025, using INMET data (https://portal.inmet.gov.br/dadoshistoricos, accessed on 1 March 2026) combined with CEPEA rice price series (https://www.cepea.org.br/br/consultas-ao-banco-de-dados-do-site.aspx, accessed on 1 March 2026) for Uruguaiana (RS), the country’s main rice-producing hub. The goal is to demonstrate the discriminative power of the topological characterization, not to claim discovery of causal drivers of rice prices.

The dataset comprises a multivariate series of 2,000 points. From the 17 available meteorological variables, all $\left({\displaystyle \genfrac{}{}{0pt}{}{17}{3}}\right)=680$ triplets were formed, each concatenated with the rice price series to produce a multivariate time series of dimension $m=4$, matching Scenario A. The pipeline was applied to each of the 680 series independently. Combinations were ranked by ${\tilde{C}}_B^{max}$ and the top 15 (Top 15) and bottom 15 (Bottom 15) were retained for comparative analysis. The separation between groups is assessed by the DTW metrics ($\overline{t}$ and $D_{max}$), which are mathematically independent of betweenness centrality and did not enter the selection criterion; convergence of the three metrics is therefore non-trivial evidence.

Figure 10a shows the centrality profiles. The separation is sharp: Top 15 curves reach ${\tilde{C}}_B^{max}\sim 0.20$, while Bottom 15 stays below ∼0.10. The value 0.20 is consistent with the saturation plateau $a(\varphi =0.0)$ in Table 2, suggesting that the most prominent meteorological combinations induce an effective causal coupling of the order $\lambda \gtrsim 1.5$ in the equivalent generating system. Figure 10b,c report $\overline{t}$ and $D_{max}$ for the two groups. The separation was assessed via two-sided Mann–Whitney U tests: for $\overline{t}$, $U=225$ ($p=2.20\times {10}^{-6}$, $\left|r\right|=1.000$); for $D_{max}$, $U=225$ ($p=3.38\times {10}^{-6}$, $\left|r\right|=1.000$). The value $U=n_1{n}_2=225$ is the maximum possible for $n_1=n_2=15$: every value of the Top 15 exceeds every value of the Bottom 15 with no overlap. Perfect separation in two metrics independent of the ranking criterion establishes that the discrimination is not an artifact of selection but a coherent structural property of the framework on real data.

Table 4 reports the frequency of meteorological variables in the two groups. The asymmetry is marked: Top 15 is dominated by minimum relative humidity, atmospheric pressure (maximum, minimum, and station-level), and wind variables (direction and speed). This set is consistent with the influence of synoptic-scale atmospheric systems, such as cold fronts, extratropical cyclones, and post-frontal anticyclones, which periodically modulate meteorological conditions during the rice growing season (October to March). In this context, atmospheric pressure should be interpreted as an indirect indicator of the passage and reorganization of these systems, while wind and relative humidity capture changes in advection, atmospheric stability, and evaporative demand. The association between high-pressure systems and evapotranspiration, however, depends on the thermal regime, humidity, radiation, and wind; therefore, their effects on irrigation costs and grain quality must be interpreted conditionally [35,36,37].

Bottom 15 is dominated by total precipitation, minimum temperature, and hourly relative humidity. The lower frequency of these variables is compatible with the specificity of flood-irrigated lowland rice in Uruguaiana, whose water availability is strongly mediated by irrigation management, reducing direct dependence on local precipitation compared to rainfed systems. Thus, the topological separation suggests a plausible regional agronomic signature [38], rather than an arbitrary methodological choice. We do not claim causality in the econometric sense; such inference would require explicit control of confounding variables—including international prices, exchange rates, agricultural policies, and management decisions—in addition to dedicated time-series modeling.

5. Discussion

The results presented above support a coherent interpretation of the proposed framework. The convergence of three mathematically independent descriptors (betweenness centrality, NetLSD, and DTW) toward a common critical threshold ${\lambda }^{*}\approx 0.9$–$1.0$ in Scenario A is the central finding: it shows that the topological reorganization of the causal VVG near the coupling-induced transition is not an artifact of a single metric but a genuine structural property captured from local, spectral, and dynamical perspectives simultaneously. The systematic leftward shift of ${\lambda }^{*}$ as $\varphi$ increases admits a physically interpretable explanation: autoregressive memory increases the effective signal-to-noise ratio of the Granger estimator at fixed $\lambda$, anticipating the topological consolidation; the same mechanism explains the monotonic growth of $c\left(\varphi \right)$ in the centrality fit and the reduction of the saturation level $a\left(\varphi \right)$, since memory distributes the causal signal more uniformly across time.

The detection of regime transitions through bimodality in the global VVG degree distribution offers a practical advantage over methods that require a priori specification of change points. Together with the chaotic-map experiments, which confirm consistency with the theoretical synchronization threshold ${\kappa }_c=1/4$, this suggests that the framework is sensitive to the underlying dynamical mechanism rather than to the specific functional form of the coupling. The Uruguaiana case study illustrates that this sensitivity extends to real data: the perfect non-parametric separation between Top 15 and Bottom 15 across three independent metrics shows that the topological discrimination is structural rather than a selection artifact, and the agronomic plausibility of the most frequent variables in each group lends external coherence to the analysis.

Several limitations should be acknowledged. First, the framework treats the causal estimator as an interchangeable input; while Granger causality is convenient and well understood, regimes with strong nonlinearity or non-Gaussian noise may require alternatives such as Transfer Entropy [8] or PCMCI [9], whose interaction with the VVG and the three descriptors deserves dedicated investigation. Second, the choice of the sliding window L is empirically tuned and exhibits a clear bias–variance trade-off: small L preserves topological discrimination but increases variability, while large L acts as a temporal low-pass filter; principled criteria for window selection in nonstationary settings remain an open direction [39]. Third, the real-data analysis is illustrative rather than confirmatory; rigorous causal inference on rice prices would require explicit control of confounders such as international prices, exchange rates, and agricultural policies. Future work includes extending the framework to higher-dimensional series, replacing Granger with nonlinear causal estimators, and exploring the convergence of additional topological descriptors (e.g., persistent homology [40,41] and graph wavelets).

6. Conclusions

We introduced a unified method for the multimodal topological characterization of dynamic causal networks built from multivariate time series, integrating sliding-window Granger estimation, the Vector Visibility Graph, and three complementary topological descriptors (betweenness centrality, NetLSD, DTW) into a coherent and interpretable pipeline. The framework does not compete with consolidated causal discovery methods such as PCMCI, Transfer Entropy, and CCM, but complements them: while these estimators aim to identify which causal relations exist, our pipeline characterizes how these relations organize topologically over time.

The central result is the convergence of three mathematically independent perspectives—shortest-path distribution (centrality), normalized Laplacian spectrum (NetLSD), and attractor-basin distribution (DTW)—to a common critical threshold ${\lambda }^{*}\approx 0.9$–$1.0$ for low $\varphi$, with a systematic leftward shift as $\varphi$ increases. This triple convergence constitutes an internal cross-validation of the framework: in the absence of established external baselines for the specific task of topological characterization of dynamic causal networks, the mutual coherence of three independent metrics is the most rigorous form of robustness evaluation available. Quantitatively, the normalized maximum betweenness centrality is well described by a saturated exponential ${\tilde{C}}_B^{max}(\lambda ;\varphi )=a\left(\varphi \right)-b\left(\varphi \right)e^{-c\left(\varphi \right)\lambda }$, whose parameters encode the competition between autoregressive memory and external coupling in an interpretable manner.

The framework also detects regime transitions passively (without partitioning the series): the abrupt $\lambda$ change at $t=2,500$ produces a bimodal degree distribution whose two modes correspond to the two coupling regimes. Scenario B (coupled chaotic logistic maps) confirms consistency with the theoretical synchronization threshold ${\kappa }_c=1/4$, showing that the framework is sensitive to the underlying dynamical mechanism (linear vs. chaotic, stochastic vs. deterministic). The case study using INMET meteorological data and CEPEA rice prices for Uruguaiana (2015–2025) demonstrated perfect topological separation between Top 15 and Bottom 15 across the three metrics, with variables associated with the Top 15 plausible in light of regional agronomic literature.