Emergence of Log-Normal Type Distributions in Avalanche Processes in Living Systems: A Network Model

Actin is the major cytoskeletal protein of mammal cells that forms microfilaments organized into higher-order structures by a dynamic assembly-disassembly mechanism with cross-linkers. These networks provide the cells with mechanical support, and allow cells to change their shape, migrate, divide and develop a mechanical communication with their environment. The quick adaptation of these networks upon stretch or compression is important for cell survival in real situations. Using atomic force microscopy to poke living cells with sharp tips, we revealed that they respond to a local and quick shear through a cascade of random and abrupt ruptures of their cytoskeleton, suggesting that they behave as a quasi-rigid random network of intertwined filaments. Surprisingly, the distribution of the strength and the size of these rupture events did not follow power-law statistics but log-normal statistics, suggesting that the mechanics of living cells would not fit into self-organized critical systems. We propose a random Gilbert network to model a cell cytoskeleton, identifying the network nodes as the actin filaments, and its links as the actin cross-linkers. We study mainly two versions of avalanches. First, we do not include the fractional visco-elasticity of living cells, assuming that the ruptures are instantaneous, and we observe three avalanche regimes, 1) a regime where avalanches are rapidly interrupted, and their size follows a distribution decaying faster than a power-law; 2) an explosive regime with avalanches of large size where the whole network is damaged and 3) an intermediate regime where the avalanche distribution goes from a power-law, at the critical point, to a distribution containing both 1) and (ii). Then, we introduce a time varying breaking probability, to include the fractional visco-elasticity of living cells, and recover an approximated log-normal distribution of avalanche sizes, similar to those observed in experiments. Our simulations show that the log-normal statistics requires two simple ingredients: a random network without characteristic length scale, and a breaking rule capturing the broadly observed visco-elasticity of living cells. This work paves the way for future applications to large populations of non-linear individual elements (brain, heart, epidemics, … ) where similar log-normal statistics have also been observed.


INTRODUCTION
Avalanche processes are very common in living systems, such as firing rates in brain [1], fractures in living cells cytoskeleton (CSK) [2], but also in amorphous [3] and random media [4], or earthquakes. Actually, all the systems that can be considered as composed by elementary threshold units, which are connected and can transfer information to each other, are subject to avalanches. Indeed, as soon as a unit reaches its threshold it can cause other units to do the same in turn.
In statistical physics, avalanches are often treated in the framework of critical systems. This is justified by the observation of critical behavior of avalanche statistics, with distributions usually approximated by power-laws, at least for some length or energy scales. Power-laws are reminiscent of self-organized criticality [5], and are ubiquitous for avalanches in solids and amorphous materials [6][7][8][9], describing for example the avalanche sizes of a granular material falling apart. Self-organized criticality has also been associated to neural dynamics (first in [10] and later in [11,12]), showing widely spread power-law statistics with a typical exponent invariant among different species [12] and in different environmental situations. Lately, sleep-wake micro-architecture and regulation have been identified as nonequilibrium processes to maintain a critical state of the brain network [13,14]. In all these examples, fluctuations, and thus their distribution, play a major role in triggering the system to the critical state. Recently, some doubts about the ubiquity of power-law distributions have been advanced [15], pointing out that many real system's data are actually better fitted by other skewed, fat-tailed distributions, such as the log-normal. This is true for network degree distributions, in social, biological or technological networks [15], but also for avalanche statistics, in neural networks [16,17] or living cells cytoskeleton [2]. Indeed, the growth of a network can be seen itself as an avalanche process, thinking for example of the preferential attachment rule [18]. In general, power-law distributions are difficult to irrefutably hold for any kind of real data. Remarkably, power-law and log-normal distributions are closely related to each other and small variations in the generative mechanisms can lead to one or the other [19,20].
Nevertheless the power-law, with some adaptations, has been the largely dominant model in the last decades for all data showing a fat-tailed distribution, either for network theory or for avalanche processes. There are mainly two reasons for the large predominance of the power-law modeling up to now: 1) because much is known about the modeling of power-laws, thanks to the statistical physics of phase transitions, leading to an ease of data treatment and interpretation. Concerning network science theory, the same is true for scale free networks [21], justifying a power-law interpretation of real data distributions, while models for other fat tail distributions do not exist so far [15,22]; 2) due to the difficulty to fit real data with fat-tailed distributions, because a power-law tail is only visible beyond a given threshold of the random variable, and not in the whole domain. Actually, any skewed distribution can be approximated by a power-law, if only a finite and small scale range is taken into account (see Figure 1C). In Figure 1, comparing the distributions of hippocampal firing rates in a linear scale ( Figure 1A), in a log-lin scale ( Figure 1B) and in a log-log scale ( Figure 1C), we realize that focusing on the tail of the distribution in a log-log representation increases the risk of missing some essential features of the underlying process [17].
The point 2) led to an ongoing debate about what is the most appropriate fitting distribution for the considered data. Some concepts such as low-degree saturation and large-degree cutoff have been developed for example in network-science theory, in order to account for important characteristics of real systems, such as their finite size [22]. This solves point 2), adapting the power-law distribution to be applied to real systems. On the other hand, point 1) reveals the need of an alternative modeling framework which could validate the choice of other fat-tail distributions, in order to understand the underlying mechanisms leading to one distribution or the other. In the same way as degree distributions of networks, so far there are no models accounting for log-normal distributions of avalanches on networks, while there is strong evidence that log-normally distributed avalanches exist. In Figure 1, we show a few examples of processes related to an avalanche dynamics showing log-normal distributions: the distribution of the spine (small dendritic protrusion crucial in the transmission of electrical signals) sizes in neural networks [23] ( Figure 1D); the distribution of new numbers of Ebola cases per week [24] ( Figure 1E); the distribution of lags between two consecutive vehicles to cross the stop line at a crossroad [25] ( Figure 1F).
Beyond power-laws and log-normals, other skewed distributions are sometimes used to model biological data. For example, it is worth mentioning the Gamma distribution used to model inter-spike intervals of neurons [26] or the inter-beat interval variation of heartbeats [27,28]. We should notice that statistically it is not always easy to distinguish log-normal from gamma distributions. In any case both are alternative to the power-law framework.
Originally motivated by the experimental observation of log-normal statistics in avalanches of fractures of the CSK of living cells [2,29], we focus in this paper on the modeling of log-normal avalanches on random networks. Our model is inspired by our previous 1-D model [2] and by works on epidemics spreading on networks [30,31], in our case the population is represented by actin filaments, being in two possible states: cross-linked or not. The network structure models the CSK structure and the avalanche is a model for rupture mechanics in living cells, when plastic deformations are considered. Given that we do not consider here active cross-linking carried out from myosin filaments, the same framework can be applied to other cross-linked polymers with glassy dynamics. Therefore similar conclusions can be applied to amorphous glassy materials, like polymers, metallic glasses or colloidal glasses, which all share slow dynamics, and long mechanical relaxation delays [32]. This can be useful for understanding all processes not showing good power-law statistics, and in general what makes a distribution shifting from power-law to log-normal, highlighting the characteristics that a process needs to have to deviate from the most common power-law modeling.

Brief Introduction on Cytoskeleton Mechanics
Among all the fascinating properties of living cells, we must emphasize their ability to constantly remodel their structural organization to withstand forces and deformations and to promptly adapt to their mechanical environment [33,34]. This versatility is fundamentally required for many vital cellular functions, such as migration, mitosis, apoptosis or wound healing, and an alteration of the cell mechanical properties can participate in pathogenesis and disease progression, such as cancer [35,36].
We focus here with more details on the actin cytoskeleton. These filaments are the most relevant for our modeling, since they cover the perinuclear zone of cells, which is the cell compartment poked by our micro-indenting tips. They present a polar structure [43], a quite fast, with respect to global active reorganization processes, polymerization rate (larger than one per minute) and a depolymerization rate a FIGURE 1 | (A) Non-normalized distribution of firing rates of hippocampal CA1 pyramidal neurons during slow-wave sleep. The firing rate is measured in Hz; (B) Distribution (non-normalized) of the logarithm of the firing rate as in (A); (C) Distribution of the firing rate as in (A) in a log-log plot. Notice that by taking into account only the tail of the distribution (black dots), the distribution could be fitted by a power-law, even though we may be missing some important information; (D) Distribution of the spine size in arbitrary units (AU) (top), and probability density of the logarithm of the spine size (bottom); (E) Distribution of the new number of Ebola cases per week (week incidence); (F) Probability density of the lag times at which 2 vehicles waiting at rest at a traffic light pass the stop line. This is what is called headway. The distributions are reported for the lag between two successive vehicles considering the second position (top) and the fourth position (bottom) of the line: for instance the bottom plot describes the distribution of the difference of times at which the third and the fourth vehicle of the queue cross the stop line. All times are given in seconds; (A) (B) and (C) are adapted from [17] (D) is adapted with permission from [23] (E) is reproduced with permission from [24] and (F) from [25].
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org January 2021 | Volume 6 | Article 613962 little slower, of the order of one per a few minutes [44]. Moreover they can build a cross-linked network whose mechanical properties depend in general on the cross-linker proteins density and on the network structure. Actin filaments networks are designed by a wide variety of actin-binding protein cross-linkers, which can be passive or active, i.e., activated by ATP (Adenosine Triphosphate) hydrolysis, the latter having a slower dynamics (time scale of tens of minutes) [37-42, 45, 46]. This cross-linked network gives to living cells some properties of soft-glassy materials [47,48], such as the weak power-law dependence of the shear relaxation modulus G with both time and frequency. This behavior was first modeled by an empirical law known as structural damping from material engineering [49] and later on associated to the fractional visco-elastic Kelvin-Voigt model (a springpot and a dashpot in parallel), see e.g., [50]. The power-law decay (in time) is quite impressively not depending much on the particular experimental technique, neither on the different type of cell: the decay exponent α mostly belongs to the interval [0.2-0.4] [47,[51][52][53][54]. Notice that α 0 corresponds to a purely elastic response and α 1 to a purely Newtonian fluid response. Other models for cell rheology have been used (see for instance [50,55]), but they all require a fractional visco-elasticity to capture the cell response, resulting in a power-law decay G(t) ∼ (t/τ) − α or in more complicated functions such as the Mittag-Leffler function [56], which can be approximated for t/τ → 0 + by the stretched exponential [57]: For our purposes the three functional forms will be considered as equivalent, since they all account for a slow relaxation dynamics given by the glassy structure and thus a memory of the past deformation.
Considering cells as soft glassy materials, we can extrapolate that they are constructed from a disordered structure of connected discrete elements by weak attractive interactions. Each of these elements would be in a metastable state [48], allowing cells to flow, and therefore prone, for instance by an external forcing, to generate avalanches of fractures, which are typically out of equilibrium processes. Indeed, by tuning the proportions of passive or active actin-binding proteins, and reorganizing the network structure, living cells can control their power-law (scalefree) CSK rheology [37,58]. Interestingly, cells exhibit both solid and liquid-like properties. Solid-like behavior is associated with strongly cross-linked actin filaments which resist sliding and accumulate tension [44,59], while weakly cross-linking proteins produce actin filaments which slide more readily, enabling the network to flow as a liquid [60].
This paradox can be solved within the theory of soft glassy materials, by considering that, upon external deformations, the CSK of a living cell can undergo deep structural transformations such as the unfolding of protein domains, the unbinding of cytoskeletal cross-linkers, and the breaking of weak sacrificial bonds. All these structural changes are inelastic (non-reversible in a strict sense), they dissipate locally the elastic energy of the CSK network (structural damping) [61,62].
The ability of cells to switch quickly from fluid-like responses to more brittle solid-like responses [2] is directly linked to the interplay between stability/rigidity and flexibility, in a way similar to what happens in neural networks [63]. This fast switch is likely driven by avalanche processes, which allow fast transfer of information. At the same time, such events reduce the connectivity of the CSK and may result in permanent plastic deformations or even more dramatic irreversible failures [38,39] which, for instance, could be at the origin of the recently observed incomplete shape recovery of living cells after repeated creep [64]. These effects are reminiscent of those in cyclically loaded solids which can lead to fatigue-induced failure [4,6,7].
The observation of universal relationships governing cell (and not only) rheology is evocative of the universality of statistical mechanics systems, such as the Ising model [65], in which individual details of the filaments and particular molecular interactions are unimportant to identify global behaviors.

Poking Living Cells With a Sharp Atomic Force Microscopy Tip
Previous experiments done in our group motivated this study [2,66,67]. We introduce here the main characteristics of these experiments useful for the modeling. Atomic force microscopy (AFM) was operated in force-spectroscopy mode, therefore giving as result the force indentation curves, i.e., the plot of the instantaneous force against the indentation length inside isolated adherent cells. The indentation was stopped at a certain set-point force, after which the motion of the AFM tip was inverted. This was achieved after some calibration steps whose details can be found in [2,66,67]. It was important to carry out this indentation with very sharp tips (pyramidal or conical) with a tip curvature radius of only a few nanometers, allowing their penetration inside the meshes of the cross-linked network. The indentation was performed at constant velocity (1 μm/s), and therefore constant strain rate, so that one indentation-retract experiment lasted for only a few seconds.
The cell compartment poked by the AFM tip was the region just above the nucleus (perinuclear), which besides being better recognizable, is also very rich of actin stress fibers and cross-linked filaments. As previously noted [2,66,67] this indentation caused locally a strain stiffening, signature of an increased tension in the network, eventually resulting in local singularities in the force indentation curves, interpreted as avalanche of fractures. We were able, thank to the wavelet transform mathematical microscope [68][69][70], to quantify these singular events and characterize them, by their force drop, indentation length and finally energy released [2].

Rupture Event Statistics From Two Primary Cell Lines
Some examples of the statistics of these singular fracture events, together with global shear relaxation moduli distributions are shown in Figure 2. We can see that both the force drop F d (Figure 2A), caused by the avalanche of failures, and the indentation length ΔZ ( Figure 2C) over which the avalanche Frontiers in Applied Mathematics and Statistics | www.frontiersin.org January 2021 | Volume 6 | Article 613962 takes place are log-normally distributed (the plots show the distribution of the logarithm of the variable, then a log-normal in a log-log representation is a parabola). For the indentation length ΔZ we could separate two different regimes of avalanches both of which approximately log-normally distributed. This lognormal statistics was found with different cell lines: in Figure 2A are shown myoblasts (red), myotubes (blue) and ATP depleted myoblasts (black), while in Figure 2C are shown immature CD34 + hematopoietic cells (blood cells) healthy (blue) and leukemic (red). It follows that the energy released E d F d ΔZ, not shown here, is also log-normally distributed. The same lognormal distribution is also observed on global quantities such as the global shear relaxation modulus G g (see Figure 2B,D), extracted by the identification of the whole force indentation curve with the Sneddon model [71]. We speculate that the log-normal statistics observed in different living cell types on macroscopic elastic moduli (see also results on breast tissue cells from [55]) is strictly connected to the microscopic processes taking place in the polymer network, i.e., avalanches of reorganization fractures. Indeed, the log-normal statistics of avalanches would reflect in a log-normal noise (skewed fluctuations are observed also e.g., in [55]) biasing the estimation of the global elastic moduli. We can thus interpret the shear-relaxation modulus as the response of a sum of microscopic avalanches, as results of a dynamical reorganization of the network structure. This interpretation is supported by Sollich's theory of soft glassy rheology [48], in which the glassy polymer is interpreted as composed of individual units in metastable energetic states (with different energy depths) and in which rearrangements are due to disordered interactions summarized by an effective temperature.

A Random Network Model for the Cell Cytoskeleton
Let us now introduce our random network model for the cell CSK. We consider a network with N nodes, being identified as the actin filaments, connected by randomly assigned (in a way explained hereafter) links, identified as the cross-linker proteins. In light of what we said previously about CSK mechanics (see section 2.1), if we want to build a random network model for the cytoskeleton we need the network to be in a percolating regime. Indeed, since a cell can be seen as a soft glassy material there needs to be a giant, percolating cluster of connected nodes, allowing for the transition between a fluid material and a glassy one. A second important condition that we ask the network to satisfy is the correct degree distribution (i.e., the distribution of cross-links per filament) as observed in real living cells CSK. Unfortunately, precise data on this distribution do not seem to be available to our knowledge. However, from the literature [72], we can deduce that the degree distribution would not be well approximated by a power-law and therefore the CSK could not be modeled as a scale-free random network. There is also another reason why the scale-free network may not be a good model for cell CSK: it has been proved [73] that on scale-free networks the critical threshold for an epidemic is exactly 0, meaning that avalanches would always affect a finite proportion of the population (in the large N limit), thus implying irreversibly a large proportion of crosslinkers, scaling as the size of the system. This is not what is observed for avalanches of fractures in cells, since their distributions are of finite size. For all these reasons, we propose a Gilbert [74] random graph to model the CSK structure. A Gilbert graph is a variant of the Erd} os-Rényi graph [75], and the two graphs are equivalent in the large N limit. The network is constructed in the following way: 1) we define a network of N isolated nodes, with N fixed; 2) we connect each possible pair of nodes with probability p l .
This process creates a network with a binomially distributed degree of connections: where the binomial factor takes into account all the possible pairs of nodes having a degree of k, out of all possible N − 1 links coming out from a node. The probability p k indicates the probability that a randomly picked node has k links. The average degree of the random network can then be computed: In the limit of N → ∞, keeping 〈k〉 fixed, the binomial degree distribution is well approximated by the simpler (one parameter) Poisson distribution, noted P(λ) k , with parameter λ 〈k〉. The average connectivity 〈k〉 will be a control parameter of our model, since its value is not well-known from experiments. The average global number of links of the network, and thus of cross-linkers available for the avalanche, can also be computed by multiplying the probability of connection of pair of nodes p l with the number of all the possible pairs (by excluding double counting): and then for a network of N 10 4 nodes, and p l 5 · 10 − 4 , the average number of links is 〈N l 〉x50000. The distribution of the number of links is also a binomial. Before moving on and describing the avalanche process, we would like to observe some characteristics of this network useful for our modeling. From probabilistic arguments it can be noticed (see [75]) that, in the limit of N ≫ 〈k〉, with fixed 〈k〉 (as in all the simulations of this article), the fraction u of nodes not belonging to the largest connected component of the network is given by the transcendental equation: This can be found by considering that the fraction of nodes u is composed by the sum of two separated events. A node i not belonging to the giant component is: i) either disconnected to an other node j because the link i − j does not exist, and this happens with probability 1 − p l , ii) or connected to a node j that also does not belong to the giant component, and this occurs with probability up l . Therefore u (1 − p l + up l ) N− 1 , because j can represent N − 1 different nodes. From Eq. 5 it can be proved [22,75] that for 〈k〉 > 1 a giant percolating cluster exists and its size scales as the network size N. At this point we can work out the probability p FC that we pick up a random link not belonging to the giant cluster: where the right hand side is the ratio of the average number of links not belonging to the giant cluster p l Nu/2(Nu − 1) and the average number of links of the network p l N/2(N − 1). On this network structure, an avalanche of ruptures is described by the following algorithm. Notice that while the probability of connection p l remains constant, the probability of breaking p(t) can depend on time, and therefore is able to capture time relaxing processes.
• we pick at random a link of the network and we break it (removing it from the network), this initiates the avalanche • we look for the links coming out from both extremities of the broken link; • we break all of them with a probability p(t), where at the first time step t 0; • we consider the links coming out from the ones broken at the previous time step from both extremities and break them with probability p(t + 1).
We repeat the last point and increase the time step by one unit, until the avalanche stops when there are no more links that break, either because the probability of breaking is too low or because of the damage caused to the network, which will have decreased the abundance of unbroken links. The time at which the avalanche stops, t * , defines the duration of the avalanche and the total size of the avalanche Z is the total number of broken links. Here, we focus on the distribution of Z, which is, up to some conversion factors, equivalent to the energy released during the avalanche. Overall the stochastic process has two sources of randomness: the random structure of the network and the avalanche process itself which is a stochastic process.
Rupture avalanches propagate on networks that remain static during the rupture process, since actin polymerization and active cross-linkers act on time scales much longer than the time scale of the experiment, which is of the order of 1 s. Nevertheless this assumption is relaxed in section 3.4, to account for the fact that cross-linkers can be restored on very short time scales even by physical contact.
Statistically speaking an avalanche is a binomial process of probability p(t) of success, with a random number of trials (the number of available links). If the number of available links were always drawn from the same Poisson distribution, with parameter λ 2Np l (the factor 2 is because we consider both extremities of a broken link), we could conclude that the distribution of broken links at time t is also a Poisson distribution with a new parameter λ 1 2p(t)Np l . The choice of the Poisson distribution is motivated by the degree distribution, which, as we said, is approximately Poissonian in the N ≫ 〈k〉 limit keeping 〈k〉 fixed, and the parameter λ 2Np l corresponds to the average number of links available at the first time step. This can be shown by applying the law of total probability to the probability of having k broken links: The probability of having k broken links is then given by the product of the Poisson distribution, giving the probability of having n links available for breaking, times the binomial distribution giving the probability of having k successes out of n trial. To include all possible disjoint events we have to sum over n > k, up to N − 1, which in the considered limit tends to infinity. The problem with this reasoning is that the number of available links is not drawn from the same distribution as far as the avalanche moves forward. Indeed the number of available links depends on the particular path of the avalanche, since the broken links disappear from the network. Therefore as far as the avalanche progresses, the number of available links may not even follow a Poisson distribution, making the analytic expression of the number of broken links impossible. The same is true if we consider the restoration of crosslinks, as in section 3.4, because the avalanche process introduces an effective and uncontrolled time dependence of p on t, by changing the actual number of available links. This effect is the main reason why we do not find perfect log-normals in our model, as in our first 1D model [2] or related models without a network structure [20]. Actually, this was checked on a random multiplicative process with a multiplicative factor drawn at each time step from the same distribution as in Eq. 7 for different parameters λ 1 proving that, as expected from [20], there is a region where the avalanche size distribution is exactly log-normal.
We considered 〈k〉 > 1 in all our simulations, in order to have a giant component in the network, as explained previously. We also checked the effect of picking the first link initiating the avalanche out of the giant component. Indeed, from Eqs 5, 6 we can compute p FC . For 〈k〉 2 for instance p FC 0.04 is already irrelevant for our results, and for 〈k〉 4, p FC x4 · 10 − 4 , is completely negligible. This was checked numerically in all the following simulations and the results do not change by running the avalanches only on the giant component. We should also mention that avalanches containing only the first randomly chosen link are not taken into account for the statistics, this making even more negligible the effect stated above.
It is conceptually interesting to note that the proposed avalanche model is actually equivalent to an epidemic model running with probability p(t) on the line graph of the original network, i.e., on the graph obtained by transposing links into nodes and nodes into links.
Simulations of avalanches were done with self-written codes using open source software R version 4.0.2 and Python 3.6. Figures and data analysis were done with Python 3.6 using basic packages and a customized version of the module powerlaw [76] with a corrected computation of the cumulative distribution function and added graphical features. For all the following simulations it was checked that changing the size of the network N ((but keeping the same values of 〈k〉), the phase diagram did not change significantly, showing always the same types of distribution for N 1000, 5000, 20000, 50000. It was also checked that running the avalanche algorithm on the same network for each stochastic realization (of course with all links restored), was equivalent, since the first initiating link is chosen randomly, to run the avalanche on a different network for each realization. The first option is computationally faster and was then used for detailed results.

Avalanche Statistics With Constant Probability of Breaking
First we implemented the avalanche model with a probability of breaking p(t) constant with time. We detected three regimes with different distributions of the avalanche size. A phase diagram is shown in Figure 3B. The three phases are pop (blue), mixed (purple), explosive (red). The pop regime presents avalanche size distributions which decrease faster than a power-law, containing then small avalanches, with a low number of broken links of less than 100. On the other hand the explosive regime is composed of large avalanches, scaling as the system size and then causing global damage to the network. In between these two regimes the mixed regime is composed of avalanches distributions made of a mixture of the two previous regimes.
A likelihood-ratio test, based on the ratio of the likelihood of the model fitted by two different distributions, can be done to compare possible outcome distributions [77]. In Figure 3D, we show the maximum of likelihood (ML) fit for the log-normal distribution and the truncated power-law distribution, with density function kx −α e −λx (a power law with an exponential cutoff with λ decay rate and k normalization constant).
When applying the likelihood-ratio test to compare the two distributions it turns out that the log-normal is always a better fit than the truncated power-law. For the example of Figure 3D, the p-value is 1 · 10 − 22 , rejecting the null hypothesis that the likelihood ratio is equal to 1, and indicating that none of the selected distributions is preferable. This result is coherent with those obtained on 1-D models in [20], where for low multiplicative factors, the size distribution is approximately log-normal, with the difference that here the process is discrete (not allowing sizes < 1). Therefore we have only the right tail of the distribution, while in the 1-D models, the multiplicative process is continuous. We can thus conclude that in the pop regime the distribution is statistically compatible with a log-normal tail.
The boundary between the pop and mixed regimes is computed by an algorithm detecting for which parameters 〈k〉 and p the avalanche size distribution is better fitted by a truncated power-law over the whole domain. We chose this criterion because power-law distributions are typical of critical points and we wanted to include the finite size cutoff (see for instance Ch. 4 of [23]). Moreover, when comparing the tails of the distributions (by using the algorithm described in [78] and implemented in [76]), at the boundary the most likely distribution is the truncated power-law, with always a very significant p-value ≪ 10 − 4 . The null hypothesis is again of a likelihood ratio equal to 1 (between log-normal, exponential, power-law). In Figure 3C is shown a typical distribution of the avalanche size at the boundary between pop and mixed. For comparison we also show the maximum of likelihood fit with a power-law.
In the mixed regime the size distribution consists of two separated parts: a small size part similar to the pop regime distribution and a large size part corresponding to a very narrow bump. In Figure 3E we show the full Z distribution in a log-log plot, in the inset is shown in detail the distribution of the large avalanche size narrow bump, represented by only one dot in the large field view. In this regime, moving toward the explosive regime causes a shift to the right of the large avalanche size bump and a decrease of importance of the small avalanche size part. Notice that the large avalanche size bump is not disappearing in the large N limit: if we take a distribution in the mixed regime for a given p and k and we increase the size of the network N (we did it for N 2000, 5000, 10000, 20000) the proportion of avalanches in the bump stays the same, even though the bump shifts to the right, because those avalanche sizes are proportional to the system size.
Finally the explosive regime, of which a typical distribution is shown in Figure 3F is detected by the criterion of less than 1% of small size avalanches. In this regime the distribution of avalanche sizes is a very narrow bump having a size comparable to the number of links of the system. The integrity of the whole network is therefore compromised and only a few sparse links are intact.
Note that even though it could be tempting to say that the bump in mixed and explosive regimes is a log-normal distribution (remember that in a log-log representation a log-normal becomes a parabola), this cannot be assessed, because the distribution is very narrow and therefore a lognormal cannot be distinguished from a normal distribution. In none of the three regimes observed so far we can recognize avalanche size distributions similar to those observed experimentally, i.e., approximately log-normal and of finite size. Nonetheless, we observe that the distributions of the pop regime (showing a log-normal tail) are strikingly similar to those observed for avalanches of firing rates in freely behaving rats in [16]. In that work, even though the firing rate distributions were also claimed to follow a lognormal tail, it was justified from data undersampling, keeping the hypothesis of a critical state (and therefore power-law distributed). Our simulation results suggest another possible reason of the observed distributions, as given either by a small connectivity of the network or by an extremely fast absorption of the avalanches stress here the fact that our network model does not have a metric, hence no concept of distance, it is only a network of interactions. It does not matter if the links are close in space (see Figure 3A), but only if they format of <<>>. However, it is possible to embed the random network model in any metric space and in any shape, considering for example nodes that are connected to each other closer than others, without impacting our results: as far as the interactions follow this structure the shape of the network does not matter. With this interpretation the small pop avalanches in the mixed regime can be thought as a boundary effect: they represent avalanches where the starting randomly chosen link is in a less than average connected region.

Avalanche Statistics by Introducing Visco-elasticity
We now introduce the exponential relaxation coming from fractional visco-elasticity (see section 2.1). We let then the probability of breaking be dependent on time in the following way: The Γ coefficient in Eq. 1, is here embedded in the constant τ. This relation imposes some temporal correlation during the avalanche, added to the interaction structure. Conceptually, Eq. 8 assumes that the local fracture mechanics follows the same mechanical response as the global network given by the shear-relaxation modulus (see section 2.1). It is important to note that the limit τ ≫ 1 is equivalent to setting α 0, since in both cases p(t) loses the dependence on time and we recover the results of section 3.2. We can thus interpret the model in section 3.2 as a purely elastic response.
In Figure 4 we show the results after introducing a time dependent p(t). The phase diagram also shows three different regimes. The pop and the explosive regimes have the same statistics as in section 3.2, and an example of them is shown in Figure 4C,D.
The important difference introduced by a fractional viscoelastic relaxation is in the mixed regime, which shows now the emergence of an approximately log-normal finite size bump, and it is thus called the log-normal bump regime, in magenta. Indeed the very narrow bump observed in the non-visco-elastic model broadens here due to the time relaxing probability and it is not at a size comparable to the size of the system, therefore not causing a global destruction of the network. Observing Figure 4B we can see that there are still some small avalanches of pop type in the distribution. These may still be an effect of picking low reticulated zones as starting point of the avalanche, but as previously said, even running the avalanches only on the giant cluster they are still present. Notice that while moving toward the red region the lognormal bumps shrinks and moves to larger sizes becoming then similar to the distribution in Figure 3E.
The boundaries between the regimes are detected with the same algorithms as described in section 3.2. Notice that here we can let the probability p 0 become larger than 1, since the stretched exponential decay will always make p(t) converge toward 0. Simulations in Figure 4 are shown for τ 4, and the influence of varying t has also been studied. Increasing τ produces a shift of the distribution to the right, since the avalanche tends to last longer. However for having a log-normal bump in the size distribution, τ must not be too large, because otherwise the distribution would tend to that one shown in Figure 3, with a larger proportion of small avalanches similar to the pop regime and a very narrow large avalanche size bump. On the other hand, decreasing τ would make the relaxation very fast and shift the distribution toward the pop region. If decreasing τ is compensated by an increase of 〈k〉 the weight of the small size avalanches in 4B is decreased, making the log-normal bump even more pronounced (see Figure 5), and broadening the magenta region. As a rule of thumb, τ < 5 (in arbitrary units) is a good compromise to find an approximately log-normal bump. The value of α, coherently to rheological experiments on living cells, is set to 0.3 in all the simulations. Therefore the best receipt to have log-normal without having small pop avalanches is then to decrease τ in such a way that if increasing p 0 or 〈k〉 the avalanche stops fast enough to escape from the red regime (and then a very narrow large avalanche size bump) possibly together with a few pop avalanches, but still p 0 or 〈k〉 has to be large enough to quit the blue region. In other words we get lognormal avalanches if the avalanche starts explosive and relaxes fast, ending gently. We showed this in Figure 5, where the same simulations as before were run, but with τ 1. Comparing Figures 4, 5 we can observe that the magenta region is broadened, and in the new area, which from red became magenta, distributions with a very low proportion of pop avalanches emerge.
We note that these new distributions are not exactly lognormal, but only approximately. It is for example possible that they could be as well approximated by other skewed distributions, as the Gamma distribution, but for sure not by power-laws. It is interesting that these skewed distributions of the log-normal bump region are originated by both the introduction of time correlation (in the form of visco-elastic relaxation) and the nonlinear dynamics of the avalanche propagation. Similar arguments (but instead of time correlations it was for phase correlations) were used in [27,28] to explain the appearance of skewed distributions for the temporal variability of the heart rate, suggesting that an avalanche process may also be latent in that case.

Avalanche Statistics With Local Restoring of Cross-Linkers
In order to test the effects of local cross-linkers restoring on the final avalanche size distribution, we consider the most extreme case, where all broken links are restored after only one time step. Therefore the broken links disappear from the network only for the first time step after their rupture. This could be interpreted as an extremely fast network restoring. In this scenario, the probability of breaking must relax with time, otherwise the avalanche would never stop as soon as λ 1 > 1 (see Eq. 7). We thus take the same law of breaking as in Eq. 8, with τ 4. The phase diagram and the resulting avalanche size distributions, when propagating avalanches with local restoring of cross-linkers, are almost the same as in section 3.3 (compare Figures 4, 6). We have chosen exactly the same points in the phase diagram as in Figure 4 for ease of comparison. The pop region is very similar to the one in Figure 4. The two other distributions have the same shape as before, but they are shifted to larger values, because the number of links available for the avalanches is more abundant than in Figure 4, this is a direct consequence of the link repair mechanism introduced for these simulations. We can thus conclude that our results about the lognormal bump regime are robust and do not seem to depend on the details of the particular avalanche algorithm, provided a viscoelastic relaxation is included in the model. However these details can be important to tune the average and the width of the resulting distributions.

DISCUSSION AND FUTURE DIRECTIONS
Our phenomenological model shows that it is possible to model approximately log-normal avalanches, as experimentally observed. We showed this on a network structure which does not have a natural metric, but it is primarily a structure of interactions. Therefore interactions are not necessarily due to physical contacts, but can for example be given by vibrational modes, or other forms of interaction. However, this characteristic should not be thought of as a limitation, but rather as a feature that makes the model more general. Indeed, the network structure can always be embedded in some metric space if desired, e.g., on a 3-D sphere, and the same results would remain valid.
Interestingly, in order to have approximately log-normal avalanches, an elastic-instantaneous breaking mechanism is not sufficient. With this first rule we only get small pop avalanches, explosive avalanches leading to a global destruction of the network, or a mixture of both regimes. Along the critical line, where large avalanches proportional to the network size start being present, the distribution looks like an exponentially truncated power-law, typical of critical systems in physics. A key ingredient for having log-normal type avalanches is the introduction of a breaking rule taking into account the viscoelastic memory of cells (or other glassy materials), thus introducing time correlations. This behavior is shown to be robust even with different avalanche propagation rules, as for example by introducing the local restoration of cross-linkers. We should also mention that even with a completely different breaking rule the qualitative avalanche behavior is the same. We did preliminary studies on avalanches with a propagation rule consisting in breaking at each time step only the weakest link (by setting randomly attributed <<strength>> to the links), supposing that the avalanche propagates along the weaknesses of the material, and we could again recognize the same three regimes. We were then able to obtain avalanche size distributions closer to a log-normal by imposing faster relaxation times as in Figure 5. In this way the log-normal bump region becomes larger, thus showing slightly new distributions with a less important proportion of small pop avalanches, while the other resulting distributions stay the same.
Our model does not give quantitative information, because we did not make the link with physical units. Accounting for this, however, may be possible by introducing an appropriate scaling factor to make the size of the avalanche a physical measurable quantity (as the released energy). At the same time the correspondence of the lasting time of the avalanche with the experimental one can be done. This would attribute a unit to τ, hence leading to an estimation of the physical relaxation time needed to have log-normal avalanches. Altogether, these results suggest that, for cells, the local fracture mechanism resembles that of the global mechanical response of the cytoskeleton network. The global shear-relaxation modulus can thus be interpreted as a sequence of local avalanches of ruptures.
We conclude by saying that, in contrast to simpler 1-D models without a network structure (see [2,20]), here the resulting distributions are not exactly log-normal. It is for example possible that they can as well be approximated by skewed distributions other than the log-normal, as the Gamma distribution, yet our model is an alternative to the power-law modeling. Furthermore our avalanches follow a skewed distribution which does not depend on the system size. Notably, these skewed distributions in the log-normal bump region are originated by both the introduction of time correlations (in the form of visco-elastic relaxation) and non-linear dynamics in the avalanche propagation. These are the same ingredients (but instead of time correlations it was for phase correlations) identified to cause the appearance of skewed distributions for the variations in the heart-rate signal in [27,28], behind which an avalanche process may also be latent. This paves the way to possible applications to physiological data, where non-linear units are organized in a networked communicating structure, such as brain seizures, heartbeats or chemical signaling inside cells.