Abstract
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 [], fractures in living cells cytoskeleton (CSK) [], but also in amorphous [] and random media [], 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 [], and are ubiquitous for avalanches in solids and amorphous materials [], describing for example the avalanche sizes of a granular material falling apart. Self-organized criticality has also been associated to neural dynamics (first in [] and later in [, ]), showing widely spread power-law statistics with a typical exponent invariant among different species [] 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 [, ]. 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 [], 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 [], but also for avalanche statistics, in neural networks [, ] or living cells cytoskeleton []. Indeed, the growth of a network can be seen itself as an avalanche process, thinking for example of the preferential attachment rule []. 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 [, ].
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 [], justifying a power-law interpretation of real data distributions, while models for other fat tail distributions do not exist so far [, ]; 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 [].
FIGURE 1
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 [
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 [
Originally motivated by the experimental observation of log-normal statistics in avalanches of fractures of the CSK of living cells [
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 [
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 [
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 [
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) [
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 [
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) [
The ability of cells to switch quickly from fluid-like responses to more brittle solid-like responses [
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 [
2.2 Poking Living Cells With a Sharp Atomic Force Microscopy Tip
Previous experiments done in our group motivated this study [
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.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 (Figure 2A), caused by the avalanche of failures, and the indentation length (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 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 , 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 (see Figure 2B,D), extracted by the identification of the whole force indentation curve with the Sneddon model [
FIGURE 2

(A) Probability distributions of the logarithm of the local force drops 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 estimated by a parabolic fit of the Force-Indentation Curves (Sneddon model [
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 [
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 [
For all these reasons, we propose a Gilbert [
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 links coming out from a node. The probability 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 , keeping fixed, the binomial degree distribution is well approximated by the simpler (one parameter) Poisson distribution, noted , with parameter . The average connectivity 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 with the number of all the possible pairs (by excluding double counting):and then for a network of nodes, and , the average number of links is . 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 [
At this point we can work out the probability 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 and the average number of links of the network .
On this network structure, an avalanche of ruptures is described by the following algorithm. Notice that while the probability of connection
remains constant, the probability of breaking
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 , where at the first time step ;
• we consider the links coming out from the ones broken at the previous time step from both extremities and break them with probability .
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, , 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 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 (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 . The choice of the Poisson distribution is motivated by the degree distribution, which, as we said, is approximately Poissonian in the limit keeping fixed, and the parameter 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 , up to , 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 [
We considered 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 . For for instance is already irrelevant for our results, and for , , 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 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 [
3.2 Avalanche Statistics With Constant Probability of Breaking
First we implemented the avalanche model with a probability of breaking 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 nodes and a probability of connection ; (B) Phase diagram of the possible resulting avalanche size distributions with respect to the model parameters (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 and the statistics is done out of repetitions. Symbols indicate the points taken to show some distribution examples in the remaining panels. Star has coordinates , triangle , diamond and doughnut ; (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 [
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 , 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 [
The boundary between the pop and mixed regimes is computed by an algorithm detecting for which parameters 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 [
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 ) 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 is equivalent to setting , since in both cases 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 . 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 and , the constant pre-factor of the breaking probability (see Eq. 8). The time constant is and the rheological exponent is . 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 ; (C) Example distribution of the pop regime, with coordinated at the triangle point . 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 .
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 become larger than 1, since the stretched exponential decay will always make converge toward 0. Simulations in Figure 4 are shown for , 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 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, (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 or 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 or 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 . 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 and , the constant pre-factor of the breaking probability (see Eq. 8). The time constant is here 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 ; (C) Example distribution of the pop regime, with coordinated at the triangle point . 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 .
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 [
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 (see Eq. 7). We thus take the same law of breaking as in Eq. 8, with .
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 and (see Eq. 8), in the case of link restoring. The time constant is and the rheological exponent is . 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 ; (C) Example distribution of the pop regime, with coordinated at the triangle point . 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 .
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 [
Funding
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.
Statements
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.
Acknowledgments
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.
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.
References
1.
BeggsJMPlenzD. Neuronal avalanches in neocortical circuits. J Neurosci (2003). 23:11167–77. 10.1523/JNEUROSCI.23-35-11167.2003
2.
PolizziSLaperrousazBPerez-RecheFNicoliniFSattaVArneodoAet alA minimal rupture cascade model for living cell plasticity. New J Phys (2018). 20:053057. 10.1088/1367-2630/aac3c7
3.
Peŕez-RecheFJTruskinovskyLZanzottoG. Driving-induced crossover: from classical criticality to self-organized criticality. Phys Rev Lett (2008). 101:230601.
4.
SaljeEKSaxenaAPlanesA. eds. Avalanches in functional materials and geophysics. Cham, Switzerland: Springer Int. Publ. (2017).
5.
BakP. How nature works: the science of self-organized criticality. New York, NY: Springer Science and Business Media (2013). p. 212.
6.
SongYChenXDabadeVShieldTWJamesRD. Enhanced reversibility and unusual microstructure of a phase-transforming material. Nature (2013). 502:85–8. 10.1038/nature12532
7.
ChuangCPYuanTDmowskiWWangGYFreelsMLiawPKet alFatigue-induced damage in zr-based bulk metallic glasses. Sci Rep (2013). 3:2578. 10.1038/srep02578
8.
Pérez‐RecheFJTrigueroCZanzottoGTruskinovskyL. Origin of scale-free intermittency in structural first‐order phase transitions. Phys Rev B. (2016). 94: 144102.
9.
Perez‐RecheFJ. Modeling avalanches in martensites. In: Avalanches in functional materials and geophysics. chap. 6 Dordrecht, Netherlands: Springer International Publishing (2017). p. 99–136.
10.
LoCCAmaralLNHavlinSIvanovPCPenzelTPeterJHet alDynamics of sleep-wake transitions during sleep. EPL (2002). 57:625. 10.1209/epl/i2002-00508-7
11.
KlausAYuSPlenzD. Statistical analyses support power law distributions found in neuronal avalanches. PloS One (2011). 6:e19779. 10.1371/journal.pone.0019779
12.
DvirHElbazIHavlinSAppelbaumLIvanovPCBartschRP. Neuronal noise as an origin of sleep arousals and its role in sudden infant death syndrome. Sci Adv (2018). 4:eaar6277. 10.1126/sciadv.aar6277
13.
WangJWLombardiFZhangXAnacletCIvanovPC. 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. 10.1371/journal.pcbi.1007268
14.
FabrizioLManuelG‐EPedroB‐GRamalingamVCliffordBSThomasESet alCritical 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. 10.1523/JNEUROSCI.1278-19.2019
15.
BroidoADClausetA. Scale-free networks are rare. Nat Commun (2019). 10:1017. 10.1038/s41467-019-08746-5
16.
RibeiroTLCopelliMCaixetaFBelchiorHChialvoDRNicolelisMALet alSpike avalanches exhibit universal dynamics across the sleep‐wake cycle. PLoS One. (2010). 5 (11):e14129. 10.1371/journal.pone.0014129
17.
BuzsákiGMizusekiK. The log-dynamic brain: how skewed distributions affect network operations. Nat Rev Neurosci (2014). 15:264–78. 10.1038/nrn3687
18.
AlbertRJeongHBarabásiAL. Diameter of the world-wide web. Nature (1999). 401:130–1. 10.1038/43601
19.
KestenH. Random difference equations and Renewal theory for products of random matrices. Acta Math (1973). 131:207–48. 10.1007/BF02392040
20.
PolizziS. 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).
21.
BarabásiALAlbertRJeongH. Mean-field theory for scale-free random networks. Phys Stat Mech Appl (1999). 272:173–87. 10.1016/S0378-4371(99)00291-5
22.
BarabásiALPósfaiM. Network science. Cambridge, United Kingdom.: Cambridge University Press (2016). p. 475.
23.
LoewensteinYKurasARumpelS. Multiplicative dynamics underlie the emergence of the log-normal distribution of spine sizes in the neocortex in vivo. J Neurosci (2011). 31:9481–8. 10.1523/JNEUROSCI.6130-10.2011
24.
LegrandJGraisRFBoellePYValleronAJFlahaultA. Understanding the dynamics of ebola epidemics. Epidemiol Infect (2007). 135:610–21. 10.1017/S0950268806007217
25.
JinXZhangYWangFLiLYaoDSuYet alDeparture headways at signalized intersections: a log-normal distribution model approach. Transport Res C Emerg Technol (2009). 17:318–27. 10.1016/j.trc.2009.01.003
26.
GersteinGLMandelbrotB. Random walk models for the spike activity of a single neuron. Biophys J (1964). 4:41–68. 10.1016/s0006-3495(64)86768-0
27.
IvanovPCRosenblumMGPengCKMietusJHavlinSStanleyHEet alScaling behaviour of heartbeat intervals obtained by wavelet-based time-series analysis. Nature (1996). 383:323–7. 10.1038/383323a0
28.
IvanovPCRosenblumMGPengCKMietusJEHavlinSStanleyHEet alScaling and universality in heart rate variability distributions. Physica A (1998). 249:587–93. 10.1016/s0378-4371(97)00522-0
29.
LaperrousazBDrillonGBerguigaLNicoliniFAuditBSattaVMet alFrom elasticity to inelasticity in cancer cell mechanics: a loss of scale-invariance. AIP Conf Proc (2016). 1760:020040. 10.1063/1.4960259
30.
NewmanMEStrogatzSHWattsDJ. Random graphs with arbitrary degree distributions and their applications. Phys Rev E–Stat Nonlinear Soft Matter Phys (2001). 64:026118–7. 10.1103/PhysRevE.64.026118
31.
NewmanME. Spread of epidemic disease on networks. Phys Rev E–Stat Nonlinear Soft Matter Phys (2002). 66:016128. 10.1103/PhysRevE.66.016128
32.
RitortFSollichP. Glassy dynamics of kinetically constrained models. Adv Phys (2003). 52:219–342. 10.1080/0001873031000093582
33.
HoffmanBDCrockerJC. Cell mechanics: dissecting the physical responses of cells to force. Annu Rev Biomed Eng (2009). 11:259–88. 10.1146/annurev.bioeng.10.061807.160511
34.
ChauvièreAPreziosiLVerdierC. Cell mechanics: from single scale-based models to multiscale modeling. Boca Raton, FL: CRC Press (2010). p. 482.
35.
BrunnerCNiendorfAKäsJA. Passive and active single-cell biomechanics: a new perspective in cancer diagnosis. Soft Matter (2009). 5:2171–8. 10.1039/B807545J
36.
KaiFLaklaiHWeaverVM. Force matters: biomechanical regulation of cell invasion and migration in disease. Trends Cell Biol (2016). 26:486–97. 10.1016/j.tcb.2016.03.007
37.
GardelMLKaszaKEBrangwynneCPLiuJWeitzDA. Chapter 19: mechanical response of cytoskeletal networks. Methods Cell Biol (2008). 89:487–519. 10.1016/S0091-679X(08)00619-5
38.
LielegOClaessensMMBauschAR. Structure and dynamics of cross-linked actin networks. Soft Matter (2010). 6:218–25. 10.1039/b912163n
39.
KollmannsbergerPFabryB. Linear and nonlinear rheology of living cells. Annu Rev Mater Res (2011). 41:75–97. 10.1146/annurev-matsci-062910-100351
40.
FletcherDAMullinsRD. Cell mechanics and the cytoskeleton. Nature (2010). 463:485–92. 10.1038/nature08908
41.
HuberFSchnaußJRönickeSRauchPMüllerKFüttererCet alEmergent complexity of the cytoskeleton: from single filaments to tissue. Adv Phys (2013). 62:1–112. 10.1080/00018732.2013.771509
42.
KäsLBoujemaa-PaterskiRSykesCPlastinoJ. Actin dynamics, architecture, and mechanics in cell motility. Physiol Rev (2014). 94:235–63. 10.1152/physrev.00018.2013
43.
LiRGundersenGG. Beyond polymer polarity: how the cytoskeleton builds a polarized cell. Nat Rev Mol Cell Biol (2008). 9:860–73. 10.1038/nrm2522
44.
TsengYKoleTPLeeJSFedorovEAlmoSCSchaferBWet alHow actin crosslinking and bundling proteins cooperate to generate an enhanced cell mechanical response. Biochem Biophys Res Commun (2005). 334:183–92. 10.1016/j.bbrc.2005.05.205
45.
WagnerBTharmannRHaaseIFischerMBauschAR. Cytoskeletal polymer networks: the molecular structure of cross-linkers determines macroscopic properties. Proc Natl Acad Sci USA (2006). 103:13974–8. 10.1073/pnas.0510190103
46.
EhrlicherAJKrishnanRGuoMBidanCMWeitzDAPollakMR. Alpha-actinin binding kinetics modulate cellular dynamics and force generation. Proc Natl Acad Sci USA (2015). 112:6619–24. 10.1073/pnas.1505652112
47.
FabryBMaksymGNButlerJPGlogauerMNavajasDFredbergJJ. Scaling the microrheology of living cells. Phys Rev Lett (2001). 87:148102. 10.1103/PhysRevLett.87.148102
48.
SollichP. Rheological constitutive equation for a model of soft glassy materials. Phys Rev (1998). 58:738. 10.1103/PhysRevE.58.738
49.
HildebrandtJ. 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. 10.1007/BF02477779
50.
BonfantiAKaplanJLCharrasGKablaA. Fractional viscoelastic models for power-law materials. Soft Matter (2020). 16:6002–20. 10.1039/d0sm00354a
51.
StreppaLRattiFGoillotEDevinASchaefferLArneodoAet alPrestressed cells are prone to cytoskeleton failures under localized shear strain: an experimental demonstration on muscle precursor cells. Sci Rep (2018). 8:8602–16. 10.1038/s41598-018-26797-4
52.
TrepatXGrabulosaMBuscemiLRicoFFarréRNavajasD. Thrombin and histamine induce stiffening of alveolar epithelial cells. J Appl Physiol (2005). 98:1567–74. 10.1152/japplphysiol.00925.2004
53.
Gerasimova-ChechkinaEStreppaLSchaefferLDevinAArgoulPArneodoAet alFractional rheology of muscle precursor cells. J Rheol (2018). 62:1347–62. 10.1122/1.5035127
54.
KhalilgharibiNFouchardJAsadipourNBarrientosRDudaMBonfantiAet alStress relaxation in epithelial monolayers is controlled by the actomyosin cortex. Nat Phys (2019). 15:839–47. 10.1038/s41567-019-0516-6
55.
CarmichaelBBabahosseiniHMahmoodiSNAgahM. The fractional viscoelastic response of human breast tissue cells. Phys Biol (2015). 12:046001. 10.1088/1478-3975/12/4/046001
56.
GorenfloRLoutchkoJLuchkoY. Computation of the mittag-leffler function eα, β (z) and its derivative. Fract Calc Appl Anal (2002). 5:491–512.
57.
BaróJDavidsenJ. Universal avalanche statistics and triggering close to failure in a mean-field model of rheological fracture. Phys Rev (2018). 97:033002. 10.1103/PhysRevE.97.033002
58.
JuelicherFKruseKProstJJoannyJF. Active behavior of the cytoskeleton. Phys Rep (2007). 449:3–28. 10.1016/j.physrep.2007.02.018
59.
Dos RemediosCGChhabraDKekicMDedovaIVTsubakiharaMBerryDAet alActin binding proteins: regulation of cytoskeletal microfilaments. Physiol Rev (2003). 83:433–73. 10.1152/physrev.00026.2002
60.
EsueOTsengYWirtzD. Alpha-actinin and filamin cooperatively enhance the stiffness of actin filament networks. PloS One (2009). 4:e4411. 10.1371/journal.pone.0004411
61.
WolffLFernandezPKroyK. Inelastic mechanics of sticky biopolymer networks. New J Phys (2010). 12:053024. 10.1088/1367-2630/12/5/053024
62.
GralkaMKroyK. Inelastic mechanics: a unifying principle in biomechanics. Biochim Biophys Acta (2015). 1853:3025–37. 10.1016/j.bbamcr.2015.06.017
63.
BadreDWagnerAD. Computational and neurobiological mechanisms underlying cognitive flexibility. Proc Natl Acad Sci USA (2006). 103:7186–91. 10.1073/pnas.0509550103
64.
BonakdarNGerumRKuhnMSpörrerMLippertASchneiderWet alMechanical plasticity of cells. Nat Mater (2016). 15:1090–4. 10.1038/nmat4689
65.
IsingE. Beitrag zur theorie des ferromagnetismus. Z Phys (1925). 31:253–8. 10.1007/BF02980577
66.
LaperrousazBBerguigaLNicoliniFEMartinez-TorresCArneodoASattaVMet alRevealing stiffening and brittling of chronic myelogenous leukemia hematopoietic primary cells through their temporal response to shear stress. Phys Biol (2016). 13:03LT01. 10.1088/1478-3975/13/3/03LT01
67.
StreppaL. 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).
68.
GrossmannAMorletJ. Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM J Math Anal (1984). 15:723–36. 10.1137/0515056
69.
DigiuniSBerne-DedieuAMartinez-TorresCSzecsiJBendahmaneMArneodoAet alSingle cell wall nonlinear mechanics revealed by a multiscale analysis of AFM force-indentation curves. Biophys J (2015). 108:2235–48. 10.1016/j.bpj.2015.02.024
70.
Martinez‐TorresCArneodoAStreppaLArgoulPArgoulF. Passive microrheology of soft materials with atomic force microscopy: a wavelet-based spectral analysis. Appl Phys Lett (2016). 108:034102. 10.1063/1.4940220
71.
SneddonIN. 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. 10.1016/0020-7225(65)90019-4
72.
EghiaianFRigatoAScheuringS. Structural, mechanical, and dynamical variability of the actin cortex in living cells. Biophys J (2015). 108:1330–40. 10.1016/j.bpj.2015.01.016
73.
Pastor-SatorrasRVespignaniA. Epidemic spreading in scale-free networks. Phys Rev Lett (2001). 86:3200. 10.1103/PhysRevLett.86.3200
74.
GilbertEN. Random graphs. Ann Math Stat (1959). 30:1141–4.
75.
ErdösPRényiA. On random graphs i. Publ Math (1959). 6:290–7.
76.
AlstottJBullmoreEPlenzD. Powerlaw: a python package for analysis of heavy-tailed distributions. PloS One (2014). 9:e85777. 10.1371/journal.pone.0085777
77.
SeveriniT. Likelihood methods in statistics. New York, NY: Oxford Univ. Press (2001). p. 105–37.
78.
ClausetAShaliziCRNewmanME. Power-law distributions in empirical data. SIAM Rev (2009). 51:661–703. 10.1137/070710111
Summary
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
Volume
6 - 2020
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
Updates

Check for updates
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, francoise.argoul@u-bordeaux.fr
This article was submitted to Dynamical Systems, a section of the journal Frontiers in Applied Mathematics and Statistics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.