Skip to main content


Front. Appl. Math. Stat., 28 January 2021
Sec. Dynamical Systems
Volume 6 - 2020 |

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

  • 1University Bordeaux, CNRS, LOMA, UMR5798, Talence, France
  • 2Institute for Complex Systems and Mathematical Biology, SUPA, University of Aberdeen, Aberdeen, United Kingdom

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.

1 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 [69], 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 non-equilibrium 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].


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].

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.

2 Evidence for Log-Normal Statistics in Cell Mechanics

2.1 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].

The mechanical properties of living cells are mediated by their CSK, a dynamic network of filamentous proteins composed of actin filaments, microtubules, and intermediate filaments [3742].

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 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) [3742, 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, 5154]. 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 (scale-free) 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.

2.2 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 [6870], to quantify these singular events and characterize them, by their force drop, indentation length and finally energy released [2].

2.3 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 Fd (Figure 2A), caused by the avalanche of failures, and the indentation length ΔZ (Figure 2C) over which the avalanche 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 log-normal 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 Ed=FdΔZ, not shown here, is also log-normally distributed. The same log-normal distribution is also observed on global quantities such as the global shear relaxation modulus Gg (see Figure 2B,D), extracted by the identification of the whole force indentation curve with the Sneddon model [71].


FIGURE 2. (A) Probability distributions of the logarithm of the local force drops Fd(nN) estimated from local disruption events collected from sets of myoblasts (red), myotubes (blue) and ATP depleted myoblasts (black); (B) Probability distributions of the logarithm of shear relaxation modulus Gg(kPa) estimated by a parabolic fit of the Force-Indentation Curves (Sneddon model [71]), from healthy (blue) and leukemic (red) immature CD34+ hematopoietic cells; (C) Distribution of the indentation lengths of the rupture events detected from healthy (blue) and leukemic (red) immature CD34+ hematopoietic cells; (D) Probability distributions of the logarithm of shear relaxation modulus Gg(kPa) on the same sets of cells as in (A). Note that plots (A) and (C) are in log-log scale and plots (B) and (D) are in semi-log scale. All the logarithms are here in base 10. Panels (A) and (D) are adapted from [67] (B) from [2] and (C) from [2].

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.

3 Random Network Presentation

3.1 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 cross-linkers, 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ős-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 pl.

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 N1 links coming out from a node. The probability pk 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 pl with the number of all the possible pairs (by excluding double counting):


and then for a network of N=104 nodes, and pl=5104, the average number of links is Nl50000. 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 Nk, 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 ij does not exist, and this happens with probability 1pl, ii) or connected to a node j that also does not belong to the giant component, and this occurs with probability upl. Therefore u=(1pl+upl)N1, because j can represent N1 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 pFC 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 plNu/2(Nu1) and the average number of links of the network plN/2(N1).

On this network structure, an avalanche of ruptures is described by the following algorithm. Notice that while the probability of connection pl 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

• 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 λ=2Npl (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)Npl. The choice of the Poisson distribution is motivated by the degree distribution, which, as we said, is approximately Poissonian in the Nk limit keeping k fixed, and the parameter λ=2Npl 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 N1, 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 pFC. For k=2 for instance pFC=0.04 is already irrelevant for our results, and for k=4, pFC4104, 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.

3.2 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.


FIGURE 3. (A) Picture of the random network model, embedded in a sphere. Notice that this is just a choice of representation, the shape of the network is not fixed: our network model is only a set of connections, it can be embedded in any metric space. Dots are the nodes of the network and lines are the links. For better visibility the network has here only N=1000 nodes and a probability of connection pl=0.003; (B) Phase diagram of the possible resulting avalanche size distributions with respect to the model parameters k (the average number of connections of the network) and p (the probability of breaking a link). The three identified regimes are: pop (blue), which has an avalanche distribution statistically compatible with a log-normal tail; mixed (purple) where a mixture of the two different behaviors, i.e., big large size avalanches and pop avalanches, are observed at the same time; explosive (red) where almost only big avalanches spanning all the network size are observed. For details about the detection of the transitions between regimes see text. The number of nodes of the network is N=10000 and the statistics is done out of n=10000 repetitions. Symbols indicate the points taken to show some distribution examples in the remaining panels. Star has coordinates (2.6,0.26), triangle (2,0.2), diamond (4,0.4) and doughnut (7,0.7); (C) Avalanche size distribution exactly at the transition between pop and mixed. At this transition the distribution is a power-law for a range of almost 3 decades, more precisely the distribution is a truncated power-law, given the finite size of the system. Dots show the simulation distribution and the dashed (resp. full) line shows the maximum of likelihood fit with a trucated power-law (resp. power-law) distribution; (D) Avalanche size distribution in the pop regime. Dots show the simulation distribution and the full (resp. dashed) line shows the maximum of likelihood fit with a log-normal (resp. truncated power-law) distribution; (E) Typical distribution in the mixed, regime, the inset shows a zoom on the distribution of large size avalanches not visible in the main panel; (F) Typical distribution in the explosive regime.

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 11022, 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 104. 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 log-normal 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 log-normal 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.

3.3 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.


FIGURE 4. Phase diagram of the resulting avalanche size distributions with respect to the model parameters k and p0, the constant pre-factor of the breaking probability (see Eq. 8). The time constant is τ=4 and the rheological exponent is α=0.3. Symbols indicate the points taken to show some distribution examples in the remaining panels. The pop (blue) and explosive (red) regimes are the same as in Figure 3, while the mixed regime changes, becoming composed of a broader bump, approximately log-normally distributed (in magenta), together with some small avalanches reminiscent of the pop regime. This regime is called here log-normal bump (magenta); (B) Example distribution of the log-normal bump regime with coordinates at the diamond point (4.5,0.55); (C) Example distribution of the pop regime, with coordinated at the triangle point (1.5,0.2). Dots show the simulation distribution and the full (resp. dashed) line shows the maximum of likelihood fit with a log-normal (resp. truncated power-law) distribution; (D) Example distribution of the explosive regime, with coordinates at the doughnut point (8,1).

The important difference introduced by a fractional visco-elastic 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 log-normal 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 p0 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 p0 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 p0 or k has to be large enough to quit the blue region. In other words we get log-normal 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.


FIGURE 5. (A) Phase diagram of the resulting avalanche size distributions with respect to the model parameters k and p0, the constant pre-factor of the breaking probability (see Eq. 8). The time constant is here τ=1 and the rheological exponent is again fixed at = 0.3. Symbols indicate the points taken to show some distribution examples in the remaining panels. The pop (blue) and explosive (red) regimes are the same as in Figure 3, while the log‐normal bump regime is composed of approximately log-normally distributed avalanche sizes (in magenta). Note that here this regime is much wider than in Figure 4 because of the choice of τ. (B) Example distribution of the log-normal bump regime with coordinates at the diamond point (4.5,0.95); (C) Example distribution of the pop regime, with coordinated at the triangle point (1.5,0.2). Dots show the simulation distribution and the full (resp. dashed) line shows the maximum of likelihood fit with a log-normal (resp. truncated power-law) distribution; (D) Example distribution of the explosive regime, with coordinates at the doughnut point (8,1).

We note that these new distributions are not exactly log-normal, 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 non-linear 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.

3.4 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 log-normal bump regime are robust and do not seem to depend on the details of the particular avalanche algorithm, provided a visco-elastic relaxation is included in the model. However these details can be important to tune the average and the width of the resulting distributions.


FIGURE 6. (A) Phase diagram of the resulting avalanche size distributions with respect to the model parameters k and p0 (see Eq. 8), in the case of link restoring. The time constant is τ=4 and the rheological exponent is α=0.3. Symbols indicate the points taken to show some distribution examples in the remaining panels. The pop (blue) and explosive (red) regimes are the same as in Figure 3, while the log‐normal bump regime is composed of approximately log‐normally distributed avalanche sizes (in magenta), together with some small avalanches reminiscent of the pop regime. Compared to Figure 3 the phase diagram does not show significant differences. (B) Example distribution of the log-normal bump regime with coordinates at the diamond point (4.5,0.55); (C) Example distribution of the pop regime, with coordinated at the triangle point (1.5,0.2). Dots show the simulation distribution and the full (resp. dashed) line shows the maximum of likelihood fit with a log-normal (resp. truncated power-law) distribution; (D) Example distribution of the explosive regime, with coordinates at the doughnut point (8,1).

4 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 visco-elastic 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.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

Conception and design: AA, FA, SP, and FP-R. Development and methodology: SP, AA, FA, and FP-R. Analysis and interpretation of data: SP, AA, FA, and FP-R. Writing, review, and/or revision of the manuscript: SP, FA, and FP-R. Administrative, technical or material support (i.e., requiring and organizing data, constructing databases): FA.


We acknowledge financial support from the Agence Nationale de la Recherche (ANR grant number ANR-18-CE45-0012-01) and from the French Research Minister (MESRI) for SP. PhD funding.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


We are very grateful for the participants of the second ISINP meeting at Lake Como for stimulating exchanges about our talk. We are thanful to P. Argoul, A. Guillet, E. Hartè, L. Delmarre for fruitful discussions. SP wishes to thank T. Matteuzzi for his inspiring considerations. We are indebted to Erika Polizzi for her graphical help.


1. Beggs, JM, and Plenz, D. Neuronal avalanches in neocortical circuits. J Neurosci (2003). 23:11167–77. doi:10.1523/JNEUROSCI.23-35-11167.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Polizzi, S, Laperrousaz, B, Perez-Reche, F, Nicolini, F, Satta, V, Arneodo, A, et al. A minimal rupture cascade model for living cell plasticity. New J Phys (2018). 20:053057. doi:10.1088/1367-2630/aac3c7

CrossRef Full Text | Google Scholar

3. Peŕez-Reche, FJ, Truskinovsky, L, and Zanzotto, G. Driving-induced crossover: from classical criticality to self-organized criticality. Phys Rev Lett (2008). 101:230601.

4. Salje, EK, Saxena, A, and Planes, A. eds. Avalanches in functional materials and geophysics. Cham, Switzerland: Springer Int. Publ. (2017).

Google Scholar

5. Bak, P. How nature works: the science of self-organized criticality. New York, NY: Springer Science and Business Media (2013). p. 212.

Google Scholar

6. Song, Y, Chen, X, Dabade, V, Shield, TW, and James, RD. Enhanced reversibility and unusual microstructure of a phase-transforming material. Nature (2013). 502:85–8. doi:10.1038/nature12532

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Chuang, CP, Yuan, T, Dmowski, W, Wang, GY, Freels, M, Liaw, PK, et al. Fatigue-induced damage in zr-based bulk metallic glasses. Sci Rep (2013). 3:2578. doi:10.1038/srep02578

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Pérez‐Reche, FJ, Triguero, C, Zanzotto, G, and Truskinovsky, L. Origin of scale-free intermittency in structural first‐order phase transitions. Phys Rev B. (2016). 94: 144102.

9. Perez‐Reche, FJ. Modeling avalanches in martensites. In: Avalanches in functional materials and geophysics. chap. 6 Dordrecht, Netherlands: Springer International Publishing (2017). p. 99–136.

Google Scholar

10. Lo, CC, Amaral, LN, Havlin, S, Ivanov, PC, Penzel, T, Peter, JH, et al. Dynamics of sleep-wake transitions during sleep. EPL (2002). 57:625. doi:10.1209/epl/i2002-00508-7

CrossRef Full Text | Google Scholar

11. Klaus, A, Yu, S, and Plenz, D. Statistical analyses support power law distributions found in neuronal avalanches. PloS One (2011). 6:e19779. doi:10.1371/journal.pone.0019779

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Dvir, H, Elbaz, I, Havlin, S, Appelbaum, L, Ivanov, PC, and Bartsch, RP. Neuronal noise as an origin of sleep arousals and its role in sudden infant death syndrome. Sci Adv (2018). 4:eaar6277. doi:10.1126/sciadv.aar6277

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wang, JW, Lombardi, F, Zhang, X, Anaclet, C, and Ivanov, PC. Non‐equilibrium critical dynamics of bursts in θ and δ rhythms as fundamental characteristic of sleep and wake micro-architecture. PLoS Comput. Biol. (2019). 15 (11): e1007268. doi:10.1371/journal.pcbi.1007268

14. Fabrizio, L, Manuel, G‐E, Pedro, B‐G, Ramalingam, V, Clifford, BS, Thomas, ES, et al. Critical dynamics and coupling in bursts of cortical rhythms indicate non-homeostatic mechanism for sleep‐stage transitions and dual role of VLPO neurons in both sleep and wake. J, Neurosci. (2020). 40 (1):171–190. doi:10.1523/JNEUROSCI.1278-19.2019

15. Broido, AD, and Clauset, A. Scale-free networks are rare. Nat Commun (2019). 10:1017. doi:10.1038/s41467-019-08746-5

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ribeiro, TL, Copelli, M, Caixeta, F, Belchior, H, Chialvo, DR, Nicolelis, MAL, et al. Spike avalanches exhibit universal dynamics across the sleep‐wake cycle. PLoS One. (2010). 5 (11):e14129. doi:10.1371/journal.pone.0014129

17. Buzsáki, G, and Mizuseki, K. The log-dynamic brain: how skewed distributions affect network operations. Nat Rev Neurosci (2014). 15:264–78. doi:10.1038/nrn3687

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Albert, R, Jeong, H, and Barabási, AL. Diameter of the world-wide web. Nature (1999). 401:130–1. doi:10.1038/43601

CrossRef Full Text | Google Scholar

19. Kesten, H. Random difference equations and Renewal theory for products of random matrices. Acta Math (1973). 131:207–48. doi:10.1007/BF02392040

CrossRef Full Text | Google Scholar

20. Polizzi, S. Emergence of log-normal distributions in avalanche processes, validation of 1-D stochastic and random network models, with an application to the characterization of cancer cells plasticity. [Ph.D. thesis]. Bordeaux (France): Univ. Bordeaux (2020).

Google Scholar

21. Barabási, AL, Albert, R, and Jeong, H. Mean-field theory for scale-free random networks. Phys Stat Mech Appl (1999). 272:173–87. doi:10.1016/S0378-4371(99)00291-5

CrossRef Full Text | Google Scholar

22. Barabási, AL, and Pósfai, M. Network science. Cambridge, United Kingdom.: Cambridge University Press (2016). p. 475.

Google Scholar

23. Loewenstein, Y, Kuras, A, and Rumpel, S. Multiplicative dynamics underlie the emergence of the log-normal distribution of spine sizes in the neocortex in vivo. J Neurosci (2011). 31:9481–8. doi:10.1523/JNEUROSCI.6130-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Legrand, J, Grais, RF, Boelle, PY, Valleron, AJ, and Flahault, A. Understanding the dynamics of ebola epidemics. Epidemiol Infect (2007). 135:610–21. doi:10.1017/S0950268806007217

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Jin, X, Zhang, Y, Wang, F, Li, L, Yao, D, Su, Y, et al. Departure headways at signalized intersections: a log-normal distribution model approach. Transport Res C Emerg Technol (2009). 17:318–27. doi:10.1016/j.trc.2009.01.003

CrossRef Full Text | Google Scholar

26. Gerstein, GL, and Mandelbrot, B. Random walk models for the spike activity of a single neuron. Biophys J (1964). 4:41–68. doi:10.1016/s0006-3495(64)86768-0

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ivanov, PC, Rosenblum, MG, Peng, CK, Mietus, J, Havlin, S, Stanley, HE, et al. Scaling behaviour of heartbeat intervals obtained by wavelet-based time-series analysis. Nature (1996). 383:323–7. doi:10.1038/383323a0

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ivanov, PC, Rosenblum, MG, Peng, CK, Mietus, JE, Havlin, S, Stanley, HE, et al. Scaling and universality in heart rate variability distributions. Physica A (1998). 249:587–93. doi:10.1016/s0378-4371(97)00522-0

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Laperrousaz, B, Drillon, G, Berguiga, L, Nicolini, F, Audit, B, Satta, VM, et al. From elasticity to inelasticity in cancer cell mechanics: a loss of scale-invariance. AIP Conf Proc (2016). 1760:020040. doi:10.1063/1.4960259

CrossRef Full Text | Google Scholar

30. Newman, ME, Strogatz, SH, and Watts, DJ. Random graphs with arbitrary degree distributions and their applications. Phys Rev E–Stat Nonlinear Soft Matter Phys (2001). 64:026118–7. doi:10.1103/PhysRevE.64.026118

CrossRef Full Text | Google Scholar

31. Newman, ME. Spread of epidemic disease on networks. Phys Rev E–Stat Nonlinear Soft Matter Phys (2002). 66:016128. doi:10.1103/PhysRevE.66.016128

CrossRef Full Text | Google Scholar

32. Ritort, F, and Sollich, P. Glassy dynamics of kinetically constrained models. Adv Phys (2003). 52:219–342. doi:10.1080/0001873031000093582

CrossRef Full Text | Google Scholar

33. Hoffman, BD, and Crocker, JC. Cell mechanics: dissecting the physical responses of cells to force. Annu Rev Biomed Eng (2009). 11:259–88. doi:10.1146/annurev.bioeng.10.061807.160511

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Chauvière, A, Preziosi, L, and Verdier, C. Cell mechanics: from single scale-based models to multiscale modeling. Boca Raton, FL: CRC Press (2010). p. 482.

Google Scholar

35. Brunner, C, Niendorf, A, and Käs, JA. Passive and active single-cell biomechanics: a new perspective in cancer diagnosis. Soft Matter (2009). 5:2171–8. doi:10.1039/B807545J

CrossRef Full Text | Google Scholar

36. Kai, F, Laklai, H, and Weaver, VM. Force matters: biomechanical regulation of cell invasion and migration in disease. Trends Cell Biol (2016). 26:486–97. doi:10.1016/j.tcb.2016.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gardel, ML, Kasza, KE, Brangwynne, CP, Liu, J, and Weitz, DA. Chapter 19: mechanical response of cytoskeletal networks. Methods Cell Biol (2008). 89:487–519. doi:10.1016/S0091-679X(08)00619-5

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lieleg, O, Claessens, MM, and Bausch, AR. Structure and dynamics of cross-linked actin networks. Soft Matter (2010). 6:218–25. doi:10.1039/b912163n

CrossRef Full Text | Google Scholar

39. Kollmannsberger, P, and Fabry, B. Linear and nonlinear rheology of living cells. Annu Rev Mater Res (2011). 41:75–97. doi:10.1146/annurev-matsci-062910-100351

CrossRef Full Text | Google Scholar

40. Fletcher, DA, and Mullins, RD. Cell mechanics and the cytoskeleton. Nature (2010). 463:485–92. doi:10.1038/nature08908

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Huber, F, Schnauß, J, Rönicke, S, Rauch, P, Müller, K, Fütterer, C, et al. Emergent complexity of the cytoskeleton: from single filaments to tissue. Adv Phys (2013). 62:1–112. doi:10.1080/00018732.2013.771509

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Käs, L, Boujemaa-Paterski, R, Sykes, C, and Plastino, J. Actin dynamics, architecture, and mechanics in cell motility. Physiol Rev (2014). 94:235–63. doi:10.1152/physrev.00018.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Li, R, and Gundersen, GG. Beyond polymer polarity: how the cytoskeleton builds a polarized cell. Nat Rev Mol Cell Biol (2008). 9:860–73. doi:10.1038/nrm2522

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Tseng, Y, Kole, TP, Lee, JS, Fedorov, E, Almo, SC, Schafer, BW, et al. How actin crosslinking and bundling proteins cooperate to generate an enhanced cell mechanical response. Biochem Biophys Res Commun (2005). 334:183–92. doi:10.1016/j.bbrc.2005.05.205

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wagner, B, Tharmann, R, Haase, I, Fischer, M, and Bausch, AR. Cytoskeletal polymer networks: the molecular structure of cross-linkers determines macroscopic properties. Proc Natl Acad Sci USA (2006). 103:13974–8. doi:10.1073/pnas.0510190103

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Ehrlicher, AJ, Krishnan, R, Guo, M, Bidan, CM, Weitz, DA, and Pollak, MR. Alpha-actinin binding kinetics modulate cellular dynamics and force generation. Proc Natl Acad Sci USA (2015). 112:6619–24. doi:10.1073/pnas.1505652112

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Fabry, B, Maksym, GN, Butler, JP, Glogauer, M, Navajas, D, and Fredberg, JJ. Scaling the microrheology of living cells. Phys Rev Lett (2001). 87:148102. doi:10.1103/PhysRevLett.87.148102

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Sollich, P. Rheological constitutive equation for a model of soft glassy materials. Phys Rev (1998). 58:738. doi:10.1103/PhysRevE.58.738

CrossRef Full Text | Google Scholar

49. Hildebrandt, J. Comparison of mathematical models for cat lung and viscoelastic balloon derived by Laplace transform methods from pressure-volume data. Bull Math Biophys (1969). 31:651–67. doi:10.1007/BF02477779

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Bonfanti, A, Kaplan, JL, Charras, G, and Kabla, A. Fractional viscoelastic models for power-law materials. Soft Matter (2020). 16:6002–20. doi:10.1039/d0sm00354a

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Streppa, L, Ratti, F, Goillot, E, Devin, A, Schaeffer, L, Arneodo, A, et al. Prestressed cells are prone to cytoskeleton failures under localized shear strain: an experimental demonstration on muscle precursor cells. Sci Rep (2018). 8:8602–16. doi:10.1038/s41598-018-26797-4

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Trepat, X, Grabulosa, M, Buscemi, L, Rico, F, Farré, R, and Navajas, D. Thrombin and histamine induce stiffening of alveolar epithelial cells. J Appl Physiol (2005). 98:1567–74. doi:10.1152/japplphysiol.00925.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Gerasimova-Chechkina, E, Streppa, L, Schaeffer, L, Devin, A, Argoul, P, Arneodo, A, et al. Fractional rheology of muscle precursor cells. J Rheol (2018). 62:1347–62. doi:10.1122/1.5035127

CrossRef Full Text | Google Scholar

54. Khalilgharibi, N, Fouchard, J, Asadipour, N, Barrientos, R, Duda, M, Bonfanti, A, et al. Stress relaxation in epithelial monolayers is controlled by the actomyosin cortex. Nat Phys (2019). 15:839–47. doi:10.1038/s41567-019-0516-6

CrossRef Full Text | Google Scholar

55. Carmichael, B, Babahosseini, H, Mahmoodi, SN, and Agah, M. The fractional viscoelastic response of human breast tissue cells. Phys Biol (2015). 12:046001. doi:10.1088/1478-3975/12/4/046001

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Gorenflo, R, Loutchko, J, and Luchko, Y. Computation of the mittag-leffler function eα, β (z) and its derivative. Fract Calc Appl Anal (2002). 5:491–512.

Google Scholar

57. Baró, J, and Davidsen, J. Universal avalanche statistics and triggering close to failure in a mean-field model of rheological fracture. Phys Rev (2018). 97:033002. doi:10.1103/PhysRevE.97.033002

CrossRef Full Text | Google Scholar

58. Juelicher, F, Kruse, K, Prost, J, and Joanny, JF. Active behavior of the cytoskeleton. Phys Rep (2007). 449:3–28. doi:10.1016/j.physrep.2007.02.018

CrossRef Full Text | Google Scholar

59. Dos Remedios, CG, Chhabra, D, Kekic, M, Dedova, IV, Tsubakihara, M, Berry, DA, et al. Actin binding proteins: regulation of cytoskeletal microfilaments. Physiol Rev (2003). 83:433–73. doi:10.1152/physrev.00026.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Esue, O, Tseng, Y, and Wirtz, D. Alpha-actinin and filamin cooperatively enhance the stiffness of actin filament networks. PloS One (2009). 4:e4411. doi:10.1371/journal.pone.0004411

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Wolff, L, Fernandez, P, and Kroy, K. Inelastic mechanics of sticky biopolymer networks. New J Phys (2010). 12:053024. doi:10.1088/1367-2630/12/5/053024

CrossRef Full Text | Google Scholar

62. Gralka, M, and Kroy, K. Inelastic mechanics: a unifying principle in biomechanics. Biochim Biophys Acta (2015). 1853:3025–37. doi:10.1016/j.bbamcr.2015.06.017

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Badre, D, and Wagner, AD. Computational and neurobiological mechanisms underlying cognitive flexibility. Proc Natl Acad Sci USA (2006). 103:7186–91. doi:10.1073/pnas.0509550103

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Bonakdar, N, Gerum, R, Kuhn, M, Spörrer, M, Lippert, A, Schneider, W, et al. Mechanical plasticity of cells. Nat Mater (2016). 15:1090–4. doi:10.1038/nmat4689

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Ising, E. Beitrag zur theorie des ferromagnetismus. Z Phys (1925). 31:253–8. doi:10.1007/BF02980577

CrossRef Full Text | Google Scholar

66. Laperrousaz, B, Berguiga, L, Nicolini, FE, Martinez-Torres, C, Arneodo, A, Satta, VM, et al. Revealing stiffening and brittling of chronic myelogenous leukemia hematopoietic primary cells through their temporal response to shear stress. Phys Biol (2016). 13:03LT01. doi:10.1088/1478-3975/13/3/03LT01

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Streppa, L. Characterizing mechanical properties of living C2C12 myoblasts with single cell indentation experiments: application to Duchenne muscular dystrophy. [Ph.D. thesis]. Lyon (France): University of Lyon (2017).

Google Scholar

68. Grossmann, A, and Morlet, J. Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM J Math Anal (1984). 15:723–36. doi:10.1137/0515056

CrossRef Full Text | Google Scholar

69. Digiuni, S, Berne-Dedieu, A, Martinez-Torres, C, Szecsi, J, Bendahmane, M, Arneodo, A, et al. Single cell wall nonlinear mechanics revealed by a multiscale analysis of AFM force-indentation curves. Biophys J (2015). 108:2235–48. doi:10.1016/j.bpj.2015.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Martinez‐Torres, C, Arneodo, A, Streppa, L, Argoul, P, and Argoul, F. Passive microrheology of soft materials with atomic force microscopy: a wavelet-based spectral analysis. Appl Phys Lett (2016). 108:034102. doi:10.1063/1.4940220

CrossRef Full Text | Google Scholar

71. Sneddon, IN. The relation between load and penetration in the axisymmetric boussinesq problem for a punch of arbitrary profile. Int J Eng Sci (1965). 3:47–57. doi:10.1016/0020-7225(65)90019-4

CrossRef Full Text | Google Scholar

72. Eghiaian, F, Rigato, A, and Scheuring, S. Structural, mechanical, and dynamical variability of the actin cortex in living cells. Biophys J (2015). 108:1330–40. doi:10.1016/j.bpj.2015.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Pastor-Satorras, R, and Vespignani, A. Epidemic spreading in scale-free networks. Phys Rev Lett (2001). 86:3200. doi:10.1103/PhysRevLett.86.3200

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Gilbert, EN. Random graphs. Ann Math Stat (1959). 30:1141–4.

CrossRef Full Text | Google Scholar

75. Erdös, P, and Rényi, A. On random graphs i. Publ Math (1959). 6:290–7.

Google Scholar

76. Alstott, J, Bullmore, E, and Plenz, D. Powerlaw: a python package for analysis of heavy-tailed distributions. PloS One (2014). 9:e85777. doi:10.1371/journal.pone.0085777

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Severini, T. Likelihood methods in statistics. New York, NY: Oxford Univ. Press (2001). p. 105–37.

Google Scholar

78. Clauset, A, Shalizi, CR, and Newman, ME. Power-law distributions in empirical data. SIAM Rev (2009). 51:661–703. doi:10.1137/070710111

CrossRef Full Text | Google Scholar

Keywords: random network, avalanches, log-normal distributions, power-law, cell plasticity, cytoskeleton ruptures

Citation: Polizzi S, Arneodo A, Pérez-Reche F-J and Argoul F (2021) Emergence of Log-Normal Type Distributions in Avalanche Processes in Living Systems: A Network Model. Front. Appl. Math. Stat. 6:613962. doi: 10.3389/fams.2020.613962

Received: 04 October 2020; Accepted: 23 December 2020;
Published: 28 January 2021.

Edited by:

Plamen Ch. Ivanov, Boston University, United States

Reviewed by:

Changgui Gu, University of Shanghai for Science and Technology, China
Chengyu Huo, Changshu Institute of Technology, China

Copyright © 2021 Polizzi, Arneodo, Pérez-Reche and Argoul. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Françoise Argoul,