Reviving a failed network through microscopic interventions

From mass extinction to cell death, complex networked systems often exhibit abrupt dynamic transitions between desirable and undesirable states. These transitions are often caused by topological perturbations (such as node or link removal, or decreasing link strengths). The problem is that reversing the topological damage, namely, retrieving lost nodes or links or reinforcing weakened interactions, does not guarantee spontaneous recovery to the desired functional state. Indeed, many of the relevant systems exhibit a hysteresis phenomenon, remaining in the dysfunctional state, despite reconstructing their damaged topology. To address this challenge, we develop a two-step recovery scheme: first, topological reconstruction to the point where the system can be revived and then dynamic interventions to reignite the system’s lost functionality. By applying this method to a range of nonlinear network dynamics, we identify the recoverable phase of a complex system, a state in which the system can be reignited by microscopic interventions, for instance, controlling just a single node. Mapping the boundaries of this dynamical phase, we obtain guidelines for our two-step recovery. Perturbations and disturbances can bring complex networks into undesirable states in which global functionality is suppressed. Now, a recovery scheme explains how to revive a damaged network by controlling only a small number of nodes.

From mass extinction to cell death, complex networked systems often exhibit abrupt dynamic transitions between desirable and undesirable states. Such transitions are often caused by topological perturbations, such as node or link removal, or decreasing link strengths. The problem is that reversing the topological damage, namely retrieving the lost nodes or links, or reinforcing the weakened interactions, does not guarantee the spontaneous recovery to the desired functional state. Indeed, many of the relevant systems exhibit a hysteresis phenomenon, remaining in the dysfunctional state, despite reconstructing their damaged topology. To address this challenge, we develop a two-step recovery scheme: first -topological reconstruction to the point where the system can be revived, then dynamic interventions, to reignite the system's lost functionality. Applying this method to a range of nonlinear network dynamics, we identify the recoverable phase of a complex system, a state in which the system can be reignited by microscopic interventions, for instance, controlling just a single node. Mapping the boundaries of this dynamical phase, we obtain guidelines for our two-step recovery.
Complex systems, biological, social or technological, often experience perturbations and disturbances, from overload failures in power systems 1-3 to species extinction in ecological networks. [4][5][6] The impact of such perturbations is often subtle, the system exhibits a minor response, but continues to sustain its global functionality. 7,8 However, in extreme cases, a large enough perturbation may lead to a major collapse, with the system abruptly transitioning from a functional to a dysfunctional dynamic state [9][10][11][12][13] (Fig. 1a-d). When such collapse occurs, the naïve instinct is to reverse the damage, retrieve the failed nodes and reconstruct the lost links. This, however, is seldom efficient, as (i) we rarely have access to all system components, 14 limiting our ability to reconstruct the perturbed network; (ii) even if we could reverse the damage, due to hysteresis, in many cases, the system will not spontaneously regain its lost functionality.
To address this challenge, we derive here a two-step recovery process: Step I. Restructuring (Fig. 1e). Retrieving the weighted topology to a point where the system can potentially regain its functionality.
Step II. Reigniting (Fig. 1f). Introducing dynamic interventions to steer the system back to its functional state. Considering the fact that in most practical scenarios we cannot control the majority of the system components, we design our reigniting around microinterventions, i.e. controlling just a small number of components. To achieve this we uncover the recoverable phase, a dynamic state in which the system can be driven towards functionality by controlling just a single node.
The challenge of irreversible collapse Consider a complex system of N components (nodes), interacting via the adjacency matrix A, a sparse, potentially directed random network with an arbitrary degree distribution P (k in , k out ) (see Supplementary Section 1). Each node is assigned an activity x i (t), whose meaning depends on context -e.g., a species abundance in a microbial network or a gene's expression level in a biological setting. We then track the evolution of x i (t) following [15][16][17] in which the interaction dynamics is characterized by three potentially nonlinear functions M 0 (x), M 1 (x) and M 2 (x). The first function, M 0 (x i ) captures node i's self-dynamics, describing mechanisms such as protein degradation 18 (cellular), individual recovery 19 (epidemic) or birth/death processes 20 (population dynamics). The product M 1 (x i )M 2 (x j ) describes the i, j interaction mechanism, e.g., genetic activation, 21 infection 19 or symbiosis. 22 The strength of the i, j interaction is governed by W ij , a random weight extracted from the density function P (w), whose average we denote by ω = ∞ 0 wP (w) dw. In the context of recoverability, we seek to revive the activity of all nodes by activating a selected set of nodes, hence we focus on cooperative interactions, in which nodes positively contribute to each others activity. This is expressed in Eq. (1) through W ij M 1 (x i )M 2 (x j ) ≥ 0 (Supplementary Section 1). Later, in our discussion of the microbiome recoverability, we relax this condition and examine the impact of mixed-sign interactions.
Setting the derivative on the l.h.s. of (1) to zero, we obtain the system's fixed-points, x α = (x α,1 , . . . , x α,N ) , which, if dynamically stable, represent different states, desirable or undesirable, in which the system can naturally reside. Transitions between these states often result from perturbations to A or W , such as node/link deletion or reduction in link weight. When this occurs, it is difficult to reverse the unwanted transition. This is because the system often avoids spontaneous recovery, even if we retrieve the lost nodes, links or weights. To illustrate this difficulty we refer to a concrete example below. Example 1. Cellular dynamics (Fig. 3). As our first example we consider the Michaelis-Menten model for gene-regulation, capturing activation interactions between genes (Fig. 3a). Here M 0 (x i ) = −Bx a i , M 1 (x i ) = 1 and M 2 (x j ) = x h j /(1 + x h j ), a switch-like function that saturates to M 2 (x j ) → 1 for large x j , representing the process of activation (see discussion in Supplementary Section 3.1).
For sufficiently connected A or large average weight ω the system exhibits an active state x 1 , in which all x 1,i > 0 -capturing a living cell. However, perturbations to A, W , such as link/node removal or weight loss can cause a sharp transition to the inactive state x 0 = (0, . . . , 0) , i.e. cell death. To track this transition systematically we measured the average activityx α = (1/N ) N i=1 x α,i , which followsx α > 0 for α = 1 andx α = 0 for α = 0. As we subject the system to increasing levels of stress -here, reducing all weights W by a factor 0 ≤ q ≤ 1, we observe a sudden transition at q = q c from x 1 (Fig. 3b, green) to x 0 (red). Next, we attempt to revive the node activities by retrieving the lost weights, finding that the system fails to recover. The reason is that while x 1 is only stable under q ≤ q c , x 0 is always stable -both below and above this threshold. This leads to a hysteresis phenomenon, where the system remains inactive despite the reversal of the perturbation. Example 1, above, while representing a specific scenario, illustrates the family of challenges that we tackle here: system's with irreversible transitions, driven by perturbations to their weighted topology. To revive such systems, we must not just restructure their lost topology, but also dynamically reignite them by exerting external control over the activities x i (t).

Recoverability
The most natural way to reignite the system is to drive all activities x i towards an initial condition from which the system naturally recovers to the desired x 1 . Namely, we must steer the system into x 1 's basin of attraction which comprises all initial conditions x(t = 0) from which Eq. (1) converges to x(t → ∞) = x 1 (Fig. 1h). The problem is that such level of control over the dynamics of all nodes is seldom attainable, hence we seek to recover the system's functionality by driving just a microscopic fraction f → 0 of forced nodes.
To achieve this, we consider the limit f ∼ 1/N , in which case our reigniting is attempted through, typically, a single, randomly selected source node s, whose activity we force to follow x s (t) = φ(t). In many practical applications our ability to exert such control, even on a single node, is restrictive, limiting the potential forms of φ(t). Hence, below we reignite (1) using an extremely simple input, φ(t) = Const. Other practically accessible forms of φ(t) are further considered in Supplementary Section 4. As simple as it is, the constant forcing itself is also constrained, as our forcing capacity is, often, bounded by φ(t) ≤ ∆. Hence, below we seek the conditions where such restricted interventions -controlling just one node and with a forcing bounded by ∆ -can push the remaining nodes into the desired B 1 .
During our intervention, the remaining N − 1 nodes continue to follow the natural system dynamics, i.e. Eq. (1), as they respond to the s-forcing. In technical terms, the failed state of the system, x 0 , captures Eq. (1)'s initial condition, and the forced node imposes a strict boundary condition at s. In a recoverable system, after some time, the activities will enter B 1 , at which point we can cease our external control and allow the system to naturally transition to x 1 , following its internal dynamics. If, however, the system is non-recoverable, such single-node reigniting is insufficient, the system remains at B 0 , and once we lift our forcing, it relaxes back to x 0 , a failed reigniting.
To analytically track the system's response to our forcing at s, we divide the rest of the network into shells K s (l) = {j | L sj = l}, comprising all nodes located at distance l from s (Fig. 2a). In this notation K s (0) = {s}, K s (1) is the group of s's nearest outgoing neighbors, K s (2) its next neighbors, and so on. Then, starting with x s (t) = ∆, we track the average activity of nodes in K s (l), via x s (l, t) = 1 |K s (l)| i∈Ks(l) where |K s (l)| represents the number of nodes in K s (l). The shells adjacent to the source (small l) will be forced to respond to s's activation ∆. The distant shells, however, at l → ∞, may be unaffected, and therefore still within B 0 . Under these conditions, upon termination of our ∆-forcing, all shells retreat back to x 0 , rendering the system unrecoverable. Successful reigniting, therefore, requires that capturing a state in which the forcing at s was able to penetrate the network and affect the activity of even the most distant nodes at K s (l → ∞). This represents a recoverable system that will naturally revert to x 1 once the forcing ∆ is deactivated. 23 To obtain the final shell states x s (l) = x s (l, t → ∞), we use the fact that, despite its potentially broad degree/weight distribution, our network is otherwise wired, and assigned link weights, at random (Supplementary Section 1). Therefore, (i) A features insignificant degree-correlations, and hence the nodes in K s (l) are statistically similar to those in K s (l ), (for l, l ≥ 1); (ii) A is locally tree-like, and therefore, asymptotically, there are almost no short-range loops surrounding the source s. Below we relax both approximations, when testing our method, numerically, against empirically constructed networks, which, indeed, feature both loops and measurable degree-correlations (Supplementary Section 5.5). However, to advance analytically, we use approximations (i) and (ii) to translate Eq. (1) into a direct set of equations for x s (l). We arrive the recurrence relation (Supplementary Section 2) where R(x) = −M 0 (x)/M 1 (x) and R −1 (x) is its inverse function. In (6) the parameterx 0 represents the mean activity of nodes in K s (l > 1) under the failed state x 0 (Supplementary Section 2.4). The remaining parameters are extracted directly from the weighted topology A, W : ω is is the average weight, represents the average neighbor's residual in-degree, 13,24 and ρ = P (A ij = 1|A ij = 1) describes the network's reciprocity, i.e. the probability to observe a link i → j, given the existence of j → i (for undirected networks ρ = 1).
Equations (5) and (6) represent our key result. They approximate the recoverability of (1), a multidimensional nonlinear dynamic equation, through a manageable first order recurrence relation. This recurrence takes the system's weighted topology (κ, ρ, ω) and its nonlinear dynamics (M 0 , M 1 , M 2 ) as input, and, together with our intervention constraints (∆), predicts the system's recoverability as output. Indeed, for any given forcing ∆, the recurrence (5) either converges to x s (l → ∞) ∈ B 1 , satisfying the recoverability condition (4), or to x s (l → ∞) ∈ B 0 , indicating a failed recovery.
To obtain x s (l → ∞) we extract the stationary states of (5), requiring which in turn provides x s (l) = x s (l − 1). 25,26 Depending on κ, ω and ρ, we observe two characteristic behaviors (Fig. 2, Supplementary section 2.5): Structurally unrecoverable (Fig.  2c,f). In case F (x) and M 2 (x) have a single intersection x ∈ B 0 , then the series in (5) inevitably converges to that point. This captures structural unrecoverability, in which regardless of ∆, single-node reigniting is prohibited. Structurally recoverable (Fig. 2d,e,g). If, on the other hand, F (x) and M 2 (x) have several intersections, then the convergence of Eq. (5) depends on the boundary condition x s (0) = ∆, whose magnitude is determined by our forcing capacity. For ∆ < ∆ c , our forcing is too small, and the system approaches B 0 , a failed reigniting. For ∆ ≥ ∆ c it will reach B 1 , capturing a successful reigniting. The recoverable-phase (Fig. 2i).
If a system is structurally recoverable and ∆ ≥ ∆ c we say that it is in the recoverable-phase, a state in which one can revive the system by forcing just one node.
Taken together, for a given dynamics M 0 (x), M 1 (x), M 2 (x), our formalism predicts a fourdimensional phase-diagram in the κ, ρ, ω, ∆ phase-space, helping us identify the boundaries of recoverability. Next we investigate this phase space on a range of relevant systems, from cellular dynamics (Example 1), to neuronal and microbial systems. In our first examples below, we consider, for simplicity, undirected networks (ρ = 1), reducing our phase-space to three relevant dimensions, κ, ω and ∆. Our final example, reviving a dysfunctional microbiome, examines recoverability under directed interactions (ρ < 1).

Applications
Example 1. Cellular dynamics ( Fig. 3; Supplementary Section 3.1). As our first application we return to Example 1, regulatory dynamics, where M 0 (x) = −x a , M 1 (x) = 1 and M 2 (x) = x h /(1 + x h ), and therefore R(x) = x a and R −1 (x) = x 1/a . Equation (8) under ρ = 1 becomes whose roots (x) determine the potential fixed-points of the reignited system. Clearly, x = 0 is a solution, capturing the fact that the failed state x 0 = (0, . . . , 0) ∈ B 0 is always stable. Hence, the question is, under what conditions do we observe a second solution x > 0, representing a potential convergence to B 1 . To answer this, in Fig. 3e,f we plot M 2 (x) vs. x (yellow) and observe its intersections with F (x) (purple) as we vary the values of κ, ω. This allows us to observe, graphically, the potential convergence of the system to B 0 or B 1 .
First we consider κ = 3, ω = 0.7 (Fig. 3e). We find that (9) exhibits only one solution, represented by the single intersection at x = 0 (red dot). This guarantees that (5) converges to x s (l → ∞) = 0, independently of ∆. Consequently, the system is structurally unrecoverable. Indeed, Fig. 3g indicates that, despite the forcing ∆ at s, the system fails to reactivate.
Increasing the network density to κ = 10, however, changes the picture, as now (9) features three intersection points (Fig. 3f): an intermediate unstable point (white) and two stable points at x = 0 (red) and at x > 0 (green), representing convergence to B 0 and B 1 , respectively. Hence now the system is structurally recoverable, with critical forcing ∆ c = 1 (vertical dashed line), above which it enters the recoverable-phase. This prediction is corroborated in Fig. 3h,i: under ∆ = 0.9 the system remains in B 0 , but for ∆ = 1.1, just above ∆ c , it successfully reignites, precisely as predicted.
This uncovers the existence of a previously overlooked dynamic phase of the Michaelis-Menten model. Indeed, the regulatory system of Fig. 3a has been previously 13 shown to follow two phases, inactive, where only x 0 is stable vs. bi-stable where both x 0 and x 1 are stable ( Fig.  3b-d). Yet we now unveil a third phase: recoverable, a subspace of the bi-stable phase, in which the system can be reignited from x 0 to x 1 by controlling a single node.
To examine this phase-space systematically, in Fig. 3j we present the recoverability phasediagram. For small κ, ω we observe the structurally unrecoverable regime, in which recoverability is unattainable even with arbitrarily large ∆ (red patch in the κ, ω plane). The remaining area in κ, ω represents the structurally recoverable regime, which is split between the unrecoverable phase, for ∆ < ∆ c (below surface), and the recoverable phase when ∆ ≥ ∆ c (above surface).
Setting ∆ constant, in Fig. 3k we construct the ω, κ phase-diagram, by numerically simulating regulatory dynamics on an ensemble of 5 × 10 4 scale-free networks, covering 2, 500 distinct combinations of ω and κ (Supplementary Section 5.3). Each data-point captures the fraction of successful recoveries η among 20 independent reigniting trails, from zero successes (η = 0, yellow) to 100% successful recovery attempts (η = 1, blue). The white solid line represents our theoretical prediction, based on analyzing the intersections of (9). We find that the boundaries of recoverability (yellow/blue) can be well-approximated by our analytical framework. We also present the ω, ∆ and κ, ∆ phase diagrams, further confirming our predicted transitions (Fig.  3l,m).
To test our predictions in an empirical setting, we collected data on two real biological networks, capturing protein interactions in Human 27 (κ = 29) and Yeast 28 (κ = 12) cells. Varying ω and ∆ we measured the reigniting success rate η, observing in both networks the sharp transition into the recoverable phase ( Fig. 3n-q), which falls precisely on the theoretically predicted phase boundary (ω c , vertical dashed lines).
Restructuring guidelines (Fig. 4). In case our system is not in the recoverable phase, we must design appropriate restructuring interventions to push ω or κ towards recoverability. In a cellular environment this can be achieved through biochemical interventions that help catalyze or inhibit specific reactions. 29 The phase-diagrams of Fig. 3j-m, mapping the boundaries of the recoverable-phase, can help us design such interventions.
To illustrate this, in Fig. 4a-d we simulate a cellular network (Yeast) that has been driven towards inactivity due to a series of node/link deletion (grey nodes/links). Some of the removed components are inaccessible (red), and hence cannot be retrieve during restructuring. With these constraints in mind, we incorporate our two-step strategy: Step I. Restructuring. We design a sequence of accessible interventions on A, W to bring the network closer to the recoverable phase. An example of such a sequence is shown in Fig. 4e (along the x-axis). On the y-axis we present the accumulating impact of these restructuring interventions on κ (blue) and ω (orange). To revive the system, we seek paths of such accessible interventions that help deliver the network into the bounds of the recoverable-phase (Fig. 4f). Our goal, we emphasize, is not to simply retrieve the lost components, but to achieve recoverability. This designates, not a single point, but rather an entire sub-space in κ, ω (Fig. 4f, blue area), affording us some level of restructuring flexibility. Indeed, despite the network's irretrievable components, we were able to design three distinct restructuring paths, leading to different destinations -Net 1,2 or 3 -all within the recoverable sub-space (blue area).
Step II. Reigniting. Once in the recoverable phase we can revive the system via single-node reigniting as shown in Fig. 4g for Net 1,2 and 3.
This example illustrates how the phase-diagrams of recoverability provide guidelines for restructuring. For example, in Fig. 4f path 1 builds mainly on controlling the interaction strengths (ω), but assumes little freedom to add nodes or links (κ). In contrast path 3 involves a significant component of adding nodes/links to A, affecting not just ω but also κ. The optimal restructuring path is, therefore, determined by the nature of our constraints, e.g., the relative difficulty in adding weights vs. adding nodes/links. While the potential space of structural interventions in Step I is incomprehensibly vast, our phase-diagrams reduce this space into just two relevant control parameters -κ, characterizing A's density, and ω, capturing its link weights (and ρ in case A is directed). This reduction allows us to asses the contribution of all potential interventions by quantifying their effect on these two (or three) parameters -providing optimal pathways towards recoverability (see further discussion in Supplementary Section 3.1.3).
Example 2. Neuronal dynamics ( Fig. 5; Supplementary Section 3.2). As our second example we consider the Wilson-Cowan neuronal dynamics, 30,31 in which (1) follows the form shown in Fig. 5a. 32,33 The system naturally exhibits two dynamic states ( Fig.5b): suppressed (x 0 , red) in which all activities are constricted; active (x 1 , green) where x i are relatively large (green).
In between these two extremes lies a bi-stable-phase, in which both x 0 and x 1 are potentially stable (grey shaded). This predicts a hysteresis phenomenon, in which a system driven to the left of the grey area will avoid spontaneous recovery.
To observe the predicted phases, in Fig. 5c we numerically analyze our ensemble of 5 × 10 4 scale-free networks. We find, indeed, the active (green) and suppressed (red) phases, separated by the strip of bi-stability (grey). Our formalism, however, predicts an additional dynamic phase -recoverable. This phase, shown in Fig. 5d,e (blue area) identifies a sub-space within the bi-stable regime, under which the system can be driven to x 1 via single-node reigniting. The theoretically predicted phase boundaries are also shown (white solid lines), precisely capturing the numerically observed transitions (yellow/blue).
In the Methods Section we further demonstrate how the brain's modular structure impacts the recoverability phase-space. We also characterize conditions under which modularity provides a fail-safe mechanism, in which one module revives the other upon failure.
Example 3. Microbial dynamics (Fig. 6, Supplementary Section 3.3). As our final example we consider the gut-microbiome, a microbial community, whose functional state has been shown to crucially impact human health. 34,35 Following perturbations, such as antibiotic treatment, the abundance of the different species may reach critical levels, potentially leading to a dysfunctional microbiome. 36 We, therefore, examine below how our two-step revival strategy can help steer a failed microbiome back to functionality.
To construct the interaction network, we collected data 37 describing the metabolic in/out flux of N = 838 microbial species, allowing us to map their complementary chemical interactions. A cooperative link i → j appears when species i produces a resource consumed by j; 38 an adversarial link i ↔ j is assigned if i and j both compete over the same resource. 39 The weight W ij of each link captures the level of inter-species reliance: for example, if i is the sole producer of j's only consumed nutrient m, then j's growth strongly depends on i. If, however, there are many alternative producers of m, or if m is just one of j's many consumed nutrients, then W ij will be small. We arrive at a directed and highly diverse network with broadly distributed weights, and most crucially, with ∼ 75% adversarial interactions, i.e. W ij < 0 (Fig. 6b,c, Supplementary Section 3.3.3). This challenges our theoretical framework, which is primarily designed around cooperative interactions, and therefore helps us test its applicability limits.
To track the microbial populations we used the dynamics shown in Fig. 6a. The self-dynamics combines logistic growth with the Allee effect, 40 and the inter-species cooperative/adversarial interactions follow the Lotka-Volterra dynamics. 22 The parameter F captures the externally introduced microbial influx. Of the original pool of N species, we find that 32% cannot be supported by the network and undergo extinction. The remaining 568 species, comprising the actual microbiome composition, exhibit two potential fixed-points ( Fig. 6d-f), a functional x 1 and a dysfunctional x 0 . Our goal is to apply our two-step recoverability strategy to drive a dysfunctional microbiome at x 0 back towards x 1 .
Selective reigniting. In the Methods Section we show how to evaluate each node's reigniting capacity R s , allowing us to rank all candidate reigniter species in the microbiome, based on their potential to revive the system (Fig. 6g,h). For a random topology, where all shells K s (l) are statistically similar, we expect minor differences in R s between the nodes. Here, however, we find that R s is highly diverse, with the top 26 species having R s > 300, and the remaining species with R s one or two orders of magnitude below (Fig. 6i). Such diversity, a consequence of the unique non-random structure of the microbiome, indicates that, in this system, the top 26 nodes represent preferred candidates for reigniting. To examine this, we attempted reigniting the microbiome with all nodes sequentially, finding that, indeed, success was by far more likely among the top ranked reigniters ( Fig. 6j-l, see also Methods).
One exception we identify among these top reigniters is P. putida, which despite having R s = 359, ranking 6 in the reigniters list, still falls short of reviving the system (Fig. 6k). This is due to the fact that P. putida's reigniting capacity is hindered by its many surrounding adversarial interactions. Indeed, reigniting, by its nature, builds on the positive activation that the source species s exerts on its neighbors at K s (l). Such positive impact is undermined by negative links. Next, we use P. putida's restricted reigniting capacity as an opportunity to examine our two-step strategy in a realistic setting.
Two-step recovery of a dysfunctional microbiome ( Fig. 6m-q). Consider a microbial network at state x 0 , which we wish to recover. Among the many practical constraints on our potential interventions, one crucial constraint is that we lack control over the majority of microbial species, and hence we must achieve recovery with a handful of accessible reigniting nodes. Specifically, let us assume that our top accessible candidate for reigniting is precisely P. putida, which, as Fig. 6k indicates, by itself, cannot revive the system. Hence we employ our two-step strategy.
Step I. Restructuring. To enhance P. putida's reigniting capacity we wish to inhibit the adversarial interactions, which stand in the way of the system's reactivation. We consider two potential interventions: (i) Removing selected nodes that have many adversarial interactions (competition hubs) by means of targeted narrow spectrum antibiotic treatment. (ii) Deleting or weakening adversarial links through nutritional or biological interventions. 41,42 For example, if i and j compete over metabolite m, we prescribe dietary supplements to ensure the availability of m, thus eliminating the i, j competitive interaction (Fig. 6m). Since antibiotic intervention on an already dysfunctional microbiome may cause additional risks, here, we implement restructuring via option (ii). First, we rank all nutrients based on their relative contribution to the adversarial weights in W (Fig. 6n, Supplementary Section 3.3.3). We then supplement the top three nutrients in the list (green) to restructure the network. Eliminating the competition over these, now freely available nutrients, we arrive at our restructured network, which, thanks to our interventions, has now a more suitable balance of cooperative vs. adversarial interactions (Fig.  6o).
Step II. Reigniting. In the microbiome, forcing can be implemented by administering probiotics, a common therapeutic practice that helps sustain, artificially, a desired abundance of a selected species. 43 The rate of the probiotic intake determines the average activation force, set to be above the reigniting threshold ∆ c . Having restructured the network, we now attempt, once again, to reignite it with P. putida. We find that the originally failed reigniter, is now, capable of reviving the inactive microbiome (Fig. 6p,q).
Applicability limits. Our theoretical analysis helps construct the precise phase-diagrams in the κ, ω, ∆ space, predicting the bounds of the recoverable-phase. These analytically tractable observations rely on a specific set of assumptions, mainly, that A features little degree-correlations and has a scarcity of short-range loops, and that the dynamics is of the form (1) and has primarily cooperative interactions. A precise description of these modeling assumptions is provided in Supplementary Section 1.
Our applications, however, provide insights that extend well beyond these analytical limits. For example, our empirical networks, cellular, neuronal and microbial, all exhibit non-negligible deviations from the above assumptions, and yet, our analysis correctly predicted their recoverability (Supplementary Section 5.5). The gut-microbiome analysis helped us derive implementable guidelines for reigniting, including also the notion of selective reigniting, despite the network's strong degree-correlations and significant share of adversarial links (Fig. 6). In Supplementary Section 3.4 we further consider diffusive interactions, in which the governing equation is generalized beyond the form of Eq. (1). We also analyze alternative forms of reigniting, using periodic activity boosts, which in certain applications may be more accessible than the time-independent forcing (∆) considered here (Supplementary Section 4).

Discussion and outlook
While the structure of complex networks has been deeply investigated over the years, our understanding of their dynamics is still emerging. The challenge is often focused on prediction, aiming to foresee a network's dynamic behavior. Here, we go a step further, and focus on influence, showing how to steer a system towards a desired behavior.
Our solution seals a crucial gap in our pursue of nonlinear system controllability. 44 Existing approaches often rely on specific system symmetries, 45-47 which do not cover complex systems of the form (1). Absent such symmetries, complex system control is frequently studied by means of linearization, 48,49 using the linear approximation to help capture the system's local behavior in the vicinity of each of its fixed-point. Such local analysis, however, is insufficient in the context of recoverability, as here we seek to drive the system outside of its current basin and hence beyond the linear regime. To overcome this restriction of locality, cross-basin control was recently developed using time-varying inputs, that adapt until the system is pushed -step by step -into the desired basin of attraction. 14 While highly effective, such interventions require a detailed control over the nodes' dynamics, first, monitoring the system's response, then updating, in real-time, the form of our intervention. Such level of observation/control is not always guaranteed.
To break this gridlock, we seek non-local control, i.e. across basins, but with simple dynamic interventions, that do not require highly detailed input signals. Our recoverablity phase-diagram addresses this by identifying unique conditions where such control is attainable. On the one hand, pushing the system across basins, but on the other hand, using an extremely crude and simple control input -a time-limited constant activation ∆ that is applied to just one or few nodes. No fine-tuning or real-time updating of the input signal is required. Indeed, all that is needed is a strong enough jolt to the system, after which it naturally relaxes to its desired target state.
The microscopic behavior of complex networks is driven by countless parameters, from the finestructure of A and W to the specific rates of each node's dynamic processes. Our analysis, however, shows, that their large-scale functionality can be traced to a manageable set of relevant parameters, i.e. κ, ω, ρ and ∆. Such dimension reduction is the fundamental premise of statistical-physics, allowing to analyze systems with endless degrees of freedom using a limited set of statistical entries. We believe, that such an approach to network dynamics, can help us understand, predict, and ultimately influence the behavior of these complex multi-dimensional systems.
Author contributions. All authors designed the research. HS and BB conducted the mathematical analysis. HS performed the numerical simulations and analyzed the data. AB supervised the microbiome analysis. BB was the lead writer of the paper.
Competing interests. The authors declare no competing interests.
Data availability. Empirical data required for constructing the real-world networks (Microbiome, Brain, Yeast PPI, Human PPI) will be uploaded to a freely accessible repository upon publication.
Code availability. All codes to reproduce, examine and improve our proposed analysis will be made freely available online upon publication.  Depending on the dynamics -e.g., cellular, neuronal or microbialthe system exhibits distinct fixed-points, active (x1, green) or failed (x0, red). Transitions between these states are driven by perturbations to A, W . (c) Unperturbed, the system resides in x1, where all xi > 0. In this presentation, here and throughout, the network nodes are laid out on the x, y plane, and their activities xi are captured by the z coordinate. Hence, an active system has all nodes spread along the positive z-axis, while a failed network is laid-out around z → 0. We also use color coding from red (small xi) to blue (large xi) as visual aid. (d) Perturbations to A, W , such as node/link removal or weight reduction, result in a collapse to the inactive x0. For this system, under x0, all activities vanish (z = 0). (e) Step I. Restructuring. To revive the failed system we first restructure A, W to a point where it can recover, namely a point where x1 is potentially stable. (f) Step II. Reigniting. After restructuring we revive the active state x1 by controlling a microscopic set of nodes, here the single node s. Externally forcing s to sustain constant activity xs(t) = ∆, we drive the network towards x1.
(g) Following reigniting (xs(t) = ∆, black node at center), the forcing signal gradually spreads, until the system's activity x1 is restored. (h) In this process we use the natural basin structure of our dynamics. To reignite x1 we steer the system from B0 (red) to any point within B1 (green). Once in B1, we cease our forcing, and the system naturally transitions to the desired x1 (dashed arrows).   (b) The average activityx vs. the weight reduction q (q = 0 -no reduction). x1 is only stable for q ≤ qc (green); x0 is always stable (red). (c) We track xi for three specific values of q. For small q (top) all xi > 0, an active state x1. As q increases xi become smaller (center), until the system collapses into the inactive x0 (bottom). (d) Retrieving the lost weights does not revive the system since x0 is always stable. To revive the failed system we apply single-node reigniting: (e) F (x) (purple) and M2(x) (yellow) as obtained from the left/right hand sides of Eq. (9) with κ = 3, ω = 0.7, ρ = 1. We observe Case 1 of Fig. 2c, with a single intersection at B0 (red). This describes a structurally unrecoverable system. (f) Under a denser A (κ = 10), F (x) changes form and the system now follows Case 3 -structurally recoverable with ∆c = 1 (vertical dashed-line). (g) Forcing s (black node) to a permanent ∆ activity, no matter how large, fails to reignite the structurally unrecoverable system of panel e. (h) -(i) Reigniting fails for ∆ = 0.9, but succeeds for ∆ = 1.1, confirming our predicted ∆c = 1 in panel f. (j) The κ, ω, ∆ phase-space. The red area in κ, ω captures the structurally unrecoverable regime. We show three specific cross-sections, fixing ∆ (blue), κ (green) and ω (red), to be expanded in panels k-m. (k) The κ, ω phase-diagram, as obtained from 5 × 10 4 scale-free networks (Supplementary Section 5.3). Each data-point represents the fraction η of successful reigniting instances across 20 independent attempts. We observe the predicted phases: inactive (red), where only x0 is stable, unrecoverable (η = 0, yellow), and recoverable (η = 1, blue), where single-node reigniting is possible. Our theoretical prediction is also shown (white solid line). (l) -(m) The ω, ∆ and κ, ∆ diagrams. For ω < 0.6 (l) and for κ < 8 (m) the system is structurally unrecoverable, and hence ∆c → ∞. The theoretical prediction is also shown (white solid lines).  Step I. Restructuring Step II. Reigniting This captures restructuring constraints that are, indeed, inevitable in realistic scenarios. Circle at center -we focus on the same sub-network shown in panel a, the deleted nodes/links appear in grey and red, the remaining unperturbed components are highlighted, . (e) Step I. Within the given constraints we restructure the network by reintroducing nodes/links or strengthening link weights. We map these interventions into their impact on the two relevant control parameters κ (blue) and ω (orange). For illustration, we show the highlighted sub-network of panel c as it restructures, acquiring nodes, links and increased weights (sub-networks along the x-axis). (f) Restructuring paths in the κ, ω phase-diagram. A successful path leads the network from the collapsed phase (red) into the recoverable-phase (blue). Using our predicted phase-diagram we design several alternative paths, affording us flexibility to, e.g., focus on increasing weights ω (Net 1, light-blue path) or also on enhancing network density κ (Net 3, purple path), all depending on the nature of our constraints. (g) Step II. Once the network is in the recoverable-phase, we can revive it via single-node reigniting, demonstrated here on each of our restructured networks, Net 1-3. Each data-point represents 20 independent realizations. As predicted, the bi-stable state is split into two distinct dynamic phases -unrecoverable (yellow) vs. recoverable (blue). Simulation results (yellow-blue transition) are in perfect agreement with our theoretical predictions (white solid lines). (f) The empirically constructed brain network, and its two hemispheres M1 (blue) and M2 (red). (g) We measured the reigniting success rate η of the modular brain network along the ωIntra, ωInter phase-diagram (see Methods). We observe three phases: unrecoverable, in which no module can be revived (η = 0, yellow); recoverable, in which both modules reactivate (η = 2, blue); modular, in which reigniting is constrained to a single module (η = 1, green). (h) Setting ωInter = 1, ωIntra = 2, in the unrecoverable phase (yellow circle), we find, indeed, that reigniting fails to revive both modules. (i) In the modular phase (green circle), M1 recovers, but fails to reactivate M2. (j) In the recoverable phase (blue circle) reigniting successfully crosses over from M1 to M2. (k) The average activityx of M1 (blue) and M2 (red) vs. t, as obtained from numerically simulating neuronal dynamics on the brain network. Adding noise, we observe sporadic fluctuations, some causing a transition to the suppressed state. The first collapse occurs to M1 at t ∼ 20. Soon after, M2 also collapses at t ∼ 40, and the entire network irreversibly fails. (l) Increasing the inter-module weights to ωInter = 5, we enter the recoverable phase. Now when one module fails, the other module reignites it. This captures a fail-safe system, whose modular structure provides internally embedded self-recoverability.

Rank
Nutrient Impact  The export (orange) import (green) network, linking microbial species to the nutrients they produce/consume. Cooperative links (blue) are assigned when i's exports are imported by j; adversarial links (red) -when i and j share mutual import/s. (c) The resulting gut-microbiome, a diverse weighted and directed network with mixed cooperative (blue) and adversarial (red) links. (d) -(f) The system exhibits two states, the functional x1 (green diamond) or the dysfunctional x0 (red diamond), with a broad range of bi-stability (grey shaded) (g) The propagation of a reigniting signal: s → i → j. In Methods we show that an effective i → j link must satisfy Wi→j > w0. (h) Counting only effective links, we obtain s's reigniting capacity Rs, as the number potentially reactivated nodes (blue). (i) We ranked all microbial species based on Rs, identifying the network's most effective reigniters. The top 26 nodes stand out with Rs orders of magnitude higher than the rest. (j) The visibly dense effective network surrounding B. bifidum (Rs = 360). As expected, B. bifidum successfully reignites the microbiome (right). (k) -(l) We examined two lower ranked species, which failed to revive the network. Specifically, P. putida, despite having high Rs, causes only a local impact but fails to revive the network. (m) To reignite with P. putida we employ our two-step strategy. In Step I, we use nutritional supplements (green) to eliminate competitive links (red dashed-link) by ensuring nutrient availability (n) We ranked nutrients based on their contribution to the adversarial links. The top three ranking nutrients are highlighted (green). (o) Supplementing these three nutrients via dietary interventions, we restructure the microbial network, reducing adversarial interactions (red). See insets for visible comparison. (p) In step 2 we administer probiotics to boost the population of a specific species. (q) After restructuring we find that P. putida, an originally failed reigniter, can now successfully revive the microbiome.

Methods
Recoverability of modular networks. Applying our formalism to an empirically constructed brain network, 1 allows us to examine its predictive power beyond our analytical assumptions of a random A. Indeed, the brain, with its two hemispheres, provides a highly structured (nonrandom) modular network, partitioned into two clearly distinctive communities M 1 and M 2 (Fig. 5f). It, therefore, offers meaningful insights on recoverability across modules. The question is, can reigniting one module, say M 1 , spillover to also revive M 2 . Our analysis indicates that this depends on the average strength of the links within the modules, ω Intra , vs. that of the links between the two modules, ω Inter (Fig. 5g). Clearly, if ω Intra is too small, then both M 1 and M 2 are, in and of themselves, unrecoverable, and reigniting will inevitably fail (Fig. 5h). If both ω Intra and ω Inter are sufficiently large, then the reactivated nodes at M 1 will further reignite their neighbors at M 2 , allowing a complete recovery of both modules via single-node reigniting ( Fig.  5j). In between these two extremes, we observe a third phase, in which reigniting is confined to M 1 , but fails to penetrate M 2 (Fig. 5i).
The result is a three state phase-space, as shown in Fig. 5g: recoverable (blue), in which reigniting at M 1 can reactivate also M 2 , unrecoverable (yellow), in which both modules cannot be revived, and modular (green), where M 1 recovers, but the reigniting signal fails to penetrate M 2 . To construct the phase-diagram of Fig. 5g we simulated neuronal dynamics on the brain network of Fig. 5f, scanning 2, 500 distinct combinations of ω Intra and ω Inter . For each combination we attempted 20 independent realizations of reigniting with randomly selected nodes. We then extracted η, as the average number of revived modules over the 20 attempts. Hence, η = 0 means that no module was revived (yellow), and η = 1 indicates that (on average) reigniting was restricted to only the source node's module (green), but did not cross over to the other module. Finally, if η = 2 then both modules were reactivated (blue).
This observation highlights the potential benefits of network modularity for self-recovery. Indeed, if one module fails, say M 2 , the other module, M 1 , if still active, can revive it. This is because M 1 's active nodes can themselves help reignite the inactive M 2 . Hence, modularity offers a fail-safe network architecture, in which, unless both modules simultaneously fail, one module can reactivate the other. This ensures a sustained activity in the face of sporadic failures. To observe this we simulated neuronal dynamics on our brain network, setting the inter-modular link weights to ω Inter = 2. Starting at x 1 , we introduce external noise, that causes the activity of both modules to fluctuate, until a sufficiently large perturbation, occurring at random around t ∼ 20, leads M 1 to irreversibly fail (Fig. 5k, blue). Shortly after, a similar fate meets M 2 (red) and the entire system collapses to x 0 .
We now repeat the exact same experiment, this time with ω Inter = 5, bringing the system into the fully recoverable-phase (η = 2, blue). Now, at every instance in which, say M 1 fails, M 2 's active nodes reignite it back into activity (Fig. 5l). Hence, modularity can afford the system a fail-safe dynamics, driven by its capacity for self-recoverability.
Selective reigniting in the microbiome. A unique aspect of the our empirically constructed microbiome network, is that it deviates significantly from a random topology. Indeed, in a random network, as degree-correlations are negligible, the statistical properties of the shells K s (l) become approximately independent of s for l > 1. In simple terms, while nodes may have diverse degrees, i.e. K s (1) = K s (1), their second or third neighbors at K s (l > 1) are statistically similar. Under these conditions, once a system is in the recoverable phase, one can reignite it using any desired node, as, indeed, all nodes have roughly identical shells in their surrounding.
If, however, the environments K s (l) vary significantly across the nodes, we expect that certain nodes become better reigniters than others. This gives rise to selective reigniting, in which the system's recoverability is not just a function of the network, but also depends on the specific source node s and its unique reigniting capacity.
Consider the reigniting signal as it propagates from the source s, and penetrates through the shells K s (1), K s (2), . . . . At a certain instance it revives a node i ∈ K s (l), then advances from i to impact its neighbor j ∈ K s (l + 1) and so on (Fig. 6g). In Supplementary Section 3.3.3 we show that such propagation across shells, namely that a revived i can, indeed, reactivate its more distant neighbor j, requires that the i → j link weight exceeds a threshold w 0 . Hence, the links whose weight W ≥ w 0 , constitute effectual links, that help propagate the reigniting signal. The remaining links with W < w 0 have but a marginal contribution to the reigniting. We, therefore, construct an effective network, comprising only the effectual links (Fig. 6h, blue links/nodes). This, more selective, network represents the relevant set of interactions that actively participate in the reigniting process. Constructing this network we obtain the effectual shells surrounding each node s, whose number of nodes captures s's reigniting capacity R s . Nodes with large R s are surrounded by many effectual links, and therefore they have a higher reigniting capacity.
In a random topology, where all shells are statistically similar, we expect minor differences in R s between the nodes. In the microbiome, however, we find that R s is highly diverse, with the top 26 species having R s > 300, and the remaining species with R s that is one or two orders of magnitude below (Fig. 6i). Such diversity, a consequence of the unique structure of the microbiome, indicates that, in this system, the top 26 nodes represent preferred candidates for reigniting. To examine this in Supplementary Fig. 7 we show η vs. the R s -ranking for all nodes in the microbiome. We clearly observe that the top ranked nodes have a much higher probability to successfully reignite the system.

Modeling framework
We consider a class of systems captured by is node i's dynamic activity (i = 1, . . . , N ) and the nonlinear functions M 0 (x), M 1 (x), M 2 (x) describe the system's intrinsic dynamics, i.e. its self-dynamics (M 0 ) and its interaction mechanisms (M 1 , M 2 ). The matrices A and W provide the network's topology and weights, respectively, tracking which components interact (A ij = 1), and how strongly (W ij ). In (1.1) we focus on cooperative interactions, in which nodes positively impact each others activity, hence, we take the interaction term to be always positive, namely W ij and the product M 1 (x)M 2 (x) are all ≥ 0. We also assume M 2 (x) is monotonous, excluding, e.g., oscillatory coupling functions. This ensures that the external activation of node s positively impacts all other activities, and thus helps reactivate all other nodes, indeed, a crucial component of the reigniting strategy.
Weighted topology. We take the network topology A, an N × N matrix with an average degree k , to be large (N → ∞), sparse ( k N ) and potentially directed. The network, we assume, is extracted from the configuration model ensemble, in which the nodes are characterized by an arbitrary degree sequence, but are otherwise randomly linked. Hence, A can follow any in/out degree distribution, P (k in , k out ), but lacks higher order structure, such as degree-correlations, clustering or modularity.
The weight matrix W captures the fully positive rates of all existing interactions, as extracted from the probability density P (w). From this density we obtain the probability P (w) dw for a randomly selected weight to follow W ij ∈ (w, w + dw). The weights are assigned at random from P (w) to all links. We denote the expectation of P (w) by where 1 = (1, . . . , 1) and ⊗ represents an element by element product (Hadamard). In the r.h.s. of (1.2) we approximate the expectation E(w) via the average weight, summing over all existing link weights A ij W ij (numerator) and dividing by the total number of links A ij (denominator).
Directionality. In (1.1) a link A ij captures a connection in which j affects i's activity x i . In that sense A ij = A i←j an incoming link from j to i. Hence i's in-degree is defined as k i,in = N j=1 A ij , counting all nodes that affect i, and its out-degree as A ij , namely all nodes affected by i. Therefore, to track the pathways of signal propagation from one node to the other we use A , as, indeed, A ij captures a link from i to j, i.e. an interaction where i impacts j. We can now use this transposed topology A to track the network paths, capturing a sequence of (outgoing) links A sj , A jn , . . . , A mi , leading from the source s to the target node i. Such paths track to potential propagation of a signal from s to i. The minimal path length between s and i defines the distance L si = L s→i between these two nodes.
In a fully random directed (sparse) topology the majority of links are strictly one-directional. This is because the probability for creating a short loop of the form A ij = A ij = 1 is tiny, scaling as ∼ 1/N . Still, in real-world networks such reciprocity is rather common. Indeed, often if i affects j, j also affects i. We quantify this by the network's reciprocity coefficient the conditional probability for a j → i link given the existence of one from i → j. In a fully random topology, with no inclination for reciprocal links we have ρ → 0, having almost no bi-directional links. Conversely in an undeirected network, where all links are reciprocal, we have ρ = 1. For example, in genetic interactions, most genes are either transcription factors (regulators) or regulated, hence genetic interactions are typically non-reciprocal, most links are one-directional and ρ is small. In contrast, protein interactions, capturing chemical relationships between proteins, are fully reciprocal, A is undirected and ρ = 1. Finally, in population dynamics, competition or symbiosis are, by their nature, bi-directional, while complementarity, an interaction in which species i produces a product needed for species j, is directional. Hence in population dynamics we have 0 < ρ < 1.
Shells. Using L si we define the shells surrounding s via comprising all nodes at distance l from s. Hence, we have K s (0) = {s}, K s (1) = s's group of direct outgoing neighbors, K s (2) its second neighbors, etc. Nodes in K s (l) can be potentially affected by s through directed pathways of length l. In a similar fashion, we denote by the incoming l shell around s. This group consists of all nodes that can affect s from a network distance l. Specifically Q s (1) is the groups of s's incoming neighbors, allowing us, e.g., to write Eq. (1.1) as in which we replace the term N j=1 A ij · · · by j∈Q i (1) · · · . We denote the size of group X by |X|, writing, e.g., k s,out = |K s (1)|.

Resulting approximations
Under these assumptions the weighted topology A ⊗ W features three characteristics, which we use below to derive our recoverability formalism: 1. Locally tree-like. The random patterns of connectivity ensure the scarcity of short range loops (barring reciprocal links of the form A ij = A ij = 1, which are technically a loop of length two). Hence, the shells K s (l) have (almost) no intra-shell links.
Consequently, a typical node i ∈ K s (l) links to a single neighbor from K s (l − 1) and its remaining k i,in − 1 neighbors are from K s (l + n), n = 1, 2, . . . .

Shell statistics.
While the (weighted) degrees of all nodes may vary, at times -across orders of magnitude, their surrounding shells K s (l), comprising a random selection of neighboring nodes are statistically similar. Hence, for example, the average degree of nodes in K s (l), l ≥ 1, is approximately independent of s or of l. In simple terms: s may have many more neighbors than s , i.e. |K s (1)| |K s (1)|, however s's |K s (1)| neighbors and s 's |K s (1)| neighbors are extracted from the same statistical pool, and therefore have similar statistical properties.
3. Cooperativity. With the interactions being of cooperative nature, the effect of reigniting is consistently positive, i.e. s activates its neighbors at K s (1), which in turn, may also activate their neighbors at K s (2) etc., capturing a sequence of recovering interactions.
In Sec. 2 we use all of the above features to derive the recurrence relation of Eq. (5) of the main text. Hence, our analytical framework is exact in the limit where the listed assumptions 1 -3 hold, and becomes approximate in case the network A or the weights W do not adhere to the configuration model framework (characteristics 1,2), or if some interactions are of adversarial nature (characteristic 3). Still, our analysis shows that our derivation is quite robust against such discrepancies, and continues to offer both quantitative and qualitative insight. Specifically, we examine: • Empirical networks, such as Yeast, Human, Brain and Microbiome, which are known to be imperfect. All of these networks exhibit a non-negligible density of short loops (e.g., clustering) and measurable degree-correlations. The Brain network also has a distinctive modular topology. Still, our results indicate that recoverability is equally relevant in these real-world network topologies (see Sec. 5.5).
• The Microbiome incorporates a mixture of cooperative and adversarial interactions, and yet, as we discuss in Sec. 3.3.3, as long as the positive links dominate, recoverability is attainable, and continues to follow our theoretical guidelines. In Sec. 3.4 we study the recoverability of ferromagnetic interactions, which violate this factorizability condition of Eq. (1.1). While we cannot apply our analytical framework on such dynamics, we can examine their behavior numerically. We find, numerically, that also in these broader forms of dynamics, the system exhibits a recoverable phase. This indicates, once again, that the insights offered by out formalism range beyond the restrictions of our modeling framework.
States of the system. The fixed-points x α = (x α,1 , . . . , x α,N ) of the system can be obtained by setting the derivative on the l.h.s. of Eq. (1.1) to zero. If dynamically stable, these fixed-points represent potential dynamic states in which the system can reside. The specific activity of a node x α,i in each of these states may generally depend on that node's weighted in/out degree s i,out = j∈K i (1) W ij , s i,in = j∈Q i (1) W ij . Therefore, since A, W may be extremely heterogeneous, the individual node activities can potentially be highly diverse. Fortunately, as explained above, in the configuration model, while individual nodes vary, the shell statistics, are approximately uniform, namely the nodes in K s (l) are drawn from a similar statistical pool as those belonging to K s (l ). We, therefore, in our derivation below, characterize the activity of the system in each of its states x α using the shell-averagē (1.7) This represents an average over the activities of all nearest neighbor nodes, first summing over the direct neighbors of s, K s (1), then covering all s from 1 to N . In (1.7) we use |X| to denote the number of nodes in group X. To denote the fixed-point we simply omit the t dependence on both sides on (1.7), writingx α instead ofx α (t). As a matter of convention, in our analysis we assume a two-state system: an undesirable state x 0 and a desirable state x 1 . This can be readily generalized to treat more states, if relevant.
Irreversible collapse. The class of systems within (1.1) that we are particularly interested in are ones that may suffer irreversible collapse. This rather common phenomenon occurs when the system exhibits multi-stability, for example, a bi-stable state where both the undesirable x 0 and the desirable x 1 are potentially stable. Such conditions arise quite frequently, as observed in the set of dynamical systems we consider here: Regulatory, Brain and Microbial. Since x 0 is stable, the system does not recover on its own, and therefore mandates dynamic intervention or reigniting, as we derive in Sec. 2.

Basins of attraction.
Under bi-stability, the state of the system, x 0 or x 1 , depends on the initial condition x i (t = 0). This splits the N -dimensional state space of the system x into basins which comprise all initial conditions from which Eq. (1.1) converges to the fixed-point x α .
Reigniting a system's failed dynamics, therefore, translates to steering it from the state x 0 to the basin B 1 .

Modeling single-node reigniting
Consider a system of the form (1.1), characterized by two stable states -an undesirable x 0 and a desirable x 1 . Let us further assume that the system exhibits a bi-stable phase, in which both states are potentially stable, and that it is presently at the undesirable x 0 . We seek dynamic interventions, preferably minimal in nature, that will help us drive the system towards x 1 . To achieve this we assign a selected set of nodes F -the forced nodes -whose dynamics we externally control. This, effectively, changes the system's dynamics (1.1) into in which nodes in F are forced to follow the external control function φ i (t), while the remaining N − |F| nodes continue to evolve via the system's natural interaction dynamics.
In a realistic reigniting scenario we require |F| N and φ i (t) to be described by preferably simple functions, capturing the fact that often we have limited access or control over the dynamic behavior of the majority of the nodes. Taken to the limit, we choose Such single node reigniting is, in principle, no different than reigniting by few nodes, since, for a large network, the immediate neighborhoods of each two randomly selected nodes have negligible overlap. Therefore, the impact of one node's forced activity has little interference with that of the other. In other words, for a randomly selected microscopic set of nodes, the group F will be, most probably, spread throughout the network, comprising an isolated set of reigniting focal points, each impacting only its local neighborhood (Fig. 1a). Under these conditions, the forcing of more than one node does not significantly contribute to the reigniting, and only begins to take effect if the recovered neighbors from one forced node overlap with those of another. Such overlap occurs only if each of the forced nodes has by itself reignited a large fraction of inactive nodes, which, in principle depicts an independent set of single-node reigniting instances. Therefore reigniting via |F| N is, asymptotically identical to reigniting via |F| = 1 -a single forced activity focal point that must penetrate the network and impact its distant neighbors towards the basin B 1 .
To evaluate the impact of s's forcing we track the response of the shells K s (l), as defined in Eq. (1.4). The average activity in K s (l) is captured by  Figure 1: How a single node impacts the network. (a) Reigniting by two nodes s 1 , s 2 (or any microscopic fraction of nodes) introduces two focal points of reactivation into the network. For a large network, the adjacent shells around s 1 and s 2 have little overlap, hence each only impacts its local neighborhood. Under these conditions s 1 's forcing does not meaningfully reinforce that of s 2 , and each behaves as an independent single-node reigniting source. In case s 1 is able to excite the distant shells, to the extent that its impact interferes with that of s 2 , this indicates that s 1 's activation had a non-local impact. This is the case where s 1 indeed successfully reignites the entire system. Therefore, we find that reigniting by a microscopic random set of nodes is, in effect, no different that reigniting with a single node. (b) The impact of our reigniting at s on a node i ∈ K s (l) is mediated by i's direct incoming neighborhood Q i (1) (red). In a large random network, whose structure is locally tree-like, this neighborhood has typically a single node in K s (l − 1), closer to the source, and κ nodes in K s (l + n), farther away from the source. This allows us to evaluate n in (2.19), and obtain the recurrence relation of Eq. (2.34).
x s (l, t) = 1 |K s (l)| i∈Ks(l) x i (t), (2.4) allowing us to evaluate the impact of the forcing x s (t) = ∆ at any distance l from s.
Being initially at the undesired state we begin with x s (l, t = 0) =x 0 for all l > 0. Indeed, x 0 in (1.7) captures precisely the average shell activity in the α = 0 state. For l = 0 we set x s (l = 0, t) = ∆, as per our forcing intervention. In a non-recoverable system the s-forcing fails to reignite the system, its impact remains local, and hence the distant shells continue to be in the undesired state's basin of attraction, namely Conversely, if the system is recoverable, s's forcing penetrates the network to impact the entire system, thus leading to in which the distant shells have been successfully reignited. Under these conditions we can terminate our forcing, to allow s's local neighborhood, whose state is frozen by our external intervention, to also transition to x 1 via the system's natural, undisturbed, dynamics. Note that, strictly speaking, B α is a group of vectors of the form x, while x s (l → ∞, t → ∞) is a scalar value, capturing the shell average. Therefore, the sign ∈ in (2.5) and (2.6) should be taken to mean that x s (l → ∞, t → ∞) represents the a shell average that is associated with an initial condition x that is within B α .
To obtain a direct set of equations for x s (l, t) in (2.4), we write which using (1.1) to express the r.h.s. derivative provides (2.8) To approximate the summations over M 0 (x i ) and M 1 (x i ) we use a mean-field approach, writing . This approximation is exact if at least one of the following two conditions applies: (i) M 0 (x) and M 1 (x) are linear; (ii) x i (t) are uniform within the shell K s (l). Clearly, these conditions are not guaranteed, however, under many practical scenarios, they represent a sufficient approximation, designed to detect the macroscale behavior of the system -as fully corroborated by our numerical examination. Indeed, while Eq. (1.1) is, generally, nonlinear, its components, M 0 (x), M 1 (x), in many of the useful models, are often sub-linear, linear or weakly super-linear, i.e. involving powers that are not much higher than unity. This satisfies, approximately, condition (i). In other cases we may observe strong nonlinearities in M 0 (x) and M 1 (x), e.g., in our microbial dynamics, but in such cases, we often have bounded activities x i (t). This ensures a narrow distribution of x i (t), roughly satisfying condition (ii). We further elaborate on the relevance of these condition in the appropriate sections, where we analyze each of our specific dynamic systems (Sec. 3).
Using approximation (2.9), we rewrite (2.8) as where we have also replaced the summation N j=1 A ij · · · by j∈Q i (1) · · · , as shown in Eq. (1.6) above. According to our random network construction the weights W ij are randomly assigned from P (w). Therefore, summing over the weights as we do in (2.10) can be replaced by extracting the average weight ω (1.2) from the summation, providing Next, we consider the nodes in Q i (1), appearing in the second summation term of (2.11). This sum includes all nodes j that are nearest incoming-neighbors of i ∈ K s (l), namely neighbors of a node that is located at distance l from s. Therefore, by definition, j is in one of the shells K s (l − 1), K s (l), K s (l + 1), . . . , allowing us to write j ∈ +∞ n=−1 K s (l + n). (2.12) We can use (2.12) to express the sum over j ∈ Q i (1) appearing in Eq. (2.11) as splitting the group Q i (1) into a disjoint union Q i (1) = +∞ n=−1 Q i (1) ∩ K s (l + n). Next we employ the mean-field approximation of (2.9), this time on M 2 (x), to write (2.14) in which we replace x j (t) from the l.h.s. with the shell average x s (l + n, t). Together, we can now use (2.14) to express the sum over j ∈ Q i (1), which appears in Eq. (2.11), as denotes the fraction of nodes within i's neighborhood Q i (1) that belongs to the shells K s (l + n) 1 .
The r.h.s. of Eq. (2.15) represents an averaging of i,n (s, l) and of Q i (1), as carried out over all nodes i ∈ K s (l). However, thanks to characteristic 2 of Sec. 1.1, stating that the shell statistics are independent of s or l (for l ≥ 1), we can replace the K s (l) average by the network average. The meaning is that, for any node characteristic X i , we can write where the l.h.s. represents an average performed selectively over nodes in K s (l), a specific shell, and the r.h.s. uses the independence on s and l to (i) average over all potential source nodes s, and (ii) replace the specific l-shell with l = 1, as indeed, on average, all shells are the same. Approximation (2.17) becomes exact under the conditions discussed in Sec. 1. In the context of i,n (s, l) and Q i (1), we can use (2.17) to further simplify the summation in (2.15) into (2.20) The two parameters n and κ represent global network characteristics, independent of the selected source node s, or of the specific shell K s (l) -solely determined by the network topology A. The first, n , captures the network expansion, seeking, for a typical node at distance l from any arbitrary source node s, how many of its direct neighbors are at greater (n > 0), equal (n = 0) or smaller (n = −1) distance from that source (Fig. 1b). Note that in (2.19) we focus specifically on l = 1, but we could have equally used l = 2, 3, . . . instead, under the assumption of the similar shell statistitics. In our discussion below, we, therefore, treat n as extracted from an arbitrary shell l > 0, not necessarily l = 1. The second parameter, κ is the average residual in-degree of the network, which quantifies, for a typical neighbor node, i.eȯne that was reached via an incoming link, the number of its remaining incoming links.
To evaluate n in (2.19) we use characteristic 1 of Sec. 1.1 -the fact that our networks are locally tree like. In such networks the shells inflate exponentially as |K s (l)| ∼ e ξl , a consequence of the tree-like topology, and hence the randomly wired links have a vanishing probability to link to the inner shells, and preferably connect to nodes from K s (l + 1), K s (l + 2), . . . . Under these conditions, a typical node i ∈ K s (l) links primarily to the most distant shells at l → ∞. Of course, it must have at least one neighbor in K s (l − 1), as indeed, absent such a neighbor i / ∈ K s (l), but then all its remaining |Q i (1)| − 1 neighbors, will be situated at the most distance shells K s (l → ∞).
Our networks, however, are not fully random, due to their reciprocity ρ in (1.3). The meaning of this is that of i's |Q i (1)| nearest incoming neighbors, a fraction ρ are also outgoing, i.e. reciprocal. To understand how this deviation impacts the inter-shell link distribution, let us consider a node j ∈ Q i (1) who is within this ρ fraction. The node j is an incoming neighbor of i, and therefore A ij = 1, but it is also an outgoing neighbor, and hence A ij = 1 as well. Now since i ∈ K s (l), there exists a path of l outgoing links from s to i, in the form A sn , A mn , . . . , A qi . We can now concatenate the link A ij to this sequence to obtain a path of length l + 1 leading from s to j. Therefore we have j ∈ K s (l + 1). This implies that out of i's |Q i (1)| nearest neighbors we have: a single neighbor in K s (l − 1), i.e. |Q i (1) ∩ K s (l − 1)| = 1; (2.21) and the remaining |Q i (1)| − 1 neighbors split into a ρ fraction in K s (l + 1) and a 1 − ρ fraction in K s (l → ∞), namely Note that, on average, none of i's neighbors is within the K s (l) shell, indeed -loops are scarce in a tree-like topology -and therefore |Q i (1) ∩ K s (l)| → 0. This predicts in (2.19) that Using (2.20) we arrive at the network's average shell expansion Equations (2.28) -(2.31) describe the average fraction of links from K s (l) that lead to K s (l − 1), K s (l), K s (l + 1) and K s (l → ∞), respectively. In the limit of ρ = 1, an undirected network, we have all links pointing either to K s (l − 1) or to K s (l + 1) (Fig. 1b). In the opposite limit of ρ = 0, a purely directed network, we have all links, but the one directed to K s (l − 1), connecting to distant nodes at l → ∞. For 0 < ρ < 1, a directed network with finite reciprocity, we have a fraction of local links, connected to nodes from K s (l + 1), and a remaining group of distant neighbors at K s (l → ∞).
We now use (2.20) and (2.28) -(2.31) to substitute (2.15) into (2.11), obtaining dx s (l, t) dt = M 0 x s (l, t) + ωM 1 x s (l, t) M 2 x s (l − 1, t) + ρκM 2 x s (l + 1, t) Beginning at an initial condition where the system is at the undesired x 0 , the shell-average, as defined in Eq. (1.7), is initiallyx 0 at all shells l ≥ 1. Of course, this state changes with t, however, for x s (l → ∞, t), the most distant shells, we set the boundary condition to x s (l → ∞, t) =x 0 . This provides us with tracking the dynamics of the average activity in the l-shell, as driven by the system's internal dynamic mechanisms M 0 (x), M 1 (x) and M 2 (x). It uses the tree-like structure of the network shells to reduce the detailed network topology A ij in (2.1) to a simplified form, focusing on the average node in K s (l), which, typically, interacts with a single node in K s (l − 1) and κ residual nodes in K s (l + 1) or K s (l → ∞). This equation is valid for all shells K s (l) at l > 0; for K s (0) we use the first equation of (2.1), providing x s (0, t) = ∆. Together we arrive at capturing the time evolution of nodes at distance l from the forced s; the initial condition is set to x s (l, t = 0) =x 0 for all l > 0. Below in our simulation of cellular and brain dynamics (Sec. 3) we consider Eq. (2.34) in the case of undirected networks. Under these conditions the equations can be simplified by setting ρ = 1, which provides . (2.35)

Steady-state analysis
To observe the system's response to our single-node reigniting we wish to obtain the fixedpoint of our recurrence relation. Hence, we set the derivative on the l.h.s. of (2.34) to zero, Note that we have now omitted the term t from x s (l, t), focusing on the steady-state i.e. the system's final, long term, activity patterns. To isolate x s (l) we rewrite the second equation in (2.36) as This, by inversion, provides a direct expression of x s (l) in function of its two neighboring terms, x s (l ± 1).
We can now substitute (2.40) into the second equation in (2.36) to transform it into a second order recurrence relation, obtaining The challenge is that the recurrence (2.41) is ill-defined, as we only have one boundary condition, x s (0) = ∆, instead of the two anchoring points required to obtain a unique solution. Hence, in and of itself, Eq. (2.41) cannot predict the final shell states x s (l), and therefore it is insufficient to determine if our reigniting is successful or not. Next, we introduce an approximate approach that allows us to track the desired fixed-points of (2.41).

Predicting the final shell states
While in general (2.41) is under-determined, we can use our prior knowledge on the states of Eq. (1.1) to constrain its potential solutions. Indeed, knowing that our system has potentially two stable fixed-points, x 0 and x 1 , we assume that our forcing at x s (0) can lead, asymptotically, to only two outcomes: successful reigniting, in which 42) or unsuccessful reigniting, where Therefore, we do not need to solve the recurrence relation fully, just to determine whether it assumes the asymptotic solution (2.42) or (2.43).
We begin by expressing x s (l + 1) in (2.41) as obtained by substituting l with l + 1 in the recurrence relation. This allows us to rewrite the recursive series as expressing the term x s (l + 1) in (2.41) via Eq. (2.44). This step provides x s (l) in terms of x s (l − 1), x s (l) and x s (l + 2). Consequently, the average activity at K s (l) is impacted by the state of the directly neighboring shell K s (l − 1), by the equidistant nodes at K s (l) itself, and by the indirectly interacting second neighbors at K s (l + 2). Our main assumption is that of these three effects -the first two, which represent K s (l)'s direct neighborhood, supersede that of the third x s (l + 2) term. Indeed, this term captures the state of the distant shell K s (l + 2) whose impact on x s (l) is marginal as compared to the other two terms. We therefore approximate this term by assuming this distant shell has not been significantly impacted by our reigniting, hence still at its initial undesired state x 0 , having on average activityx 0 . This discrepancy, we expect, will have little impact on the accuracy of x s (l), as it only applies to l's distant neighbors. Using (2.46) to rewrite (2.45) we now have from which we can extract x s (l) as recovering Eqs. (5) and (6) of the main text. Under ρ = 1, an undirected network, F (x) takes the form

Evaluatingx α
The components of (2.49) and (2.50) can all be extracted directly from our given parameters: the functions M 1 (x), M 2 (x), M 3 (x) and hence also R(x) are provided within Eq. (1.1), κ and ω are extracted from A and W via (2.20) and (1.2), and ∆ is determined by our forcing capacity. This leaves the shell-averagex 0 in (2.50), as the only unknown in this recurrence. We, therefore, now discuss how to calculatex 0 andx 1 -the natural states of the system, ans specifically how to extractx 0 in (2.50). We offer three different options by which to evaluate these parameters:

Option 1. Settingx 0 to zero
The most crude approach is to setx 0 = 0. Indeed, in the typical scenario we have the undesired state representing a system under suppressed activity in which all x i are small. Under these conditionsx 0 x 1 , and can be approximated by a zero state. This approach represents a rough approximation for, e.g., Brain or Microbial dynamics, but becomes exact for Regulatory, where the undesired state is, indeed, x 0 = (0, . . . , 0) . Such approximation allows us to solve the recurrence relation directly, predicting whether single-node reigniting will, indeed, revive the system. Yet, it offers no information on the recovered statex 1 .

Option 2. Mean-field approximation
We can approximate analytically the dynamic states of the system using a mean-field approach, which was originally presented in Ref. 2 . This approach is exact under the modeling framework presented in Sec. 1, and specifically builds on characteristic 2 in Sec. 1.1 -the statistical similarity of all shells K s (l). Our presented results indicate that it also holds quite well even when considering realistic networks, which show some discrepancies with regards Sec. 1's assumptions.
Starting from the definition ofx α in (1.7) we can use (1.6) to write a direct equation for the shell-average as (2.52) Let us first treat the summation over Q i (1) on the r.h.s., by writing it as (2.53) Using the mean-field approximation of (2.9) we obtain where ω i = |Q i (1)| −1 j∈Q i (1) W ij is the average weight of i's incoming links, and x i← (t) is the average activity of i's incoming neighborhood Q i (1) (analogous to x s (l, t) in (2.4), defined for outgoing neighborhoods). Substituting (2.54) into (2.52) we arrive at (2.55) Averaging over all i ∈ K s (1), as we do in (2.55), allows to substitute ω i , the average weight surrounding i, by ω(s, 1), the average weight surrounding the K s (1) shell. Further, as we average over all s, we can extract the global weight ω, as defined in (1.2), outside of the sum. |Q i (1)| = κ + 1. (2.56) Finally, absent degree-correlations i's incoming neighbors Q i (1) are, on average, no different than the incoming neighbors of any other node, allowing us to express the average over all x i← (t) as (2.57) the global shell-average of (1.7). This, once again, is a direct consequence of our configuration model network ensemble, and specifically of characteristic 2 in Sec. 1.1. We arrive at (2.58) on which we apply, once again, the mean-field approximation of (2.9) to obtain To construct a specific equation forx α (t) we use the averaging over all s on the r.h.s. of (2.59) to substitute the s-shell average, x s (1, t), by the network shell-averagex α (t). This provides us with is a parameter fully determined by the weighted network topology A ⊗ W .
Equation (2.60) provides us with a single equation for the shell-average, which we can solve to directly obtainx α (t). The equation is based on characteristic 2 of Sec. 1.1, which predicts that the shell-averages around all nodes, e.g., in K s (l) or in Q i (1), are all narrowly distributed and can, therefore be approximated by their ensemble average. This reduction is a direct consequence of the random nature of A and W , in which K s (l) and K s (l ) are statistically indistinguishable. While, generally, (2.60) represents an approximation, it is exact in the limit of a large (N → ∞) and perfectly random configuration model network, i.e. the framework discussed in Sec. 1. A detailed derivation of Eq. (2.60), including a discussion on its validity limits appears in Ref. 2 , and further generalized in Refs. [3][4][5][6][7] .
Extractingx α . Setting the derivative on the l.h.s. of (2.60) to zero, we obtain an algebraic equation forx α . This allows us to obtain the different fixed-points in function of the system's network parameters κ and ω. The stability of each potential fixed-point is then evaluated using linear stability analysis. This requires that a negative derivative, ensuring that deviations fromx α decay exponentially. Satisfying (2.62), for a given statex α , depends on the value of κ and ω through β, and hence on A and W . Equation (2.62), therefore, exposes the range of κ, ω values (if any) under which each fixed-pointx α is stable. Hence it exposes the naturally occurring states of the system, providing us with the state-space on which we apply our recoverability formalism. Specifically, it also helps us extract the undesiredx 0 , to be used in Eq. (2.50).

Option 3. Numerical analysis
The most accurate evaluation ofx 0 can be obtained from the actual observed state of the failed system. Indeed, before we apply our reigniting, we can use the detailed unperturbed state of the system x 0 = (x 0,1 , . . . , x 0,N ) , to accurately constructx 0 via (1.7). Such direct calculation provides an exact alternative to the mean-filed approximation offered in Option 2 above. It is especially useful if the network shows significant deviations from the configuration model. Yet, we emphasize, that our result clearly indicate that even imperfect networks continue to be well-approximated by Option 2.
Throughout our analysis we used Option 2, the mean-field approximation, to evaluatē x 0 analytically via Eqs. (2.60) and (2.62), and solve the recurrence relation to predict recoverability. The actual results forx 0 andx 1 , presented, e.g., in Fig. 5c of the main text were obtained from numerical simulations of a large scale-free network ensemble.  analyzing the intersection/s of F (x) and M 2 (x) we can predict the system's recoverability, and if recoverable, the required critical forcing.

Solving the recurrence relation
Cobweb plots (Fig. 2). We solve the recurrence relation of (2.49) using cobweb plots 9 . Starting from an initial setting of x s (0) = ∆, our maximal forcing, we track the evolution of the recurrence. First obtaining M s (∆) (vertical path), then shifting horizontally to F (x), extracting x s (1). Continuing the process we observe weather the recurrence converges to B 0 or to B 1 . In case the function F (x) is non-monotonic, this process may lead to convergence ambiguity, as illustrated in Case 3 of Fig. 2c, with both red and green pathways enabled. Such ambiguity is, of course, a mathematical artifact, as the real system in (1.1) will indeed follow only one of the potential tracks, not both. In reality, since we employ our reigniting from an initial condition which is in B 0 , the system, under any such instance of ambiguity, will converge back to B 0 , namely, it will select the red path and not recover. To remove this duality we use the construction of Fig. 2d, in which we introduce a plateau along the non-monotonic range in F (x). The corrected F (x) (solid purple line) is now monotonic and we can unambiguously analyze it via the proposed cobweb plots.

Dynamic models
To demonstrate our framework we examined the recoverability of four distinct dynamic systems, three within the form of Eq. (1.1), and on that extends beyond this form. Below we detail the analytical/numerical treatment of each of these systems, starting from the free system, in which we examine the states of the system absent our forcing ∆, then treating the reignited system, in which we introduce our single-node activation.

Cellular dynamics
We consider gene-regulatory dynamics, as captured by the Michaelis-Menten model 10,11 , for which (1.1) takes the form Under this framework M 0 (x i ) = −Bx a i , describing degradation (a = 1), dimerization (a = 2) or a more complex bio-chemical depletion process (fractional a), occurring at a rate B 12 . For simplicity, in our simulations we set B = 1. The activation interaction is captured by a Hill function of the form , a switch-like function that saturates to M 2 (x j ) → 1 for large x j , representing j's positive, albeit bounded, contribution to i's activity x i (t).
Next we test the dynamic stability of each of these three solutions via Eq. from which we can directly obtain the system's natural state-space. For the inactivex 0 = 0 we have the left hand side of (3.5) being −1, independently of β, and hence this state is unconditionally stable. The stability of the other two states is obtained by substitutingx 1 andx 2 into (3.5), finding thatx 2 is never stable, and hence it does not represent a potential state of the system, while the desirablex 1 is stable under the condition that β > 2 (Fig.  3b). Therefore β < 2 forces a collapse on to x 0 , while β > 2 represents a bi-stable state, where the system can reside in both x 0 and x 1 , depending on its initial conditions. (e) Increasing the average weight to ω = 1 the system transitions to Case 2 -recoverable if ∆ ≥ 0.67 (right), and unrecoverable otherwise (left). (f) -(i) The same analysis for the case h → ∞. In this limit the activation function M 2 (x) behaves as a step function. We find that here the structure of F (x), and hence the recoverability of the system, is independent of κ, affected solely by ω, ∆.

Reigniting
To examine the behavior of our cellular dynamics (3.1) under reigniting we seek to construct the recurrence relation (2.49), and specifically the function F (x) in (2.50). First we write where we used (2.39) to obtain R(x). Setting a = 1, we can now collect all the terms to construct F (x), providing us with a function whose shape depends on the topological parameters κ and ω. Here we used an undirected interaction pattern, setting ρ = 1 in (2.50). Also note that herex 0 = 0, and hence M 2 (x 0 ) = 0 on the r.h.s. of (2.50). Equation (3.9) maps to F (x) in Eq. (9) of the main text.
Varying κ and ω we can observe the recoverability of our system as F (x) transitions from Case 1 to Case 3 of Fig. 2. Two specific examples are presented in Fig. 3d and e. The first with κ = 5 and ω = 0.7 falls under Case 1, and therefore it is structurally unrecoverable, and the second, in which the average weight is increased to ω = 1 follows Case 3, i.e. recoverable if ∆ ≥ ∆ c .
To obtain the complete phase diagram and hence the boundaries of the recoverable phase, as we do in Fig. 3j-m of the main text, we systematically plot F (x) in (3.9) for a range of κ, ω values, seeking for each κ the critical ω in which F (x) transitions to the form of Case 3. These critical transition points provide the theoretical phase boundaries (Fig. 3km, white solid lines). At the same time we tested numerically, for each κ, ω combination whether single-node reigniting indeed reactivates the system. To achieve this, for each κ, ω combination we conducted 20 reigniting attempts with randomly selected source nodes, and plotted η, the fraction of successful recoveries (η = 0 yellow; η = 1 blue, see details in Sec. 5.3).

The role of the Hill coefficient
The Hill coefficient h in (3.1) determines the saturation rate of the activation function M 2 (x). A small h captures a mild activation, in which M 2 (x) increases gradually with x, while h → ∞ describes an effective step-function of the form being M 2 (x) = 1 (activation) if x > 1 and M 2 (x) = 0 otherwise; θ(x) is the Heaviside step-function. Taking this limit in (3.9) we obtain For ω < 1 we have a linear function whose slope is 1/ω > 1. This function has a single intersection with M 2 (x) = θ(x − 1) at x = 0, therefore rendering the system unrecoverable (Fig. 3h, left). In case ω ≥ 1 Eq. (3.11) becomes having two intersections with M 2 (x), describing a recoverable system for ∆ ≥ 1 (Fig. 3i, right).
This describes a limit in which κ plays no role in recoverability, and reigniting is driven solely by ω and ∆, as discussed in the main text under Restructuring guidelines.

Neuronal dynamics
We consider a modified Wilson-Cowan model 13,14 for excitation in neuronal networks, writing which we examine under µ = 10 and δ = 1. In this version of the model, adapted to the form of Eq. (1.1), the summation is extracted outside of the exponential function, namely we write N j=1 A ij W ij (1 + e µ−δx j ) −1 instead of (1 + e µ−δ N j=1 A ij W ij x j ) −1 . As a result, each node receives accumulating inputs from all its neighbors, and hence higher in-degree nodes gain a stronger overall activation signal from their surrounding Q i (1). This is in contrast to the standard version of the model, with the summation appears inside the exponent, and hence the effect of the multiple incoming signals quickly reaches saturation. In the present context, where we wish to observe the role of degree heterogeneity on reigniting, e.g., through κ, the form of Eq. (3.14) provides a more relevant testing ground.

Free system
To obtain the fixed-points of the system we use the mapping of (2.60) to reduce (3.14) into whose fixed-points follow Plottingx α vs. β (Fig. 4b) we obtain the dynamic phases of the system -the suppressed state x 0 (red), in which all activities are constricted, is obtained when the network is extremely sparse, i.e. small β; the active x 1 (green), in which x i are relatively high, is observed when β is large. In between these two extremes the system features a bi-stable state, in which both x 0 and x 1 are potentially stable. These phases are separated by two critical points β c,1 < β c,2 , predicting a hysteresis phenomenon: if β was driven below β c,1 and the system has failed, it will not spontaneously recover unless we retrieve β to be above β c,2 . Hence, we seek the sub-space within the bi-stable regime, β c,1 < β < β c,2 , in which the system can be reignited.

Reigniting
The relevant functions to construct F (x) (2.50) are from which we obtain, for an undirected network (ρ = 1) Once again we arrive at a function which depends on κ and ω, sometimes following Case 1, in which the system is structurally unrecoverable (Fig. 4d), and sometimes following the Case 3, in which it is potentially recoverable if ∆ ≥ ∆ c (Fig. 4e).

Microbial dynamics
We capture cooperative population dynamics via The self dynamics describes migration at a rate F coupled with logistic growth at rate B, with the system carrying capacity set to C, and the Allee effect 15 with strength K. The mutualistic interaction follows the Lotka-Volterra form x i x j . 16 In our simulations on scalefree networks (Fig. 5g-i) we set the parameters of (3.21) to F = 1/4, B = 1, C = 3 and K = 10; when simulating the microbiome we used F = 5, B = 3, C = 3 and K = 10. The reason for this difference is because for reigniting to be observed, we require a sufficiently broad bi-stable regime (Fig. 5b, grey shaded), prompting us to adjust the parameters, accordingly, as we implement (3.21) on vastly distinct networks. Specifically, setting F = 5, B = 3 in the case of our scale-free networks led to the system recovering spontaneously, as the active state was broadly stable, while the bi-stable window was narrow. This prohibited us from observing, in practice, the partitioning of the bi-stable state into two separate phases -nonrecoverable and recoverable. We therefore set these parameters lower to F = 1/4, B = 1, conditions under which the system has a sufficiently broad bi-stable state, for us to observe, with satisfactory resolution, its two sub-phases.

Reigniting
Here the functions comprising F (x) in (2.50) are and The challenge is that in order to construct F (x) we must invert R(x), which as indicated in Fig. 5c is non-monotonic and hence, in principle, non-invertible (purple solid line). The result is that R −1 (x) is ill-defined for a certain range of x, matching potentially two values for the same x. This ambiguity is directly related to the bi-stability ofx 0 andx 1 in the range β < β c . Indeed, Eq. (3.23) can be written in the form from which it follows thatx (3.27) in which the ambiguous value of R −1 (βx α ) for β < β c is precisely to root of the observed bi-stability. It, therefore, follows that of the two branches in R −1 (x) the relevant branch is the one associated withx 0 , as, indeed, the reigniting is applied on the failed state. Hence, to correct for this ambiguity we use a similar construction to the one shown in Fig. 2d when treating the non-monotonic F (x), namely, we eliminate the non-monotonic range of the function by replacing it with a constant plateau. This leads to the correctedR(x) shown in Fig. 5d, in which the original function (dashed-line) is replaced by a corrected monotonicR(x) (purple solid line). Its inverseR −1 (x) (yellow) is now well-defined, and, most importantly, suitable to predict the system's response to reigniting from a failed initial condition. We can now use this corrected function to construct F (x) in (2.50) providing This allows us to systematically use our cobweb plots to assess the system recoverability for different combinations of κ and ω (Fig. 5e,f). The resulting phase-diagrams in the κ, ω plane under different ∆ are also shown ( Fig. 5g-i). Note that for sufficiently large kappa, ω the system exits the bi-stable regime and enters a state in which x 1 is the only stable fixed-point (green phase at top right of phase-diagrams). In this state it spontaneously recovers even without reigniting.

Recoverability of the gut microbiome
As a realistic application of Eq. (3.21) we considered the recoverability of a failed microbial network, using the interaction topology of the human gut microbiome. The construction of this network and its link weights is discussed in detail below. The resulting topology is highly diverse, both in terms of P (k i,in , k i,out ) and in terms of P (w). Consequently, its shell statistics are non-uniform, with some shells characterized by higher degrees/weights than others. This represents a non-negligible discrepancy with respect to the configuration model ensemble, specifically violating characteristic 2 of Sec. 1.1. Moreover, the network comprises a combination of cooperative (W ij > 0) and adversarial links (W ij < 0), violating also characteristic 3. Hence, it offers an important testing ground for the broader applicability of our proposed framework.
Constructing the microbiome interaction network. To construct the gut microbiome we collected data on 766 types of bacteria, 53 archaea and 19 eukaryotes, together comprising N = 838 microbial species. 17 For each species the data include a description of all its consumed and produced nutrients among a set of M = 283 molecular compounds. This allows to construct two bipartite networks: the export network In case species n produces nutrient m 0 Otherwise , (3.29) and the import network These two networks encode information on the inter-species interactions. For example, if j produces m (E jm = 1) and i consumes m (I im = 1) then j positively contributes to i's growth rate, capturing a directed cooperative link j → i. In contrast if both i, j consume the same nutrient m (I im I jm = 1) then these two species compete over a shared resource, a bi-directional adversarial interaction i ↔ j. To quantify this we first normalize both matrices as The normalized E jm quantifies the relative contribution of j to the availability of the nutrient m, as a fraction of all other m-exports in E, namely N n=1 E nm . Similarly, I im captures the reliance of i on m, as one of i's total of M n=1 I in potential imports. Both terms become maximal in case j is the sole exporter of m, or if m represents a unique import of i. In contrast, if there are many alternative paths to produce m, or if m is but one of i's many consumed imports, the normalized weights in E jm and I im will be reduced, accordingly.
We can now use (3.31) to construct two N × N interaction networks. The first is the complementarity network (3.32) whose weights capture the positive contribution of species j to species i, namely P ij = P i←j . Indeed, j contributes to E jm percent of m's availability, which, in turn contributes to I im percent of i's growth rate. The second network captures the i, j competition level via quantifying i and j's mutual consumption of the shared resources. Note that, as opposed to P , Q is symmetric, hence Q ij = Q ji = Q i↔j .
Together P and Q allow us to construct the complete interaction network A and its weights W . For the topology we take A ij = 1 if either P ij or Q ij are non-zero. Namely A comprises all interactions, cooperative and adversarial. For the weights we assign a linear combination of the complementarity and competition networks. The rate constants ω P and ω Q represent a degree of freedom to set the average complementarity and competition rates. Indeed, P and Q are only designed to capture the relative complementarity or competition between all node pairs, namely which pairs have a stronger and which have a weaker interaction. They are not, however, suited to provide the specific rates, which depend on units (time −1 ), as well as on external conditions. For example, increasing the abundance of all nutrients is expected to weaken all competitive weights, as resources become affluent. Still it will preserve the relative strength of these weights, as expressed via (3.42), since these relationships are embedded in the intrinsic structure of the export and import networks of (3.29) and (3.30) -indeed, a characteristic of the chemical nature of all species and nutrients. Therefore, ω P and ω Q offer a set of tunable parameters, by which to rescale all weights in P and Q.
As a guideline for setting ω P , ω Q , we consider the survival of the network's N species. Indeed, as ω Q is increased, competition becomes dominant and many species undergo extinction. We, therefore, simulated (3.21), taking W from (3.34), under different values of ω P and ω Q , seeking the window of relevant complementarity/competition strengths under which the majority of species survive. This represents conditions, e.g., a level of resource availability, in which the system can potentially function. Of course, even under these conditions the system may still fail, precisely the scenario for our recoverability framework. Still, the parameters in (3.34) must, at least in principle, enable a viable active state, which constrains us to a bounded range of ω P vs. ω Q . We assume that these are the relevant conditions in which to examine the gut microbiome, as, indeed, under normal conditions, the majority of microbial species have non vanishing x i . In the simulations of the microbiome recoverability we used these guidelines to set ω P = 30 and ω Q = 1. The remaining parameters in (3.21) were set to F = 5, B = 3, C = 3 and K = 10. Under these conditions, the microbial network sustains a steady-state in which there are 568 active species. The remaining 270 species succumb to the adversarial links and undergo extinction.
Selective reigniting. Within the configuration model framework of Sec. 1 most nodes are surrounded by statistically similar shells. As a consequence the specific source node s is of little significance, and reigniting can be done with any selected node. In that sense, recoverability captures a characteristic of the system, e.g., the density/weights in the shells surrounding each node, and not of the specific choice of s itself. This is rooted in characteristic 2 of Sec. 1.1, stating that the K s (l) statistics are independent of s. Such uniformity is clearly observed in our simulations, for example in Fig. 3n-q of the main text, where recoverability rapidly jumps from η = 0 to η = 1. This indicates that, typically, reigniting either succeeds for all nodes, or fails for all nodes, with only a narrow window in which there is a preferred fraction of nodes capable of reigniting. Hence, unless the system is only borderline recoverable, there is no significant preference for one node or the other.
In the microbiome, however, K s (l) are diverse, and hence some nodes are better reigniters than others. To rank the reigniting potential of a specific node s, we characterized the shells K s (l) surrounding it. The idea is that for s to successfully activate the network its forcing signal must penetrate not just its immediate neighborhood, but also reach the distant shells. Therefore s's ranking is determined by all K s (l), rather than just by its degree.
To evaluate this, consider the process of reigniting. The reactivating signal originates at s, and propagates along the shells K s (l). A typical node i ∈ K s (l) receives this reactivating signal through its immediate neighbor j at K s (l − 1). Let us assume that j has been successfully reactivated, i.e., that x j is in a revived state x j = x Rev ∈ B 1 , then ask what are the conditions for this reactivation to sweep also through i. Note that we do not assume that x Rev =x 1 , sincex 1 represents the average state in K s (l) after the entire network has been revived. In the current scenario, reigniting is in the midst of its propagation, and most of j's neighbors, all in shells more distant than K s (l − 1), are still around B 0 .
Taking i to represent a typical node, its residual degree can be extracted from (2.20), and hence in addition to j, which is closer to s, node i is exposed to κ other nodes that are all farther away at K s (l + 1), K s (l + 2), . . . . At this time these other nodes have not yet been affected by the reigniting signal, and hence they are all still at the x 0 state. We can now approximate i's response to j's incoming signal via where we split i's κ + 1 interactions into two forms: κ neighbors with an average weight ω at state x 0 , and a single neighbor j at the state x Rev and with weight W ij . The question is, therefore, what is the minimal weight W ij that is sufficient to propagate the reigniting signal from j to revive i. where is the effective reviving force acting on x i . The resulting Eq. (3.36) tracks x i (t), as it interacts with a single activated neighbor j, and κ suppressed ones.
Node i will successfully reactivate if j's reviving forceW is sufficient to transition x i towards x 1 , namely if Eq. (3.36) no longer has a stable state around x i =x 0 , but only around x i =x 1 .
In Fig. 6b we display f (x) for different values ofW . The fixed-points of the equation are captured by the intersections at f (x) = 0, and their stability by the condition f (x) < 0, i.e. intersections where f (x) has a negative slope. For smallW we observe two stable states (red), one around x 0 and the other around x 1 . Under this condition, since x i 's initial state is in B 0 , i will not be revived, and revert to its pre-activation state. For largeW the only stable fixed-point is around the active state x 1 . This represents a successful activation, in which j was able to transfer the reigniting signal to i, i.e. the reviving force was sufficient (blue). In between, we observe a critical forceW = W 0 (grey), in which f (x) transitions from the bi-stable state to the single active fixed-point. This represents the minimal reviving force required for the reigniting at s to spillover from one shell (K s (l − 1)) to the next (K s (l)).
If, indeedW is at the critical W 0 , Eq. (3.36) predicts that x i will depart the inactive state, which is no longer stable, and transition to the revived state, x Rev . It will, in a sense, become the revived node itself, now propagating the reigniting signal further on to its neighbors in K s (l + 1). This new state of i is precisely captured by the intersection of f (x) (grey solid line) with the x-axis in Fig. 6b (green dot). Hence, this intersection point represents the revived state x Rev , which is an a priori unknown variable of (3.35). We can now graphically extract x Rev from the f (x) intersection, and together with the fact thatW = W 0 , obtain the critical weight as For j to successfully revive i, their link weight must satisfy W ij ≥ w 0 . This helps assess the reigniting capacity of a given node s: in case most links weights between K s (l − 1) and K s (l) are below w 0 , reigniting may fail to propagate to K s (l). If, however, there are many links with W ij ≥ w 0 , then K s (l) will be revived by s, as it received the activation signal from K s (l − 1).
The analysis above indicates that the prime contributors to reigniting are the links whose weights are greater than w 0 . If all shells are statistically identical, than they all include a similar density of such links, and hence they are all equally capable of propagating the reigniting signal (up to statistical discrepancies, which mainly occur at the bounds of recoverability, as seen, e.g., in Fig. 3n-q of the main text). If, however, the shells differ significantly, as is the case for the microbiome, then we must rank them by the number of nodes they have that can be reach via links whose weight W ij ≥ w 0 . To achieve this we first ⇓ (c) Step I 0 = 0 Rev ⇓ Figure 6: Extracting the critical weight w 0 and the revived state x Rev . (a) We track the propagation of the reigniting signal from node j ∈ K s (l − 1) to i ∈ K s (l). This propagation is mediated by the reviving forceW = W ij x Rev . (b) We seek a solution in which x i is revived, and hence must exit B 0 . This translates to two conditions: (i) f (x) has a single intersection with the x-axis in the range of B 1 (green); (ii) this intersection occurs at x = x Rev , as i is activated towards the state of its revived neighbor j.
Step I. First we seekW = W 0 that satisfies condition (i). This occurs atW ≥ 13.1, predicting this to be the critical reviving force W 0 (grey). (c) Step II. Once we set W 0 , we extract x Rev from the intersection of f (x) with the x-axis (green dot), fulfilling condition (ii). Here we find x Rev = 11.7, which using Eq. (3.39) predicts the critical reviving weight at w 0 = 1.12.
generate the weight-selective network in which we only keep the links with sufficient (positive) weight, equal or above w 0 . The selective network B includes only the links that actively contribute to spreading the reigniting signal -the effectual links. We can now use B to obtain the effective shells K B s (l). These shells, unlike K s (l), contain only nodes that can be reached via paths of length l form s, along the links of B. Such paths are constructed from sequences of effectual links, whose weights are all above w 0 , and hence nodes in K B s (l) represent potential target nodes that are prone to reigniting by s. This allows us to define, for each source node s, the reach which tracks the total number of nodes that s's reigniting signal can, in fact, reach. We predict that the ideal reigniters are nodes with the highest R s .
Simulating the microbial dynamics with the inter-species interaction network constructed above, we find that the critical reviving force isW = 13.1 (Fig. 6b, grey). The resulting f (x) intersects the x-axis at x Rev = 11.7 (green dot), which, using (3.39) predicts w 0 = 1.12 0 100 200 300 400 500 Rank 0 1 Figure 7: Selective reigniting in the microbiome. The reigniting success rate η vs. the R s -ranking as obtained from numerically simulating microbial dynamics. The top ranked reigniters (left) have a significantly higher probability to successfully revive the system. Highlighted (grey shade) are the top 26 reigniters, of whom 18 (∼ 70%) successfully revived the network. Among the remaining 543 species, only 2 (∼ 0.4%) were successful.
( Fig. 6c). Constructing the selective network B, we calculated the reach R s of all nodes and ranked them for their reigniting potential (Table 1, Supplementary dataset 1). We find, indeed, a highly non random distribution, in which the top 26 nodes have a significantly higher R s than all remaining nodes -a direct consequence of the uneven distribution of weights in the microbial interaction network. Therefore, we designate these top-reigniters to be our potential source nodes.
Restructuring. Our reigniting attempts with the top-ranked reigniters were mostly successful. We find that the top ranked 26 species, all characterized by a large R s , of order ∼ 300, (Table 1) were, by and large, capable of steering the network back to B 1 , as predicted (success rate: 18 out of 26, i.e. ∼ 70%). The remaining 543 surviving species have R s < 33, namely more than an order of magnitude below the top 26, and, therefore, quite expectantly, failed to reignite the microbiome (success rate: 2 out of 543, i.e. ∼ 0.4%). In Fig. 7 we present the reigniting success rate η vs. a node's R s ranking, showing, indeed, that all successful recovery instances are concentrated within the top ranking reigniters.
One of the failed nodes among the top ranked reigniters is the microbial species Pseudomonas putida (P. putida below), which, despite being ranked among the top 26, failed to revive the system. We, therefore, used this species as grounds for testing our restructuring strategy.
The observed reigniting failure of P. putida is rooted in the abundance of adversarial interactions, in which W ij < 0, that occur within its surrounding shells. These negative interactions quench the impact of the forcing, and hinder recoverability. It is, therefore, the natural directive to restructure the network by weakening or removing the adversarial links. This can be achieved via two types of intervention: 1. Node removal -via antibiotic treatment. Using targeted antibiotics to reduce or remove altogether the population of selected nodes. Targeting the nodes which contribute most to the competition, namely that strengthen Q in (3.42), will help systematically quench Q in (3.34), thus rendering the system more recoverable.
2. Link deletion -via nutritional intervention. Here we reduce the weight, or even delete, selected competitive links, by eliminating competition over specific nutrients. For example, if, using nutritional supplements, we generate an abundance of nutrient n, then the competition over this resource is quenched. In technical terms, such intervention replaces Q in (3.42) by where the contribution of n to the weight of the competitive (i, j) link is eliminated from the sum.
In our implementation we used option 2, which represents a more nuanced intervention, that we believe is more likely to be used in the face of a dysfunctional microbiome. Indeed, antibiotics, are, often the root cause of microbiome failure, [18][19][20] and removing additional species, from an already diluted microbiome may be deemed risky. Hence we ranked all nutrients m = 1, . . . , M based on their overall contribution to the competition matrix, as evaluated by capturing the percentage reduction in the overall competition that originates in nutrient m.
The resulting ranking is shown in Table 2, with the top three nutrients contributing to Q identified as Glucose, Pantothenic acid (vitamin B5) and Nitrate. 21 Ensuring an abundance of these three nutrients we reconstructed Q and W in (3.34), eliminating the competition over these selected modes. We then attempted again to reignite the failed microbiome by forcing P. putida. We observe that following these, quite accessible, restructuring interventions, the microbiome, indeed, becomes recoverable through the P. putida-reigniting.

Diffusive dynamics
We consider a system of coupled units (spins) subject to an external force (magnetic field ), as modeled through the Ginzburg-Landau Equation.
Here, the spin dynamics are captured via a continuum approximation, following 22 The first term on the r.h.s. represents the spin self-dynamics, alternating between the two spin-states x i = −Q and x i = P . We focus on the case P > Q, capturing the breaking of symmetry due to the external field, which we arbitrarily take to be in the positive direction. The interaction term, diffusive in form, captures the magnetic pull between neighboring spins, seeking to align them towards x i = x j . The weighted topology A ⊗ W describes the adjacency patterns and strength of the spin-spin interactions. We consider initial conditions where all the spins are in the negative state x i = −Q. Under a small external field, this state remains stable. We then select a random source spin s and force it to sustain a positive state x s (t) = ∆, aiming to flip the entire network towards x i = P . This is analogous to single-node reigniting, transitioning the state of the entire network by pinning the activity of just one node.
Note that Eq. (3.44) does not fall under the framework of Eq. (1.1), as, indeed, the interaction term cannot be factorized into the form M 1 (x i )M 2 (x j ). We, therefore, explore this dynamics numerically, aiming to examine the potential breadth of our proposed analysis, and its applicability beyond the limits of our analytical framework of Sec. 1.
In Fig. 8 we show the numerically obtained phase diagram in the κ, ω space, under a fixed ∆ = 5. Similarly to all other examined systems, also here we find a recoverable phase (blue) in which the s-forcing was, indeed, able to propagate to the entire network. Interestingly, in this system, recoverability favors lower κ, and hence the sparser the network A, the more recoverable it becomes. This is a unique feature of the diffusive dynamics in (3.44), which is attractive (x j − x i ) rather than cooperative (e.g., x i x j ). In such dynamics x i tends to match with x j , rather than simply grow with j's activation. Therefore, having many neighbors in the Q-state makes flipping the network towards P more difficult, as they all tend to oppose our reigniting force at s.

Practical reigniting
Our mathematical formulation focuses on a time-independent forcing of the form x s (t) = ∆. This allows a tractable analysis, in which the forcing acts as an effective boundary condition, allowing us to analytically predict the long term response of the system via the recurrence (2.49). In reality, it is often difficult to sustain such a constant reigniting signal, especially as the source node s continues to interact with all other nodes via (1.1). The crucial point is, however, that recoverability is, by and large, insensitive to the specific form of the forcing, as long as the reigniting node is kept at sufficiently high activity for a sufficiently long period. Indeed, an active s reactivates its surrounding shells K s (l) even if its activity x s (t) fluctuates, as long as it maintains an average high activity for enough time to drive all K s (l) into B 1 .
To examine this we consider a practical forcing, in which at given time intervals T we externally administer jolts to the activity of s as with n = 0, 1, . . . . This can capture, for example, the periodic intake of medication or probiotics, which creates an instantaneous boost in the concentration of a specific protein, metabolite or microbial species. In neuronal or ferromagnetic networks, such forcing represents repetitive excitation, which can help gradually reactivate the entire system. Following each boost, x s (t) will, at first, be highly active, but then undergo relaxation back towards B 0 , driven by the natural system dynamics (1.1). Let us consider x s (t)'s natural relaxation time τ , after which the effect of the ∆-boost wares off. While this is not guaranteed, we can treat, for simplicity, x s (t)'s relaxation following the nth reigniting boost as an exponential decay of the form Other forms of decay are also possible, depending on the structure of Eq. (1.1), but the analysis that follows remains equally relevant. For the periodic reigniting (4.1) to succeed we require that x s (t) remains on average above the critical forcing ∆ c . Using (4.2) this translates to which in turn predicts  Figure 9: Practical reigniting. (a) We applied reigniting through periodic boosts a là Eq. (4.1) to our cellular dynamics. The network has κ = 5 and ω = 0.8, for which our theory predicts a critical reigniting force of ∆ c = 1. (b) Using boosts of magnitude ∆ = 2 and periodicity of T = 2 we observe a successful reigniting. The source node (black) exhibits a chainsaw dynamics, with each boost followed by the natural system relaxation, hence it cannot sustain a constant forcing ∆. Still, the boosts are sufficiently frequent to gradually reactivate all other node activities x i (colored solid lines), successfully steering the system into B 1 . (c) Increasing the lag between successive boosts to T = 4, we now witness a failed reigniting, as predicted.
required activation boosts (large T ). On the other hand, if ∆ c is large, i.e. the system is difficult to reignite, we need more frequent instances of activation. This analysis can be readily generalized for different relaxation functions, other than the exponential decay of (4.2).
In Fig. 9 we examine this reigniting strategy on our Cellular dynamics of Sec. 3.1. We evaluate the relaxation time τ via the nodes' self-dynamics, writing dx s dt ≈ −Bx a s (t) + . . . .

(4.5)
The idea is that the interaction term in (1.1) contributes little to s's relaxation as compares to its own internal dynamics. Setting B = 1 and a = 1, as we do in Sec. 3.1 we have in Eq. Setting ∆ = 2, i.e. boosts that are roughly double in size of the critical forcing, this predicts T 2. Indeed, we find that reigniting via (4.1) with ∆ = 2 and intervals of T = 2 successfully revives the system (Fig. 9b). However, similar conditions with T = 4 fail to reignite the system, as the boosts taper off before the system exhibits a sufficient response (Fig. 9d).
In practice, Eq. (4.4) provides an underestimation of the maximal lag T . This is because, as the reigniting boosts accumulate, the surrounding nodes increase their activity, providing feedback that further reinforces the system's recovery. Therefore, while a lag significantly above T = 2 may fail, as shown in Fig. 9d, if T is just slightly above the predicted threshold, we may still witness a successful recovery. Indeed, in Fig. 9c we observe a revived system with T = 3, slightly above the predicted threshold of 2.

Numerical integration
To numerically test our predictions we constructed Eq. (1.1) for each of the systems in Sec. 3, using the appropriate network A (Scale-free, Erdős-Rényi, empirical, etc.) and weights (random P (w) or empirically constructed). We then used a fourth-order Runge-Kutta stepper (Matlab's ode45) to numerically solve the resulting equations. Starting from a pre-selected initial condition x i (t = 0), i = 1, . . . , N we allowed the system to reach steady-state by waiting forẋ i → 0. To numerically realize this limit we implemented the termination condition where t n is the time stamp of the nth Runge-Kutta step and ∆t n = t n − t n−1 . As the system approaches a steady-state, the activities x i (t n ) become almost independent of time, and the numerical derivativeẋ i = (x i (t n ) − x i (t n−1 ))/∆t n becomes small. The condition (5.1) guarantees that the maximum ofẋ i over all activities x i (t n ) is smaller than the pre-defined termination variable ε. In our simulations, across the different dynamics we tested, we set ε ≤ 10 −2 to ensure that our system is sufficiently close to its true steady-state.
In case of bi-stability we examined the convergence of the system from multiple initial conditions. For example, setting x(t = 0) to a low value in B 0 ensures convergence to x 0 , in case x 0 is stable; setting it in B 1 ensures convergence to x 1 , in case x 1 is also stable. If only one of the states is stable -all initial conditions will converge to that single fixed-point.

Reigniting
To simulate reigniting we set the initial condition of the system to x(t = 0) = x 0 . We then select a random node s, decouple it from the remaining N equations and set its state to x s (t) = ∆. Together with the remaining N − 1 equations of (1.1) we arrive at Eq. (2.1), which takes the form Integrating this equation until reaching steady-state, i.e. condition (5.1), we find the final state x Forced of the forced system. We then relax our forcing, re-couple x s (t) to the remaining N − 1 equations, and allow the system to relax to its final (unforced) state. This is achieved by setting the new initial condition to x(t = 0) = x Forced , and numerically solving Eq. (1.1) until reaching steady-state. In case x Forced ∈ B 1 , a successful reigniting, the system will reach x 1 . If, however, our reigniting failed, and x Forced remains in B 0 , the system, after forcing ceases, will revert to the undesired x 0 .
To simulate the forcing via periodic boosts, as considered in Sec. 4, we simply replace the equation x s (t) = ∆ in (5.2) with the dynamic forcing of Eq. (4.1). Sustaining this dynamic activation for ∼ 20 − 30 cycles, we once again cease our forcing and observe the system convergence to its natural state -x 1 , successful, or x 0 , unsuccessful.

Constructing the phase diagrams in Figs. 3 and 5
We used the configuration model framework, described in Sec. 1, to construct undirected networks with a pre-assigned κ and omega. The networks all included N = 10 4 nodes, had a scale-free degree distribution P (k) ∼ (k − k 0 ) −γ and a bounded weight distribution P (w). This allowed us to control κ by varying k 0 within the range [1,3] and γ within [2.5, 4]. Together we generated networks corresponding to 50 distinct values of κ, ranging from ∼ 5 to ∼ 50 (Fig. 3) or ∼ 3 to ∼ 30 (Fig. 5). Next we matched each of these networks network with weights, representing 50 distinct values of ω, providing us, together, with 2, 500 independent combinations of κ and ω, upon which to examine our reigniting analysis. Each such κ, ω combination represents a single pixel in the κ, ω phase diagram. For each κ, ω combination we constructed 20 independent network realizations, i.e. 2, 500 × 20 = 5 × 10 4 weighted scale-free networks in each phase-diagram, comprising 20 networks per pixel.
For each network, we obtained the steady-state/s as explained in Sec. 5.1, observing whether they are in the active, inactive or bi-stable regime. We then also tested whether they are recoverable via single-node reigniting, following Sec. 5.2. Counting the successful reigniting attempts C among the 20 networks with a given κ, ω, we extracted η = C/20, the fraction of successful reigniting instances in each pixel. In a similar fashion we also selected 50 values of the forcing ∆, providing us, once again with 2, 500 combinations of κ, ∆ and an additional 2, 500 of ω, ∆, upon which to test reigniting. As above, each pixel, represents 20 realizations of reigniting, here selecting a random source node.

Model and empirical networks
We used model and real networks, as summarized below: Scale-free. Scale-free networks, constructed via the configuration model, with N = 10 4 nodes, degree distribution P (k) ∼ (k − k 0 ) −γ . Varying the parameters k 0 and γ we obtain an ensemble of scale-free networks, which we used to construct the phase-diagrams of Figs. 3 and 5 of the main text.
Brain (Neuronal). Mapping the physical fiber bundle connections between 998 brain regions, as measured using diffusion tensor imaging techniques 25 . The empirical network has a very broad weight distribution P (w), ranging over several orders of magnitude, as 3 × 10 −5 < w < 0.23. To extract the meaningful connections, we set a threshold at w 0 = 0.03, and removed all links (i, j) with W ij < w 0 , resulting in a network constructed from the top 12.5% strongest links. In this network the link weight W ij does not represent the interaction rate between i and j, but rather our confidence in the observed (i, j) link. Therefore, once we selcted the top-scoring links, their weight is irrelevant in Eq. (1.1), and hence we set it to be uniform.
This empirical network has a natural modular structure owing to the brain's two hemispheresthis allowed us to investigate the impact of molecularity on reigniting, by setting distinct weights for the commiseral vs. the intra-hemispheral links. The results are presented in Fig Table 4: Robustness of our analytical framework. For each of our empirical networks we present the clustering C, the normalized clustering c and the degree-correlation coefficient r. Note that while the gutmicrobiome is a directed network, here, for simplicity, we extracted C and r from its undirected equivalent A UD ij = max(A ij , A ij ). This is because our only goal is to evaluate the level of discrepancy of our real-world A from the configuration model, not to assess a specific characteristic of this or another networks.

of the main text.
Gut microbiome (Microbial). This network construction is detailed in Sec. 3.3.3.

Deviations from the configuration model
While our analytical derivations are suited for the configuration model ensemble of Sec. 1, in practice, we find that they offer rather accurate predictions and qualitative insights also on real-world networks, all of which exhibit deviations from the configuration model. Specifically, as explained in Sec. 1.1, there are two features of the configuration model that are especially relevant in our analysis: (i) The tree-like structure, implying a scarcity of loops; (ii) the shell statistics, indicating the random nature of the connectivity patterns in A.
We can quantify the extent to which a network violates (i), by measuring its clustering coefficient C 26 . This measure, indeed, tracks the abundance of short range loops in the network. For a randomly wired network we expect C 0 = k /N , which approaches zero in the limit N → ∞. In Table 4 we show, for each network its actual clustering coefficient C, and also its normalized clustering c = C/C 0 . The latter helps us quantify the level of deviation from our assumed tree-like network. We find that our empirical networks show significant discrepancies from the random wiring assumption, with C often an order of magnitude above C 0 , and yet our predictions continue to hold. This indicates that our theory is, by and large, robust against the presence of loops. In the case of Brain (c = 75.6), the deviation is truly consequential. And, as we have shown in the main text, it has a distinctive modular structure, which, indeed, requires specific treatment (Fig. 5 of main text).
Characteristic (ii) can be examined via the network's degree-correlation coefficient r, designed to capture the distinct statistical properties of the shells surrounding different degree nodes. Once again, we expect, in the configuration model ensemble, to have r → 0. In reality, our networks exhibit non-negligible r, ranging from −0.48 to +0.37, an array of non-negligible deviations from characteristic 2 of Sec. 1.1. This, again, is a testament to our framework's robustness against such microscopic and meso-scopic discrepancies.