Dynamical robustness of collective neuronal activity upon targeted damage in interdependent networks

In the last decades, the availability of data about the structure of social, technological and biological systems has provided important insights on the mechanisms governing their correct functioning and robustness. These mechanisms are grounded on the complex backbone of interactions among the constituents of the system, which include both topological and dynamical aspects. Here, we analyze interdependent networks composed of two layers of interacting neuronal units and explore their robustness when these synthetic cultures are subjected to damage in the form of either targeted attack or failure. Our results show that the functionality of these networks does not decrease monotonically with damage but, on the contrary, they are able to increase their level of activity when the experienced damage is sufficiently strong.


Introduction
Naturally-formed neuronal circuits, from the worm C. elegans up to the human brain, have an anatomical structure characterized by a modular organization, in the sense that groups of neurons (from few dozens to millions) tend to connect more strongly within their group than with the rest of the network [1,2]. In the mammalian brain, for instance, the visual cortex and the sensory-motor area shape two important modular units, which operate as specialized circuits and that share information when required with other modules and the rest of the brain. In addition, these modules exist in a physical environment that combines two-dimensional and three-dimensional characteristics, with inherited developmental cues and wiring cost constraints [3,4]. Although this modular organization is central for information processing and functional performance, the topology of the modules themselves and their interconnectivity is pivotal to shape a flexible yet robust neuronal system [5]. Flexibility is required to balance specialization and whole-brain integration, while robustness is required to ensure that perturbaa e-mail: sergiofacilazaro@gmail.com (corresponding author) tions or damage do not critically compromise the operability of the circuits.
To understand the robustness of modular neuronal circuits, some important aspects to take into account are node centrality and network topology as well as their interrelation with intrinsic neuronal dynamics and noise. Given the inherent complexity of neuronal systems, studies in the brain [5,6] have been combined with in vitro experiments in neuronal cultures [7][8][9][10] and numerical simulations [11,12]. These explorations aim at investigating the effect of damage under two main scenarios: (i) targeted attacks, i.e., the deletion of central neurons, and (ii) failures, i.e., the deletion of randomly selected neurons. With these approaches it is possible to explore the behavior of neuronal networks as they experience the removal of units, i.e., as damage progresses, and provide important information to comprehend and model neurological disorders [13].
Recently, Faci-Lázaro et al. implemented a biologically-realistic spatial model for the construction of synthetic neuronal cultures and explored numerically their response to attacks and failures [12]. One of the main results observed is that metric correlations between the neurons, caused by the spatial embedding of the culture, are closely linked to noise-amplification mechanisms that trigger spontaneous activity [14]. Moreover, as the damage in the culture progressed, they found a transition in the network dynamics from wholeculture to fragmented. This transition was related to a fall-and-rise behavior of network activity level with damage, a surprising and apparently counter-intuitive phenomenon.
In the present work, we aim to acquire a better understanding of the mechanisms that underlie the former transition and the influence of node centrality in the behavior of a synthetic modular culture. Our work is inspired by the intricacy of brain modularity and by recent experiments in vitro by Yamamoto et al. [7] and Barral et al. [8], who engineered small modular or layered neuronal cultures and studied the collective dynamics of the system for different coupling strengths among neuronal populations. To this aim, instead of a biologically realistic model, we build much simpler cultures by using standard synthetic graph models, in particular the Erdös-Rényi (ER) and the Barabási-Albert (BA) models, to represent the interconnections between neurons. Thus, by partially renouncing the accuracy of biologically inspired models, we gain the physical intuition provided by simple graph models with reduced complexity, and whose properties are well known and characterized in the literature [15,16].
The article continues as follows. In Sect. 2 we introduce the structure of the synthetic modular networks and the way these graphs are combined to form simple interdependent networks with tunable correlations. In Sect. 3 we describe the neuronal dynamic model which, as in Ref. [12], follows the Izhikevich description. In Sect. 4 we study the behavior of the synthetic cultures when subjected to different damage processes and, finally, in Sect. 5 we summarize the obtained results.

Structural model
Our in silico interdependent neuronal cultures were designed following a multilayer model with L = 2 layers, where each layer represents a module in a neuronal network. Each layer contained N = 250 nodes and was built as a directed graph, following either an Erdös-Rényi (ER) [17,18] or a Bárabasi-Albert (BA) [19] topology. We considered the situation in which both layers in the multilayer construction had the same topology, so that the resulting multilayer networks were either ER-ER or BA-BA, with a total of 2N nodes. To achieve this, the connection probability p of the ER model, as well as the number of initial connections for a new node m of the BA model, were set to p = k /N and m = k . Here k is the average connectivity in each layer. The layers were interconnected by adding f · N 2 directed links between them, where f is the desired fraction of interlayer links.
The nodes from which each interlayer link originated and terminated were set following a probability distribution, constructed as proportional to the intralayer degree of each node elevated to a certain exponent β , The two networks have the same intralayer connectivity, but different interlayer one. The one on the left is constructed with a random interlayer connectivity scenario, with β1 = β2 = 0, whereas the one on the right is constructed with the hub-to-hub interlayer connectivity scenario, with β1 = β2 = 1. Circles represent nodes (neurons), and their darkness is proportional to their degree. Intralayer links are plotted in black and interlayer links in red. In this representation it can be observed that for the β1 = β2 = 0 the generation of an interlayer link between two nodes is independent of nodes' degree, whereas for β1 = β2 = 1 the interlayer links are preferentially formed between high degree nodes with indicating the layer, i.e. = 1, 2. This probability distribution then reads: Thus, β = 0 would correspond to the case in which the nodes from the -layer are chosen at random as sources or destinations of links, whereas β > 0 (alternatively β < 0) would give preference to nodes with higher (lower) values of intralayer degree as sources or destinations. Two possible scenarios for the interconnectivity ER-ER or BA-BA can be considered depending on the values of β 1 and β 2 : (i) β 1 = β 2 = 0, so that the nodes that shape the links between layers are chosen at random in both layers. (ii) β 1 = β 2 = 1, so that the nodes of interlayer links are most likely those with the highest degree, i.e. the hubs in each layer will tend to connect with each other.
A schematic representation of the networks obtained with this model and the above discussed values of β is shown in Fig. (1). The details of the parameters used in the construction of the structural networks are provided in Table 1.

Dynamical model
The dynamics of the neurons were simulated through the Izhikevich model [20], which is an efficient and biologically plausible model for cortical spiking neurons [21]. Following the implementation described in Refs. [12,14], the generation of action potentials in a neuron soma is constructed as an integrate-and-fire system with adaptation and governed by: where v j and u j are the membrane potential and the recovery current of the j-th neuron; τ c is the leaky capacitance; v r and v t are the resting and threshold potentials, respectively; I j is the synaptic current the j−th neuron received from the rest of the network; τ a is the main time scale of the inhibitory currents; v p is the value of the membrane potential at which a synaptic pulse is emitted; v c is the reset value of the membrane potential after spike emission; d accounts for the adaptation and recovery of the neurons; R and b are, respectively, the growth rate of the membrane potential and the neurons' sensitivity to subthreshold fluctuations; and ξ is the neuronal network noise. Two contributions are taken into account: a Gaussian noise η with η(t)η(t ) = 2g s δ(t − t ), associated to fluctuations in the membrane potential; and a shot noise which accounts for random synaptic pulses elicited by the neurons without external excitation. These pulses are called minis, occur with a frequency λ and have a synaptic strength g m that decays with a time constant τ m . In this model the spikes generated by each neuron are transmitted to their neighbors as synaptic currents of the form: where A ij is the connectivity matrix of the network and E i (t, t m ) is the current induced by the ith neuron at a time t as a result of the spike generated at time t m . The variables g A and τ A are the strength and decay time of the synaptic current, respectively. D i (t) is the depression term that describes the efficacy of the neuron presynaptic terminals [20,23]. The depression has a resting value of 1 and relaxes exponentially with a decay time of τ D . β is a coefficient related to the loss in efficiency that occurs whenever a synaptic pulse is generated. Table 1 summarizes the parameters used in the construction of the numerical model. A given network with 2N neurons was simulated for about 30 s, which provided sufficient statistics for spontaneous neuronal activity, which was typically characterized by episodes of whole-network activity combined with quiescent intervals. For a given activity episode we quantified the avalanche size S and the degree of activity A [12]. The avalanche size S was defined as the fraction of the neurons in the network that co-activated together. An avalanche was accepted as such whenever any two neurons had at most a delay of 50 ms in their activations. A larger time automatically separated an avalanche from another. The degree of activity A was the number of observed avalanches in the 30 s of simulation. A total of 10 network realizations of either ER-ER or BA-BA were explored. The variability among realizations was quantified through the standard deviation of the mean for both S and A.

Results: dynamical response to damage
Once the cultures were constructed and their dynamics characterized, we adjusted the parameters to ensure that any culture -pair of interconnected layersoperated in a synchronous manner, i.e., that all 2N neurons activated together in a short time window. In this situation, we studied their response to damage caused by sequential neuronal deletion, and where the connectivity k of each neuron, given by k = k intra + k inter , was the important measure to decide whether a neuron was a central node or not. Indeed, to assess the dynamical response of the cultures to damage, we compared the results obtained when removing nodes in decreasing order of degree k (attacks) with those obtained when removing the neurons at random (failures). We also studied the effects of the degree correlation between the two network layers and the influence of the degree distributions of each layer in the response to damage.  . 2 Size of the giant component G, avalanche size S and degree of activity A upon targeted attack for single layer cultures with ER and BA topologies. The panels show the evolution of G, S and A as a function of the fraction of removed nodes in the network q following a node deletion from higher to lower connectivity k, comparing BA networks (blue) and ER (red) ones. To better compare the ER and BA configurations, the avalanche size S is given as the fraction of neurons in the network that activate together, while A is scaled relative to the activity value at q = 0 (17 avalanches in 30 s for ER, 20 avalanches in 30 s for BA). Data is averaged over n = 10 numerical realizations of each topology, with the shading corresponding to the standard deviation of the mean. Both networks show an increase in activity for q ≈ 0.8. For the ER networks, this increase in activity is accompanied by a slight drop in avalanche size while the giant component remains G ≈ 1. Instead, for BA networks, the avalanche size remains large while the activity increases even though the size of the giant component starts decreasing for q ≈ 0.85, eventually reaching G = 0 for q ≈ 0.95 For the sake of clarity, before analysing the dynamical response to damage in two-layer cultures let us show the behavior of monolayer cultures, i.e. single ER or BA networks. In Fig. 2 we show the evolution as a function of the fraction of removed neurons, q, of three different metrics that capture the structural and dynamical evolution, namely the giant component of the network G, the avalanche size S and the degree of network activity A. The two latter (dynamical) metrics appear normalized to their corresponding values at q = 0, i.e. to those of the original culture, while in the case of G we compute the number of nodes in the giant connected component divided by the total number of nodes with degree k i > 0. In this way we avoid the trivial decrease of G due to the removal of nodes while, at the same time, we easily capture potential giant component break down events.
From the panels it is clear that the dynamical response of both networks remains unaltered as q increases until large values of q are reached. In particular, a fraction of neurons q ≈ 0.7 should be removed to observe a substantial change in the activity of the network that, for both ER and BA networks, experiments an increase while maintaining the same average avalanche size S. It is important to note that the boost of the dynamical activity of ER and BA networks occurs for G ≈ 1, i.e. it is not related to a network fragmentation that, in its turn, occurs for much larger values of q (q ≈ 0.9). Summing up, for both ER and BA monolayer networks, the impact of structural damage on network behavior appears before the giant component breaks down; an impact that is reflected in the increase of network activity up to 50% of the original network (q = 0). This increase has a non-trivial nature associated to the complex interplay between connectivity, intrinsic neuronal activity and noise. Now we put our focus on the two-layer cultures. Let us start with the case in which there are no correlation between the degrees of the neurons that are coupled through inter-layer links. In this case, the probability of forming an interlayer link is characterized by β 1 = β 2 = 0. In Fig. 3a we show again the evolution of normalized values of the giant component size G (top), the avalanche size S (middle), and the degree of activity A (bottom) as a function of the fraction of removed neurons, q, for ER-ER (red) and BA-BA (blue) cultures. Again, as in the case of monolayer cultures, the dynamical changes observed in S and A are not correlated with a break down of the network in different connected components, since G remains close to 1 until q reaches very large values (q ≈ 0.95). However, the dynamical response of BA and ER bilayers is different, and richer, than in the case of monolayers. On the one hand, for ER-ER we observe a similar pattern to that shown in Fig. 2, i.e., an increase of network activity, although here we observe two well defined peaks instead of a single maximum. On the other, the BA-BA behavior is clearly different form the case of the single BA network, for it shows a significant drop in the activity pinpointed by two local minima of A(q) and that take place at q ≈ 0.5 and q ≈ 0.8. After these two drops, and before the BA-BA duplex breaks apart, the level of activity recovers to values close to the original undamaged network, although the avalanche size is much smaller (S ≈ 0.7). It is important to note that, in this regime of large values of q, the large variability between BA-BA culture realizations indicates that this system is very sensitive to the wiring among neurons. Now we turn our attention to the case of correlated two-layer cultures. When the correlation between interlayer nodes' degrees is introduced by using β 1 = β 2 = 1 [see Fig. 3b], the former observations become more noticeable for large values of q, specially for the ER-ER cultures. In particular, when the correlated cultures experience damage through targeted attacks, a significant burst in activity is found at q ≈ 0.75 for ER-ER cultures and at q ≈ 0.8 for BA-BA ones. The only difference between these two network designs is that, in the ER-ER cultures, this burst is accompanied by a substantial decrease in the size of the avalanches, whereas this does not occur for the BA-BA cultures. The same behavior can also be found for the case of failures (figure not shown), but they are not as noticeable and take place at higher values of the fraction of removed neurons, q ≈ 0.95.
We note that, as in the case of monolayer and uncorrelated two-layer cultures, the dynamical effects in response to damage emerge much before the duplex breaks down. In fact, the change in the degree of activity of the ER-ER cultures pinpoints that, as these cultures experience damage, delivered through either targeted attacks or failures, they eventually reach a point in which their activity transitions from encompassing the whole-network to a fragmented dynamical state, in which there are more events of coherent activity but involving a smaller fraction of the remaining neurons of the culture, as also observed in Ref. [12].
In the case of BA-BA cultures, this transition to coherent but fragmented dynamics is not observed, since the size of the avalanches remains close to 100% of the network up to q → 1. However, this is to be expected, since the main characteristics of BA networks are the presence of hubs and a high clustering among these highly connected nodes. These topological features of BA networks act as amplification mechanisms of the activity of the culture [22] by facilitating the existence of feed-forward and feed-backward loops [14,24].
The observation that BA-BA cultures do not fracture regardless the level of damage points out that these cultures are extremely resilient and robust. This is, BA-BA cultures are able to keep the main topological properties responsible for the presence of wholenetwork coherent activations, as opposed to the ER-ER cultures, which show a behavior closer to that of cultures with biologically realistic connectivity patterns [12]. This can be understood by taking into account that small enough realistic cultures display a structure of connections similar to those present in ER graphs [24].
It is also important to note that, unlike what is observed in realistic neuronal models [12], Fig. 3 shows that the degree of activity does not approach zero at any point, especially in the ER-ER cultures. For this to occur, the structure would have to be such that it would be possible to eliminate (almost) all avalanche precursor nodes, and without any other node replacing them and taking up their role before the transition to fractured dynamics occurred. In other words, the structure would have to be hierarchical, which is a property that the simple BA and ER graph models lack. Hierarchical organization is a feature that has already been proposed as a trait of great importance in biological systems [7,25]. Thus, topological hierarchy would be a To better compare the ER-ER and BA-BA configurations, the avalanche size S is given as the fraction of neurons in the network that activate together, while A is scaled relative to the activity value at q = 0. a Interdependent network built using a random interlayer connectivity (β1 = β2 = 0). The panels show the evolution of the different metrics as a function of the fraction of removed nodes in the network q following a node deletion from higher to lower connectivity k, comparing BA-BA networks (blue) and ER-ER (red) ones. Activity at q = 0 was 17 avalanches in 30 s for ER-ER and 20 avalanches in 30 s for BA-BA. b Equivalent exploration for a network with hub-to-hub interlayer connectivity (β1 = β2 = 1). Activity at q = 0 was 16 avalanches in 30 s for ER-ER and 19 avalanches in 30 s for BA-BA. In all panels data is averaged over n = 10 numerical realizations of each topology, with the shading corresponding to the standard deviation of the mean. ER-ER networks show a substantial increase in activity for q ≈ 0.8 in both interlayer configurations, which is accompanied by a drop in avalanche size. BA-BA networks show strong fluctuations in the random interlayer configuration, but activity remains similar along q within error bars. These networks substantially increase in activity for in q ≈ 0.8 in the hub-to-hub interlayer configuration, but this increase is not accompanied by a decrease in avalanche size. It is important to note that all these phenomena occur with G ≈ 1 fundamental ingredient when trying to develop accurate models aimed at simulating or predicting the response to damage in biological neuronal networks.

Discussion
Recognizing the basic structural properties of neuronal networks is fundamental to fully understand their behavior. To this aim, in this work we have used two of the most studied and well characterized models to built our synthetic cultures, namely the Erdös-Rényi and the Barabási-Albert models. With them, we created synthetic bilayer cultures of ER-ER or BA-BA, intercon-nected at random or by considering degree correlations, and studied their response to attacks and failures. We observed that once the culture has experienced enough damage (q ≈ 0.8 for targeted attacks and q ≈ 0.95 for failures) it experiments a burst in its activity as well as a decrease in the average avalanche size. This indicates that, at that level of damage, the dynamics of the culture transitions from whole-network to fractured, a behavior that is consistent with previous results [12]. In addition, we also observed that spontaneous activity never reached zero in the synthetic networks studied, i.e., they are always capable of sustaining a substantial level of activity. This means that, regardless the amount of avalanche initiating neurons that were removed from the culture upon damage, there was always other nodes able to take up their role.
With this in mind, a first ingredient that a numerical model should have to capture the dynamics of living neuronal systems -particularly in the context of damage-is a broad degree distribution, so that it is sufficiently heterogeneous for the culture to have a fraction of neurons with the necessary topological properties to act as initiators of coherent activity. Some of these properties include the presence of feed-forward and feed-backward loops [14], which are strongly related to high clustering coefficients and connectivity. A second ingredient is that the model should shape a hierarchical structure, as otherwise the culture activity would remain somewhat constant until the transition to a dynamic fracture state. Interestingly, experiments in neuronal cultures in which nodes were randomly removed (the equivalent to failure), showed an important drop of activity, by 60%, following damage on just a single node [9]. This illustrates how important is the development of numerical models to understand which ingredients are necessary for modeling experimental behavior.
In a nutshell, the results obtained in this work point out that, although some of the properties of simple graphs can explain the dynamic behavior and robustness observed in real neural cultures, the role of more sophisticated features, such as hierarchical patterns and metric correlations, are necessary for a correct modeling of these neural systems.