The overall execution logic and data flow of the proposed physical computation framework are illustrated in Figure 1. As the framework executes, the microscopic thermal characteristics of local overflow load redistribution and the nonlinear evolution of node degree distributions are captured (overviewed in Figure 2 and Figure 3, with dynamic details described below. The subsequent sections map the simulation outcomes directly to the protocol stages.
Temporal topology evolution via TERGM
Figure 4 visually deconstructs the spatial topology and community structure of the ICT-DE500 core network, highlighting the distribution of high-betweenness nodes targeted in deliberate attack simulations. The goodness-of-fit test confirms that the generated network topology effectively models the temporal evolution of real ecosystems, successfully avoiding gradient explosion or traps of local optimality after the 10,000-step burn-in period. Figure 5 presents the MCMC-MLE parameter convergence diagnostic trajectories and the geodesic distance goodness-of-fit. Figure 5A shows that the three core parameters representing edge density θ₁, mutuality θ₂, and geometrically weighted edgewise shared partnership θ₃ all terminate their large directional drift after crossing the 10,000-step burn-in threshold, with expected means converging and stabilizing near the horizontal baseline. Figure 5B demonstrates that the empirical observations of shortest geodesic distances fall robustly within the overall distribution boundaries of 1,000 independent network realizations. These realizations were extracted using a thinning interval of 10 from the 10,000 sequential Gibbs sampling iterations to minimize autocorrelation, confirming the reliability of the topology generation base. The specific MCMC-MLE parameter estimates, standard errors, and statistical significance for endogenous structural effects across the observation years are detailed in Table 2.
The temporal evolution of macroscopic topological structure is quantified in Figure 6. The network density steadily increased from 0.015 to 0.035 between 2018 and 2023, while the average clustering coefficient rose from 0.22 to 0.37, highlighting a significant density-cluster divergence phenomenon. The density rises most sharply between 2020 and 2021, fluctuating from 0.021 to 0.029, whereas the clustering coefficient reaches a local peak of approximately 0.31 in 2020 and subsequently drops to approximately 0.29 despite rapid density growth in 2021. This divergence reveals the adaptive evolution mechanism under macroeconomic cycle fluctuations, where risk-avoidance grouping in 2020 drives local clustering, while massive new cross-border connections in 2021 temporarily dilute the close-knit community structure.
Physical capacity calibration and load initialization
The bivariate sensitivity analysis in Figure 7 examines the joint effect of physical capacity redundancy and load polarization on the phase transition trajectory of the largest connected component. Across the nine orthogonal combinations of capacity tolerance α and load heterogeneity β, the array of panels shows that increasing α and decreasing β both delay network collapse. Under the polarized load scenario with β = 1.5 and minimal redundancy α = 0.1 in Figure 7A, the critical collapse threshold is approximately fc = 0.08. Raising α to 0.5 in Figure 7C shifts the inflection point rightward to fc ≈ 0.23. Under the balanced load scenario with β = 1.0 and α = 0.1 in Figure 7G, the threshold remains robust at fc ≈ 0.18, and at the optimal combination α = 0.5 and β = 1.0 in Figure 7I, the threshold extends significantly to fc ≈ 0.38. These results demonstrate that load balancing produces a greater marginal gain in resilience than capacity expansion alone.
Nonlinear cascading dynamics under deliberate attacks
As introduced in the framework overview, the microscopic thermal characteristics of local overflow load redistribution following the initial cascade failure are shown in Figure 2, and the nonlinear evolution of node degree distribution across three typical micro-time steps is presented in Figure 3.
System resilience evaluation and threshold identification
The critical transition threshold for global transmission efficiency is located at 12.41% ± 0.63% node removal during a targeted attack. In the context of the 500-node network, this proportion corresponds to the targeted removal of approximately 62 core hub nodes. This threshold marks the efficiency-collapse point (i.e., the onset of the steepest decline in transmission efficiency) rather than a complete topological disconnection. Figure 8 presents the three-dimensional evolution surface of global efficiency E(τ) over removal ratio and peak network load in Figure 8A, and the two-dimensional cross-section with numerical differentiation in Figure 8B. When the removal ratio f is below 0.10, E(τ) remains above 0.8 and the first derivative oscillates in a shallow region. The minimum of the first derivative trajectory identifies the critical transition threshold, with Figure 8B displaying a single-run cross-section at fc = 12.0% that is highly consistent with the statistical mean from 100 independent Monte Carlo simulations.
Simulation configurations and baseline implementations
The proposed model significantly outperforms the static and deep learning baselines under targeted attack scenarios. However, under random failure conditions, the Static BA-ML baseline exhibits a higher survival threshold (49.12%) compared to the TERGM-ML model (46.28%). It is important to note that the comparison with the Static BA-ML model serves as a distinct topological benchmark rather than a strictly controlled ablation, as the Barabási-Albert generative mechanism fundamentally differs from the ERGM framework. Figure 9 displays the raincloud plot of cascade propagation peak velocities across four model architectures. The Static BA-ML baseline shows a median peak velocity of approximately 49.7 nodes per step, with extreme batches approaching 140. The SNA-Cascading and GCN-Attack models have medians of approximately 35.6 and 23.9, respectively. The TERGM-ML model exhibits the strongest convergence with a median of 13.2 nodes per step, almost eliminating extreme collapses exceeding 40. Table 3 summarizes the critical thresholds and statistical significance for all models. These comparisons clarify that while the TERGM-ML framework demonstrates superior structural resilience to targeted asymmetric shocks, the observed performance differences reflect the joint influence of distinct underlying topologies, temporal evolution, and capacity-aware redistribution, rather than being solely attributable to isolated mechanism ablations.
DATA AVAILABILITY:
The raw data used in this study are derived from the Crunchbase global corporate investment and merger database, publicly accessible via the Kaggle platform at https://www.kaggle.com/datasets/justinas/startup-investments. The processed ICT-DE500 subset, consisting of 500 entities with annual edge matrices for the period 2018–2023 and node attribute data including degree and betweenness centrality, along with the TERGM estimation and convergence diagnostic scripts, the cascading failure simulation code with the improved Motter-Lai algorithm, and the complete dependency specifications, have been deposited in a publicly accessible GitHub repository at https://github.com/moonmoon1189/digital-economy-resilience-complex-networks.

Figure 1: Temporal topology evolution and nonlinear cascade physical computation framework. This figure illustrates the overall execution logic and data flow, including topology evolution, physical capacity calibration, nonlinear cascading, and resilience evaluation stages for identifying the critical transition threshold. Please click here to view a larger version of this figure.

Figure 2: Thermal diagram of the microscopic evolution of nonlinear load redistribution and local cascade overload. The figure illustrates the dynamic thermal characteristics of the local redistribution of overflow load after the initial cascade failure from micro-time step 0 to step 5. Please click here to view a larger version of this figure.

Figure 3: Scatter plot of nonlinear evolution of node degree distribution during cascading failure. The figure presents the evolution trajectory of the node degree distribution of the system at three typical micro-time steps (0, 3, 6). Please click here to view a larger version of this figure.

Figure 4: Spatial topology, community structure, and deliberate attack target distribution of the ICT-DE500 core network. The figure visually deconstructs the highly nonlinear macroscopic topology and microscopic node attributes of the core network, highlighting star hubs and vulnerable sources. Please click here to view a larger version of this figure.

Figure 5: Goodness-of-fit test of Markov Chain Monte Carlo parameter diagnostic trajectory and geodesic distance. (A) This panel shows the diagnostic trajectory of MCMC-MLE parameter estimation across iterations, while panel (B) displays the goodness-of-fit test of the shortest geodesic distance. Please click here to view a larger version of this figure.

Figure 6: Joint distribution of macro-topological feature parameters in time series evolution. The figure shows the changing trend of macro-topological parameters, specifically network density and average clustering coefficient, for the digital economy ecosystem between 2018 and 2023. The shaded regions around the trend lines represent the 95% confidence intervals derived from 100 independent Monte Carlo simulations. Please click here to view a larger version of this figure.

Figure 7: Array plot of phase transition decay for bivariate sensitivity to capacity tolerance and load heterogeneity. (A–I) These panels display the phase transition trajectories across various orthogonal combinations of capacity tolerance and load heterogeneity. Please click here to view a larger version of this figure.

Figure 8: Joint profile of global transmission efficiency in three dimensions and cross-sectional critical transition threshold. Panel (A) constructs the three-dimensional spatial evolution of global transmission efficiency, and panel (B) extracts the two-dimensional cross-sectional critical transition threshold (i.e., the efficiency-collapse point) localization using numerical differentiation. Please click here to view a larger version of this figure.

Figure 9: Heterogeneous distribution of cascade propagation peak evolution velocity in raincloud plot. The figure comprehensively displays the heterogeneous probability density distribution of the cascade propagation peak velocities of four models during secondary disaster outbreaks. Please click here to view a larger version of this figure.
| Parameter | Variable & Range | Constraint & Justification |
| Network Scale | N=500 | Fixed size of the extracted ICT-DE500 empirical network subset. |
| Centrality Weight | λ=0.5 | Anchored baseline ensuring equal weighting for degree and betweenness. |
| Load Allocation Index | β∈[1.0,1.5] | Non-linear exponent controlling initial business load heterogeneity. |
| Capacity Tolerance | α∈[1.0,1.5] | System-level redundancy ratio defining the upper bound of physical capacity. |
| MCMC Burn-in Period | >10,000 steps | Mandatory initial iterations discarded to achieve stationary distribution. |
| Gibbs Sampling Iterations | >10,000 steps | Sequential sampling steps to generate time-aligned network topologies. |
Table 1: Core parameter settings for time-series evolution and cascading failure physical simulations. The table defines the core parameters, including network scale, centrality weight, load allocation index, capacity tolerance, and Markov Chain Monte Carlo (MCMC) iterations.
| Observation Year | Edge Density (θ1) Estimate | Edge Density (θ1) Std.Error | Mutuality (θ2) Estimate | Mutuality (θ2) Std.Error | GWESP (θ3) Estimate | GWESP (θ3) Std.Error | Significance |
| 2018 | -6.35 | 0.12 | 1.95 | 0.08 | 1.25 | 0.05 | *** |
| 2019 | -6.15 | 0.11 | 1.98 | 0.09 | 1.35 | 0.06 | *** |
| 2020 | -6.05 | 0.14 | 2.15 | 0.1 | 1.55 | 0.07 | *** |
| 2021 | -5.85 | 0.13 | 1.9 | 0.09 | 1.3 | 0.06 | *** |
| 2022 | -5.75 | 0.12 | 2.05 | 0.08 | 1.42 | 0.05 | *** |
| 2023 | -5.65 | 0.11 | 2.1 | 0.07 | 1.48 | 0.05 | *** |
Table 2: TERGM parameter estimates for endogenous structural effects across observation years (2018–2023). The table reports the Markov Chain Monte Carlo Maximum Likelihood Estimation (MCMC-MLE) parameter estimates, standard errors, and statistical significance for edge density, mutuality, and geometrically weighted edgewise shared partnership across the six annual observation periods. ***p < 0.001. Standard errors are reported next to parameter estimates.
| Model Architecture | Temporal Evolution Mechanism | Dynamic Reallocation Mechanism | Critical Threshold (Targeted Attack) | Critical Threshold (Random Failure) | Statistical Significance (p-value) |
| TERGM-ML | Yes | Yes | 12.41% ± 0.63% | 46.28% ± 1.75% | Reference Baseline |
| Static BA-ML | No | Yes | 7.85% ± 0.42% | 49.12% ± 1.88% | p < 0.001 *** |
| SNA-Cascading | Yes | No | 8.93% ± 0.55% | 37.54% ± 1.42% | p = 0.003 ** |
| GCN-Attack (SOTA Baseline) | Implicit | Implicit | 10.76% ± 0.81% | 43.15% ± 2.05% | p = 0.021 * |
Table 3: Comparison of core mechanism ablation and multimodel architecture resilience boundary. The table details the critical thresholds and statistical test results for system resilience of time-series graphs and improved Motter-Lai (TERGM-ML) and three baseline models under both deliberate attack and random failure scenarios. Values are reported as mean ± standard deviation based on 100 independent Monte Carlo simulations. Targeted attack refers to the sequential removal of nodes based on descending betweenness centrality. Statistical significance evaluates the difference in the Targeted attack threshold between the respective baseline model and the proposed framework using an independent two-sample t-test (*p < 0.05, **p < 0.01, ***p < 0.001).