Percolation in Biological Systems

Sahimi, Muhammad

doi:10.1007/978-3-031-20386-2_18

Muhammad Sahimi¹⁸

Part of the book series: Applied Mathematical Sciences ((AMS,volume 213))

801 Accesses
2 Citations

Abstract

Most biological systems are so complex that they preclude any reasonable description in terms of the exact interactions among their fundamental constituents. We do know, however, that many of the most fascinating biological phenomena, and in particular those involving life, are the result of interactions between many elements.

Download chapter PDF

18.1 Introduction

Most biological systems are so complex that they preclude any reasonable description in terms of the exact interactions among their fundamental constituents. We do know, however, that many of most fascinating biological phenomena, and in particular those involving life, are the result of interactions between many elements. For example (Mora and Bialek 2011), many amino acids determine the structure of a single protein; the fate of a cell is determined by many genes, and shaping our thoughts and memories involve many neurons. But, for quite some time, application of statistical mechanics of disordered media, and in particular critical phenomena and percolation theory, to biological materials and phenomena seemed intractable, even though, due to their nature and the fact that they are driven by collective behaviors, the application seemed natural. But, considerable progress has been made over the past two decades, as it has been shown, through experiments and simulations, that a surprising number of such phenomena manifest behavior and properties that are akin to systems near a critical point. For example, as described by Krause and Ruxton (2002), groups of animals, such as flocking birds, move with fascinating coordination in that, instead of being dictated by a leader or responding to a common stimulus, the collective patterns of flock dynamics is self-organized, arising from local interactions between individuals who propagate information throughout the entire group. Flocks of birds are also highly responsive and cohesive in the face of predatory threat. The balance between order and high susceptibility indicates that a flock of birds may represent a critical state. Indeed, this has been demonstrated (Bialek et al. 2012).

Percolation or connectivity threshold also represents a critical point. There are many biological systems that are seemingly stochastic in nature and in which the role of connectivity of the microscopic elements or constituents is prominent. Examples include self-assembly of tobacco mosaic and other simple viruses (see, for example, Hohn and Hohn 1970), actin filaments (Poglazov et al. 1967) and flagella (Asakura et al. 1968), lymphocyte patch and cap formation (Karnovsky et al. 1972), many precipitation and agglutination phenomena, and many other biological phenomena (see below). Some of such phenomena, such as precipitation, occur spontaneously if the functional groups are sufficiently reactive. Thus, they depend on their level of chemical complexity and that of the solvent in which they occur. Other factors are not directly related to the solvent, but have a great influence on the outcome of biological processes. For example, in antigen–antibody reactions, clusters of all sizes react with one another, forming complex branched networks that grow in size as time progresses. We may also have reactions that can proceed by rapid addition of monomers to growing chains after a slow initiating event. Such phenomena are, therefore, similar to percolation processes. This chapter discusses applications of percolation concepts to such phenomena. We start with a discussion of antigen–antibody reactions and aggregations, and then describe several other biological processes to which percolation may be relevant.

18.2 Antigen–Antibody Reactions and Aggregation

Under suitable conditions, mixtures of antigen and antibodies—a system comprised of bifunctional and multifunctional monomers—form large networks or aggregates that contain many molecules of both. In general, the networks or branched structures react reversibly with one another. They are often insoluble and, hence, will precipitate. Since we are interested in the concentrations of antibodies and antigens as a function of time, it is important to examine the conditions under which precipitation may be expected to occur.

The first theory of antigen–antibody aggregation was developed by Goldberg (1952) who used combinatorial mathematics to derive results that are equivalent to those obtained for various percolation properties on a Bethe lattice. His work was extended by several others, including Aladjem and Palmiter (1965), Bell (1971), and Delisi (1974). Delisi and Perelson (1976) analyzed the problem extensively and derived several analytical results for the various properties of interest.

Consider aggregation of two types of monomers, S and L, of functionality 2 and $Z\ge 2$, respectively. The S monomers can be any of a number of substrates, such as immunoglobulin G (IgG), while the L monomers might be some appropriate ligand or antigen. The functional units on the monomers are equivalent to the sites in percolation on random networks described and utilized so far in this book. To analyze the problem, Goldberg (1952) made three basic assumptions: (a) all the free sites are equivalent; (b) no intramolecular reactions can occur; and (c) during the aggregation process cyclization reactions do not occur. Since closed loops are not formed, it immediately points to the similarity with percolation on Bethe lattices. Delisi and Perelson (1976) relaxed the second assumption to allow for intramolecular reactions. In the immunological literature, such reactions are often referred to as monogamous bivalency.

Several parameters were then introduced: (a) $p_b$, the probability that an antibody site picked at random is bound; (b) $p_g$, the probability that an antigen site picked at random is bound; (c) $p_{br}$, the probability that an antibody site is bound and has not reacted intramolecularly; and (d) $p_{gr}$, the probability that an antigen site is bound and has not reacted intramolecularly. Consider, then, an aggregate composed of x Z-functional antigens and $(x-1+n)$ bifunctional antibodies, where n is the number of antibodies that do not serve as connectors. Of these antibodies, $n_1$ are bound univalently, and $n_2$ bivalently. An aggregate with compositions x, $n_1$, and $n_2$ is called an $(x,n_1,n_2)$-mer.

The main quantity of interest is $C_{n_1,n_2}(a, p)$, the concentration of such aggregates. To determine it, we define $P(x,n_1,n_2)$ as the probability that a free site picked at random is part of an $(x,n_,n_2)$-mer. If closed loops are not formed, then

$$\begin{aligned} P(x,n_1,n_2)=\frac{\text{ the } \text{ number } \text{ of } \text{ free } \text{ site } \text{ on } (x,n_1,n_2)\text{-mers }}{\text{ Total } \text{ number } \text{ of } \text{ free } \text{ sites }}\;, \end{aligned}$$

(18.1)

which implies that

$$\begin{aligned} P(x,n_1,n_2)=\frac{(Zx-2x+2-2n_2)C_{n_1,n_2}(x,p)}{L+S}\;. \end{aligned}$$

(18.2)

We write $P(x,n_1,n_2)$ as

$$\begin{aligned} P(x,n_1,n_2)=\rho _l\omega _l\Omega _l+\rho _s\Omega _s\Omega _s\;, \end{aligned}$$

(18.3)

where l and s label, respectively, free antigens and antibody sites. Here, $\rho _l$ is the probability that the root is a free antigen site, $\omega _l$ the probability that the root is on an $(x,n_1,n_2)$-mer, given that it is a free L site, and $\Omega _l$ is the number of ways an $(x,n_1,n_2)$-mer can form, given that the root is a free L site, with analogous definitions for $\rho _s$, $\omega _s$, and $\Omega _s$. Each of these quantities can be calculated using the methods developed by Flory and Stockmayer for their theory of gelation described in Chap. 13, or by the methods developed for percolation on Bethe lattices. If closed loops are forbidden, then Goldberg (1952) showed that for

$$\begin{aligned} p=p_g =\frac{1}{Z-1}\;, \end{aligned}$$

(18.4)

there is a finite probability that an infinitely large aggregate forms. This is exactly the same as the site (or bond) percolation threshold of a Bethe lattice of coordination number Z (see Chap. 2).

Similar to the gelation and polymerization, the problem with such formulation is that the formation of closed loops is forbidden. Thus, the results that were derived represent the mean-field limit of percolation, i.e., a system whose spatial dimension is six or higher. In reality, antigen–antibody aggregation is a three-dimensional (3D) phenomenon and, therefore, we must use percolation on finite-dimensional lattices, and in particular 3D ones, in order to explain antigen–antibody aggregation phenomena.

18.3 Network Formation on Lymphocyte Membranes

According to modern theory of the clonal selection, the cells of the immune systems arise from the division and differentiation of the stem cells in the bone marrow (Burnet 1959). The cells that are potentially able to secrete antibodies are known as B-lymphocytes, which insert a homogeneous set of antibody-like receptor molecules on their membrane’s surface, where they are used for recognition of complementary antigens. The population of B-cells is itself heterogeneous, i.e., different cells may have different types of membrane-bound receptor molecules. Under appropriate conditions, binding of complementary antigen to the receptors on a B-cell activates that cell to differentiate and/or proliferate into antibody-secreting cells. Moreover, the antibody secreted by the progeny of a particular B-cell is assumed to have specificity for the antigen, identical to that of the receptor molecules on the B-cell progenitor. Therefore, the role of antigen is cellular selection and amplification, and its mission is mediated by interaction with cellular receptors.

The process of proliferation and secreting large amounts of antibody is preceded by cross-linking of antigens to the receptors until a macroscopically large patch is first formed on the surface. The patch or network is the 2D analogue of 3D aggregation, described in the preceding section, a process called network or lattice formation. Antigen-stimulated lymphocyte lattice formation may play a key role in triggering immunocompetent cells. For immunogens with repeating arrays of antigenic determinants, e.g., polysaccharide flagellin, triggering is possible without the aid of helper cells, whereas other immunogens require the presence of T-cells or their products. It has been suggested that the T-cell requirement results from having to present the immunogen to the B-cells in an aggregated or network form, thus increasing its valence.

The relationship between network formation and biological activity is not completely clear yet. It is presumed that lattice formation provides some sort of stimulatory signal, but this may not be enough to trigger a cell. In any event, it is not unreasonable to assume that the triggering is a function of the strength of antibody–receptor interaction. Delisi and Perelson (1976) assumed that lattice formation is the 2D analogue of the 3D precipitation reaction described above and analyzed the problem in a similar way. To simplify the problem, they assumed that the effect of intramolecular reactions is unimportant. A well-understood characteristic of immune response to antigens is an increase in antigen–antibody affinity with time after immunization. The implication is that some aspects of the antigen–receptor interaction are equilibrium-controlled. Since antigens bound by more than one receptor dissociate very slowly, however, such binding may be effectively irreversible on the time scale of lattice formation. Delisi and Perelson (1976) derived several analytical results for the problem, using (implicitly) percolation on Bethe lattices, by a method similar to that of Fisher and Essam (1961) for percolation on the same lattice, including the cluster-size distribution, the behavior of the system near the point at which a large lattice is formed, i.e., near the percolation threshold, and membrane transport. Their results cannot be directly compared with any experimental information, however, as they were obtained with a Bethe lattice, which may not be directly relevant to 2D systems. Their work did indicate definitely the relevance of percolation to such biological phenomena.

18.4 Percolation in Immunological Systems

A relatively old, but still quite useful review of the relationship between immunology and critical phenomena was given by Perelson and Weisbuch (1997). Although immunology is an old research field, the application of percolation theory to it is relatively new, and some of our discussions here may not be upheld by new research in the years to come.^{Footnote 1}

The deadly disease of smallpox has been eradicated from the Earth for many years, but against AIDS (acquired immunodeficiency syndrome) no cure has been found yet, although considerable progress has been made for controlling it and reducing its devastating effects, as of the time of writing this book (2021). The two examples are perhaps the extreme cases success, and the lack of it, in application of immunology. So, how does the immune system work?

If we get the flu in winter, we will not normally get the same thickness shortly thereafter. Our body’s white blood cells or lymphocytes produce antibodies or other cells that are able to neutralize the virus or other foreign antigens. Such antibodies fit the antigens the way a lock fits a key. Apart from small inaccuracies, one type of antigen fits only specific antibodies and vice versa. If we get a new virus, some antibodies accidentally—through mutations—may fit this antigen, our immune system notices this fit, produces more of such specific antibodies, and in this way combats the disease. When we get healthy again, a few very long-lived memory cells survive and allow a quicker response from our body, if the same virus returns. Vaccination produces a controlled amount of sickness so that antibodies and memory cells are formed for the specific disease. The immunodeficiency virus HIV seems to escape destruction by the immune system, and instead destroys slowly the T4 white blood cells, which are crucial as “helper” cells for the functioning of the immune response.

According to the now widely believed ideas of Jerne^{Footnote 2} (1974) antibodies can be treated in turn as antigens by the immune system, which then produces antiantibodies, antiantiantibodies, etc. Moreover, the inaccuracies of the lock-and-key relation between the various types of molecules let the same antibody fit slightly different antigens. In this way, all possible antibody shapes are connected, directly or indirectly, by antibody–antigen relations in a Jerne network, spanning the entire shape space.

What does this have to do with percolation? When we are ill, we take a medicine specific for that illness; not all possible types of medicines have been invented yet. Similarly, our immune system would destroy us if a flu infection could trigger all possible shapes of antibodies. Thus, the immune response should be very specific, i.e., the antibodies triggered by one specific antigen should form a finite cluster in the topological network of possible shapes. They should not percolate throughout this shape space. In other words,

Warning: Immunology has determined that percolation is dangerous to your health!!

de Boer and co-workers (1989, 1992) and Stewart and Varela (1991) investigated the conditions under which the immune response remains limited to a few types of antibodies, once a foreign antigen enters our body, i.e., the conditions under which percolation does not occur. Inherent in percolation theory, as described throughout this book, is the clear distinction between occupied and empty sites. Thus, percolation is the easiest to apply to immunology, if antibody molecules are either there or not there, without a more quantitative and more realistic distinction between such molecules according to their concentration. Indeed, “cellular automata” approximations—discrete systems in which each site can be in one of only two states, sick or healthy, and nothing in between—were developed in immunology by Kaufman et al. (1985) before percolation aspects were investigated. The combination of percolation and cellular automata methods then allowed computer simulation of five- to ten-dimensional shape spaces (Stauffer and Sahimi 1993, 1994; Sahimi and Stauffer 1993), as required for natural immune systems according to Perelson and Oster (1979).

Such studies of localization—not percolating—versus percolation can be done in two ways. We may look at the immune response, if no element of the immune system has yet been triggered. Then, our clusters are simply the sets of activated antibody types connected by the Jerne network bonds (Neumann and Weisbuch 1992). Alternatively, we may look at an immune system that has already evolved into some stationary equilibrium of present and absent antibody types and check for the changes made to this dynamic equilibrium by one specific type of antigen (Stauffer and Weisbuch 1992, Stauffer and Sahimi 1993, Stauffer and Sahimi 1994). This is analogous to the studies of chaos in a dynamical system, and has been used before in genetics (Kauffman 1969) and Ising magnets (Creutz 1986). Physicists often call it damage spreading, whereas experts on pattern recognition and neural networks refer to it as the Hamming distance.^{Footnote 3} Damage is the number of sites that differ in their spin or other characteristic values in a site-by-site comparison of two lattices, if the initial configurations of the two lattices differ by only one site. Hamming distance is used in neural networks to describe the number of pixels that differ in a comparison of two pictures. In particular, if one picture is the original aim and the other picture is an attempt to restore that original picture from a blurred or noisy version of it. Thus, in this sense damage spreading and the Hamming distance are the same.

Various studies have indicated that in the immunological models investigated, the question of whether the immune response remains localized or percolates through the entire system depends on the parameters that one selects, which is hardly a surprise to anyone familiar with the concept of a percolation threshold.

18.5 Percolation Conductivity in Biological Materials

Interactions between water and globular proteins influence their folding, enzymatic activity, and other properties. One way of studying various properties of such biological systems is through measuring their hydration-dependent dielectric losses and electrical conductivity in the hydration range critical for the onset of enzymatic functions. Behi et al. (1982) showed that hydrated protein powders exhibit dielectric dispersion at three distinct frequencies. The first one, $\Delta \epsilon _1$, occurred at frequencies near 1 HZ, with a corresponding change in the DC conductivity of the system, which was attributed to an isotope effect. The second and third dispersions, $\Delta \epsilon _2$ and $\Delta \epsilon _3$, occurred at much higher frequencies, $10^5$ and $10^{10}$ Hz, respectively. Bone et al. (1981) and Behi et al. (1982) proposed that these two dielectric dispersions are due to the Debye relaxation of water bound on the surface of the macromolecule.

Careri et al. (1985) made careful measurements of dielectric losses of lysozyme powders of various hydration levels in the frequency range 10 kHz–10 MHz. The powders were prepared from the Worthington thrice crystalized and salt-free proteins. The dielectric permittivity of the protein is not measured directly, because the sample is only one part of the composite condenser that consists of a layer of the powder included between layers of dry air and glass. In the frequency range of interest to Careri et al., however, the glass did not display appreciable dielectric loss. Thus, their system was essentially a capacitor consisting of two layers of similar thicknesses, one of which had a vanishingly small conductivity and a dielectric constant close to that of the vacuum. The dielectric relaxation time $t_d$ is predicted to be given by

$$\begin{aligned} t_d=\left( \frac{\epsilon _0}{g_e}\right) \sqrt{1+\epsilon }\;, \end{aligned}$$

(18.5)

where $\epsilon _0$ is the vacuum permittivity, $\epsilon $ is the relative dielectric constant, and $g_e$ is the effective conductivity of the hydrated protein. Careri et al. confirmed that their system was insensitive to slight changes in the thickness of air or powder layers. In the frequency range of interest, the hydration dependence of $\epsilon $ was much weaker than that of $g_e$. Thus, a plot of $t_d$ versus h, the hydration level, is essentially equivalent to a plot of $g_e$ versus h.

Figure 18.1 displays the measured conductivity as a function of the hydration level h. To remove from the data the nearly negligible contribution of the non-percolative processes and systematic errors in the evaluation of capacitor geometry to the effective conductivity, the value of $g_e$ at the percolation threshold was subtracted from $g_e$ at all other points for any hydration level h. Careri et al. (1985) established that protonic conduction was the dominating contribution to the dielectric relaxation $\Delta \epsilon _2$ in the frequency range in which the data were collected. The relaxation was attributed to proton displacement on a single macromolecule. Since h is proportional to p, the usual occupation probability in percolation, we can write

$$\begin{aligned} g_e(h)-g_e(h_c)\sim (h-h_c)^\mu \;, \end{aligned}$$

(18.6)

where $h_c$ is the hydration level at the percolation threshold. The critical exponent $\mu $ was evaluated by Careri et al. for three distinct samples, namely, native lysozyme hydrated with H$_2$O, native lysozyme hydrated with D$_2$O, and 1 : 1 complex with (GlcNAC)$_4$ hydrated with H$_2$O. The results are shown in Fig. 18.2. There are two distinct regions: (a) a region for which we obtain a critical exponent $\mu \approx 1.3$ and (b) a region that yields $\mu \approx 2.0$. The first region is interpreted as being an essentially 2D system, as it is relatively thin. In the second region at higher hydration levels, intermolecular water bridges are established, giving rise to a higher dimensional system. Both exponents are in perfect agreement with the values of $\mu $ for 2D and 3D percolation conductivity (see Table 2.3).

What is the mechanism for the observed percolative conduction? Careri et al. (1988) argued that the observed percolation process consisted of proton transfer along a thread of hydrogen-bonded water molecules adsorbed on the protein surface. In their interpretation, water molecules are equivalent to the conducting elements in percolating system. The mean-free path of the protons at or above $h_c$ is the distance between the poles of the macromolecules set by the boundaries of the molecule. The local structure of the protein itself is not important, only the structure of the clusters made of water molecules acting as interconnected conducting elements is significant.

The data have important biological consequences. The protein system used in the studies, with its surface sparsely covered by water or conducting elements, is similar to a protein membrane whose internal surface is sparsely populated, although the membrane itself is immersed in a solvent of near-unit water activity. Percolation theory suggests that membrane conduction is possible with channels only partly filled with conducting elements. Moreover, conduction can be turned off or on by adding or subtracting a few elements, without changing the basic structure of the membrane. Finally, since percolation focuses on the effect of disorder on conduction, we can bypass the need for a high level of structure extending over the full thickness of the membrane or its hydrocarbon core.

Rupley et al. (1988) extended the above studies to lightly hydrated purple membranes, which are more complex than the protein system. The conduction paths may be predominantly within one of several regions of the membrane, such as the lipid surface, the lipid-protein interface, or entirely within the protein. The existence of such preferred paths may give a 2D character to the conduction process in such membranes. Indeed, the measured critical exponent, $\mu \approx 1.23$ was close to that of 2D percolation, although the morphology of the system may be 3D. Bruni et al. (1989) extended such studies to a dry tissue where the conductivity process is integrated into complex, living systems. They used tissues of maize seeds, where water-induced effects lead to the onset integrated metabolism and, thus, to germination. The protonic nature of conduction process in such tissues was established by deuterium substitution. The critical exponent for the conduction was found to be, $\mu \approx 1.23$, close to that of 2D percolation.

18.6 Neuromorphic Computing

Neuromorphic computing is a computing paradigm that utilizes very large-scale integration (VLSI) systems with electronic analog circuits, in order to mimic neuro-biological architectures in the nerve system. More specifically, “neuromorphic” is used to describe analog, digital, and mixed-mode analog/digital VLSI, as well as softwares that implement models of neural systems. Although the research field has at least a 30-year history (see, for example, Mead 1990), it was perhaps fabrication of atomic switches (Terabe et al. 2005) and memristors (Strukov et al. 2008) that led to a new paradigm for constructing architectures that make it possible to carry out neuromorphic computations.

Atomic switches (for a review see Hino et al. 2011) are nanoionic devices that control diffusion of metal cations and their reduction/oxidation reactions in the switching in order to either form or annihilate a metal atomic bridge, which is a conductive path between two electrodes in the on-state. Unlike conventional semiconductor devices, atomic switches provide highly conductive channels, even if their size is on the order of nanometer, which together with their low on-resistance and nonvolatility has made it possible to develop new types of programmable devices.

The existence of memristors, on the other hand, was suggested by Chua (1971) who, based on symmetry considerations, advanced the idea that in addition to the standard resistors R, capacitors C, and inductors L, there must be a fourth fundamental element in electrical circuits. Chua reasoned that six basic relations connect pairs of four fundamental variables in a circuit, namely, the current I, voltage V, charge q, and magnetic flux J. One relation is determined from the definition, $q=\int I(t)dt$, and another from Faraday’s law of induction, hence reducing the number of independent relations to four. Thus, there must also be a fourth circuit element. Chua called the new hypothetical element memory resistor, or memristor for short.

A memristor with memristance M is an element of an electrical circuit characterized by a link between the charge and magnetic flux, $dJ=Mdq$. If M is constant, one recovers the standard linear relation between J and q. If, however, M depends on q, no combination of nonlinear R, C, and L mimics the behavior of such a memristor. Then, a current-controlled memristor for circuit analysis is expressed by

$$\begin{aligned}{} & {} V=\mathcal{R}(w)I\;, \end{aligned}$$

(18.7)

$$\begin{aligned}{} & {} \frac{dw}{dt}=I\;, \end{aligned}$$

(18.8)

where w is the state variable of the device (in this case coinciding with the charge q) and $\mathcal{R}$ a generalized resistance that depends on the system’s internal state. Chua and Kang (1976) considered a general memristive system described by

$$\begin{aligned}{} & {} V=\mathcal{R}(I,w,t)I\;, \end{aligned}$$

(18.9)

$$\begin{aligned}{} & {} \frac{dw}{dt}=f(I,w,t)\;, \end{aligned}$$

(18.10)

and, thus, the flux is no longer defined uniquely by q. Equation (18.9) implies that a memristive system is a particular nonlinear dynamical system, since there would be no current in the memristive system, if the voltage drop across it is zero. Therefore, although a memristor is a passive two-terminal element that is similar to a resistor, it is a particular type of nonlinear resistor that provides controllable resistance. As an electrical device, it is not yet available commercially.

Chua’s proposal was not tested experimentally until Strukov et al. (2008) showed how memristive behavior arises in physical systems. They considered a thin semiconducting film of thickness $\ell $, sandwiched between two metal contacts and divided into two parts: one with a high concentration of dopants (positive ions) and a low resistance $\mathcal{R}_\textrm{ON}$, in series with an undoped material with a high resistance $\mathcal{R}_\textrm{OFF}$. If an external voltage V(t) is applied to the film, due to the dopants’ drift the boundary between the doped and undoped parts moves forward. Thus, even in the simplest case—an ohmic conductor (OC) and linear ionic drift with average ionic mobility $\mu $—one obtains

$$\begin{aligned} V(t)=\left[ \mathcal{R}_\textrm{ON}\frac{w(t)}{\ell }+\mathcal{R}_\textrm{OFF}\left( 1- \frac{w(t)}{\ell }\right) \right] I(t)\;, \end{aligned}$$

(18.11)

$$\begin{aligned} \frac{dw(t)}{dt}=\mu \frac{\mathcal{R}_\textrm{ON}}{\ell }I(t)\;, \end{aligned}$$

(18.12)

$$\begin{aligned} w(t)=\mu \frac{\mathcal{R}_\textrm{ON}}{\ell }q(t)\;. \end{aligned}$$

(18.13)

In this case, w(t) represents the doped area’s width. Combining Eqs. (18.11) and (18.13) yields an expression for the system’s memristance M(t). Since, $\mathcal{R}_\textrm{ON}\ll \mathcal{R}_\textrm{OFF}$, one obtains

$$\begin{aligned} M(t)=\mathcal{R}_\textrm{OFF}\left[ 1-\frac{\mu \mathcal{R}_\textrm{ON}}{\ell ^2}q(t) \right] \;. \end{aligned}$$

(18.14)

The second term in the brackets is responsible for the memristive effect, since it is proportional to $\mu $ and $1/l^2$, which is very important at nanoscale.

Just as an Ohmic heterogeneous material with constant linear conductivity is too complex to be represented by resistors in series and/or parallel configurations, so also is a memristive material and, therefore, cannot be represented as memristors in parallel or series. Random resistor networks and their generalizations, described throughout this book, were developed precisely for this reason. Likewise, a network of interconnected memristive elements, or one that contains a mixture of memristors and OCs, is more realistic and versatile, and was fabricated by Borghetti et al. (2009) to which transistors were connected, and was capable of self-programming. The first neuromorphic computations involving a network of interconnected memristors were carried out by Nedaaee Oskoee and Sahimi (2011). In practice, neuromorphic computing on the hardware level can be realized by oxide-based memristors, spintronic memories, threshold switches, and transistors.

Using 2D models, simulations, and experiments, Brown and co-workers (Fostner et al. 2014; Fostner and Brown 2015; Mallinson et al. 2019; Pike et al. 2020) showed that one can obtain neuromorphic behavior in thin percolating films of nanoparticles. Experimentally, Pike et al. (2020) formed percolating networks of nanoparticles by depositing conducting nanoparticles onto silicon nitride substrates, which brought the particles into contact and formed interconnected clusters, made of particles in Ohmic contact, which are separated by tunnel gaps that have a distribution of sizes with fractal geometries. The tunnel gaps act as switching sites such that upon application of an external voltage stimulus (see below), atomic-scale filaments can be formed, and subsequently broken, in the tunnel gaps, changing the network conductance G. Deposition was terminated when the fraction p of the surface area covered with conducting particles approached the percolation threshold $p_c$. After deposition the overall structure of the network was held fixed.

To simulate the same system, Fostner and Brown (2015) used a model consisting of overlapping disks with $200\times 200$ particle diameters, which resembled a continuum percolation model of the type described in Chap. 3 with a percolation threshold, $p_c\approx 0.676$. Tunneling was allowed explicitly to occur in the system. It was shown (Ambrosetti et al. 2010) that, with tunneling present, the effective conductivity of the system below $p_c$ depends exponentially on surface coverage p. Simulating the system below $p_c$, where the conductance G of the system is due to tunneling across gaps between groups of connected particles, Fostner and Brown (2015) assumed that the connections between overlapping particles to have negligible resistance, and that each gap has a conductance, $G_g=A\exp (-\beta L)$, with A and $\beta $ being two constants, and L the size of the gap, measured in units of the disks’ diameter.

To study the properties of the network, the voltage across the system was increased in a triangular ramp from $V=0$ to $V_\textrm{max}$ and back to $V=0$, with a constant step size of 0.025 V, and a number $N_V$ of voltage steps was used to characterize both the voltage at any time and the time since the start of the ramp. Voltage ramping was done for consistency with the experimental protocol used in Ref. Fostner et al. (2014). Thus, all the quantities of interest may be studied as a function of $N_V$.

The switching process was simulated by two models, one stochastic, and one deterministic. First, the smallest gaps between each pair of clusters were identified, and the gaps with electric fields larger than a preset threshold $E_t$ were replaced with a large conductance, $G_O=10\;\Omega ^{-1}$ with probability $p_s$, which accounted for the stochastic nature of the switching process. Pike et al. (2020) refined the model by allowing for two switching probabilities, $p_s^+$ and $p_s^-$, such that when the electric field in a gap or the current in a filament was greater than $E_t$, the switch was allowed to change state—switch on with probability $p_s^+$, or off with probability $p_s^-$. Doing so simulated the formation in the tunnel gap of an atomic wire, generated by surface diffusion or evaporation induced by the electric field. The effective conductance G of the network was then recalculated. The process, as well as the voltage ramps, continued until the system conductance converged to within 0.1% of the value on the previous step. For all the simulations, the cumulative number $N_R$ of the replacements, or the switching events, as well as G and the electric current I, were recoded.

Figure 18.3 presents the macroscopic current for disk coverages $p=0.55$ and 0.65, and the corresponding voltage ramp V, versus the number of voltage steps $N_V$. They indicate that, for $p=0.55<p_c$, the initial conductance is small, because the current flows through small tunneling conductors. It then increases in response to the increasing applied voltage during the second voltage cycle, because the tunneling conductors are replaced (switched) by large conductances that represent atomic-scale wires. The switching event also reduces the potential difference between the two clusters of nanoparticles. Further increase in the current is caused by the subsequent cycles of the voltage ramp. After the fifth cycle, however, the network conductance is dominated by the high conductances $G_O$ and saturates.

In contrast, there is large initial current for $p=0.65$, and the conductance G saturates much more rapidly, since the tunneling gaps are both fewer and smaller and, therefore, the applied voltage causes a stronger increase in the electric field, increasing the number of switching events. Note that if the switching probability $p_s$ is small, a larger number of voltage cycles would be needed, but eventually the current saturates again. The same type of patterns were observed in the experiments (Pike et al. 2020). These processes are qualitatively similar to leaky integration and fire mechanisms in biological neurons (see, for example, Burkitt 2006).

Figure 18.4 presents the dynamic evolution of $N_R$ and G with $N_V$ (time) for a range of $V_\textrm{max}$ and $p_s$, indicating that increasing $p_s^+$ leads to a faster increase in the number of replacements (switching) and, thus, in the conductance. But, it also increases the fluctuations in G. Although, as expected, the maximum conductance for each set of parameters is the highest for the largest $V_\textrm{max}$, the effect of $V_\textrm{max}$ is more subtle. To see this, consider Fig. 18.4a and d, for which the threshold $E_t$ was exceeded by only a relatively small number of tunnel gaps. Thus, the number of replacements was not large and the conductance saturated at small values, but still in the tunneling regime. For large $V_\textrm{max}$, once $E_t$ was exceeded, both $N_R$ and G increased nearly exponentially until the last tunneling conductance on the conduction path between the contacts was replaced (marked by colored dots), hence defining a critical number of voltage steps. $N_R$ eventually saturates because most of the relevant tunnel gaps are replaced, and the existence of an Ohmic conductance path implies that the voltage distribution is relatively uniform, so that the number of places that can exceed the threshold $F_t$ is very few.

Some of such features were qualitatively similar to what happens in biological systems. Moreover, potentiation of the film happens for all surface coverages p, where potentiation refers to the sensitization of a conduction pathway that corresponds to a memory response in the system of interest, and results from an avalanche of neuronal connections. The size of the neuromorphic effect—the change in the conductance G—is more striking for systems with low particle coverages, because the films are more complex. As $p\rightarrow p_c$, the size of the connected groups of particles increases and approaches that of the simulated system, implying that the number of tunneling connections that must be traversed in order to span the system becomes smaller. This implies that the film contains fewer locations at which atomic-scale wires could be formed and, hence, that potentiation involves fewer switching events.

To obtain closer agreement with the experiments, Pike et al. (2020) developed a second model. They replaced the stochastic switching with a deterministic one in which the size $d_i$ of each tunnel gap was changed in response to the electric field $E_i$ in the gap according to

$$\begin{aligned} \Delta d_i=\left\{ \begin{array}{cc} r_d(E_i-E_t)\;, &{} \textrm{if}\; E_i\ge E_t\\ 0\;, &{} \textrm{otherwise} \end{array}\right. \end{aligned}$$

(18.15)

The current flow $I_j$ in each existing filament causes electromigration effects that decrease its width $w_j$ according to

$$\begin{aligned} \Delta w_j=\left\{ \begin{array}{cc} r_w(I_j-I_t)\;, &{} \textrm{if}\; I_j\ge I_t\\ 0\;, &{} \textrm{otherwise} \end{array}\right. \end{aligned}$$

(18.16)

Here, $r_d$ and $r_w$ are parameters that control the rates at which d and w change when, respectively, threshold fields $E_t$ and currents $I_t$ are exceeded. Experimentally, $E_t=10$ V and $I_t=0.01$ A, which were also used in the simulations.

The simulation results with the deterministic model and their comparison with the experimental data are shown in Fig. 18.5. There is a qualitative agreement between the data and the simulation results. Most important is the fact that the switching produces avalanches in the system as a result of the change in the macroscopic conductance G, which are qualitatively similar to observed neuronal avalanches in cortical tissue (see, for example, Fontenele et al. 2019). Both experiments and simulations indicated that the probability density functions of the size s and duration $T_a$ of the avalanches follow power laws (Mallinson et al. 2019)

$$\begin{aligned} f(s)\sim & {} s^{-\tau }\;, \end{aligned}$$

(18.17)

$$\begin{aligned} h(T_a)\sim & {} T_a^{-\alpha }\;, \end{aligned}$$

(18.18)

both with a long tail at large times. Note the similarity between Eq. (18.17) and the percolation cluster-size distribution, described in Chap. 2. Equations (18.17) and (18.18) are the signature of a scale-free critical state. Experiments yielded, $\tau \simeq 1.6\pm 0.1$ and $\alpha \simeq 2.0\pm 0.1$. The deterministic model described above produced, $\tau \simeq 2.1\pm 0.1$, and $\alpha \simeq 2.6\pm 0.1$, in qualitative agreement with experiment. In addition, the autocorrelation function also followed a power law

The similarities and differences with a network of interconnected memristors should also be pointed out. The tunnel junctions, formed at gaps in the percolating film, which switch to a highly conducting state on the formation of an atomic scale wire, behave differently than memristors. But the tunnel gaps are memristive in the sense that, similar to a memristor, it is the device history that determines its state. Switching occurs at a well-defined threshold value of the electric field, but the memristors in the junction cannot be switched back to the low conduction state by simply reversing the polarity of the bias voltage.

18.7 Sensory Transmission and Loss of Consciousness

How dynamic communications between neurons lead to consciousness, or its loss, is a fundamental question that has been studied for decades. In addition to clinical studies, many attempts have been made to model individual electroencephalographic (EEG) features that are associated with loss of consciousness during general anesthesia (see, for example, Ching et al. 2012; Wang et al. 2014, and references therein), including the development of a stochastic model to describe general anesthesia as a thermodynamic phase transition (Steyn-Ross et al. 2001).

Zhou et al. (2015) developed a percolation model to determine the flow of information between neurons. The model produced typical EEG features under general anesthesia, as well as dose–response characteristics for loss of consciousness. Their model was based on a layered hierarchical fractal structure, from an input node to multiple output ones. The layers, meant to represent laminar and divergent organization seen (Sherman 2012) in mammalian thalamocortical structures, were generated by scale-invariant fractal expansion. Within each layer small-world (SW) network (see Chap. 17) properties among the nodes were generated using the Watts and Strogatz (1998) algorithm. As described in Chap. 17, an SW network is a graph in which most nodes are not neighbors of one another, but any node’s neighbors are likely to be neighbors of each other, and most nodes can be reached from every other node by a small number of steps. Then, the typical distance L between two randomly selected nodes scales as, $L\propto \ln N$, where N is the number of nodes in the network. Such networks have found a large number of diverse applications, from geoscience (Yang 2001) to brain (Sporns et al. 2004; see also Chap. 17).

A directional weight $W_{ij}$ from node j to node i was attributed to each edge of the network, with $W_{ij}\ne W_{ji}$ to represent the counterstream architecture of human brain. The anterior nodes were differentiated from the posterior ones by attributing different edge weights in the feedforward and feedback directions. Suppose that $A_i(t)$ and $P_j(t)$ are, respectively, the neural activity of node i and the preceding input from node j at time t. $A_i(t)$ is simply the weighted average of activities from all input nodes,

$$\begin{aligned} A_i(t)=\frac{\sum _iW_{ij}P_j(t)}{\sum _jW_{ij}}=\frac{W_{ii}P_i(t)+\sum _{j\ne i}W_{ij}P_j(t)}{W_{ii}+\sum _{j\ne j}W_{ij}} \end{aligned}$$

(18.19)

where the input function $P_j(t)$ represents the accumulated history of neural activity from the preceding m time steps, weighted by exponentially decaying memory, i.e.,

$$\begin{aligned} P_j(t)=\frac{\sum _{\tau =1}^m\exp (-\tau )A_j(t-\tau )}{\sum _{\tau =1}^m\exp (-\tau )} \;, \end{aligned}$$

(18.20)

Percolation is then used to assign randomly the weights $W_{ij}$ to the edges from a sampling function, with probability p that represented the likelihood of activity transmission, according to

$$\begin{aligned} W_{ij}=\left\{ \begin{array}{cc} \mathcal{G}_\textrm{CDF}(U[0,1]) &{} i\ne j\\ c\exp (-\lambda \sum _{k\ne i}W_{ki}) &{} i=j \end{array}\right. \end{aligned}$$

(18.21)

Here, $\mathcal{G}_\textrm{CDF}$ is the Gaussian cumulative distribution function with its center at $(1-p)$ and a standard deviation of 0.05, U[0, 1] represents the uniform distribution, and c and $\lambda $ are constant.

With decreasing p the probability of activity transmission along individual edges also decreases, representing the inhibition of information flow under anesthesia. Although various anesthetic classes act differently at molecular and cellular levels, the net effect is represented as a global inhibition of arousal (Brown et al. 2011). Sampling $W_{ij}$ is done independently by p, when $i\ne j$, whereas for $j=j$ $W_{ii}$ represents the memory of the past activity of the same node, which is influenced by the incoming connection strength. When all non-self-connections diminish, i.e., when $\sum _{k\ne i} W_{ki}\rightarrow 0$, self-connection is the dominating factor. Neuronal transmission, including axonal propagation and synaptic events, involves cycles of receptor inactivation, activation, and deactivation. To account for such dynamic behavior, edge weights were periodically resampled using the above sampling algorithm, with a periodicity proportional to $\exp (-ap)$, where a is a constant.

The percolation model was then validated by reproducing four key features of global EEG responses: (a) characteristic EEG waveforms, including burst suppression under deep anesthesia. A burst–suppression pattern is a discontinuous EEG that has periods of marked suppression or isoelectric intervals, alternating with “bursts” of activity, with or without embedded epileptiform features, and is commonly seen following hypoxic-ischemic injury. (b) An EEG power shift to lower frequencies with increasing anesthetic concentrations. (c) Synchronization of cortical nodes and (d) shift of $\alpha $ and $\delta $ rithms to the anterior, often called anteriorization. An $\alpha $ rithm is a uniform rhythm of waves in the normal EEG, exhibiting an average frequency of 10 per second, which is typical of a normal person awake in a quiet resting state. A $\delta $ rithm is an EEG wave with a frequency less than 3.5 per second, which is typical in deep sleep, in infancy, and in serious brain disorders.

As documented by Zhou et al. (2015), their percolation model reproduced all the four key features. The model indicated that as p decreases, a steep divergence of the transmitted signal from the original is developed, together with a loss of signal synchrony, as well as a sharp increase in information entropy in a critical manner, which resembles the precipitous loss of consciousness during anesthesia.

18.8 Actomyosin Networks

Actomyosin is the actin–myosin complex that forms within the cytoskeleton, and is inherently contractile, with the myosin motor protein able to pull on actin filaments, giving rise to contractile fibers that form the basis of skeletal muscle, and enabling cell motility and force generation at the sub-cellular level. An actomyosin network consists of self-assembled actin filaments—two-stranded helical polymers of the protein actin, with a diameter of around 6 nm—which have variable lengths, and are connected by cross-linkers and myosins.

Alvarado et al. (2013) (see also Sect. 18.15 below) fabricated a quasi-2D network of actin filaments with the size $2.5\times 2.5$ mm$^2$, which were connected by fascin cross-linkers. Myosin motors pull the filaments together and exert contractile forces on the network, leading to cross-linker unbinding and actin filament movement that fundamentally altered the underlying actin network. Alvarado et al. traced the trajectories of the actin filaments during the contraction, in order to reconstruct the connectivity of the actin network throughout the dynamical process. Assuming that the number of cross-linkers in the system was $M_c$, three distinct regimes were observed: (a) For small $M_c$, the actin network was contracted into many foci, when the myosins were activated. Thus, if one retraces the actin molecules’ position back to time $t=0$, one observes that the resulting in foci originated from actin clusters with areas of similar sizes. (b) For intermediate values of $M_c$, the retracing to the original areas indicated a power-law distribution in sizes with an exponent of −1.91, which was similar to the exponent $\tau $ of the cluster-size distribution at the percolation threshold (see Chap. 2), $\tau =187/91\simeq 2.05$. The power-law cluster size distribution existed for a wide range of cross-linker concentrations, hence indicating that a critical state can be reached without fine-tuning, which is the hallmark of self-organized criticality. (c) For large $M_c$, the actin network contracted to a single cluster.

Lee and Pruessner (2016) proposed a 2D percolation model for the experimental observations that reproduced their essential features. They used a $L\times L$ square network with $N=L^2$ nodes. A strength $s_k$ was attributed to a site k, selected from the Poisson’s distribution with a mean value, $\langle s \rangle =M_c/N$, intended to mimic the random number of cross-linkers connected at each node. It was then assumed that at $t=0$ the pulling force $f_k$ at site k, which are generated by the molecular motor myosins, were also distributed according to Poisson’s distribution with a mean value, $\langle f\rangle = M_m/N$, with $M_m$ being the number of myosin motors in the system. In addition to $f_k$, adhesion connects the actin network to the boundary, and to mimic that Lee and Pruessner assumed that the strength of the boundary adhesion was $s_A$, and that the actin networks adhered to the top and bottom boundaries only.

Next, one must introduce the myosin-induced dynamics in the actin network. Lee and Pruessner assumed that when the pulling force $f_k$ begins to act on nodes k, the sites rip and, therefore, disappear sequentially according to the reverse order of their unbinding propensity function $B(s_k,f_k)$, with B decreasing with $s_k$, but increasing with $f_k$. The disappearance of a node simulates unbinding of the cross-linkers from the actin filament due to loading. When a node disappeared, the pulling forces were adjusted to achieve force balance in the network, which was tantamount to assuming that the time scale $t_u$ for unbinding was much larger than the corresponding time $t_b$ to reach force balance.

In addition to the myosin unbindings, connected clusters of actin filaments may also contract. To include this in the model, Lee and Pruessner assumed that they occur if a connected cluster is detached from either boundary, and that the contraction time scale $t_c$ is small compared to $t_u$, so that when a cluster is free to contract, it contracts to a focal point before the next unbinding event occurs. When a cluster was detached from the boundary, it could freely contract, and all forces due to mechanical tension ceased. Thus, the contracted cluster was removed from dynamics of the network and, thus, no node was deleted from the cluster.

It should be clear to the reader that the model resembles an invasion percolation (IP) process with trapping. The sequential node disappearance according to the reverse order of their unbinding propensity is an IP algorithm. That no node was removed from a detached cluster implies that such nods were trapped and could not be affected by the dynamics of the network.

The percolation model could also explain the aforementioned three regimes in the experiments of Alvarado et al. (2013). If $M_c$ is small, many nodes will be empty and, thus, the actin network could not percolate, and fragments into many small foci. If, on the other hand, $M_c$ is large, then the average strength of a node may be larger than the adhesion strength, $\langle s\rangle >s_A$, resulting in the detachment of the entire network from the adhesion boundary due to internal contraction. In between, one has networks with intermediate values of $M_c$, which is the regime in which the actin network percolates and the cross-linker unbinding occurs predominately in the interior.

The simulation of Lee and Pruessner confirmed the picture. However, the exponent $\tau $ of the cluster-size distribution was determined to be very close to $\tau =187/91\simeq 2.05$, indicating that the trapping effect was minimal and did not influence the results significantly. As a result, the model is close to the standard percolation model.

Although they did not mention percolation explicitly, a model that was developed by Wang and Wolynes (2012) for active contractility in actomyosin networks represents a percolation model in high-dimensional space, or on random graphs described in Chap. 17. Their model accounted for two key features of actomyosin self-organization, namely, the asymmetric load response of individual actin filaments, and the correlated motor-driven events that mimicked myosin-induced filament sliding. Wang and Wolynes (2012) studied the effect of the concentration and susceptibility of the motors on their collective behavior, as well as their interplay with the network connectivity, and how they regulate macroscopic contractility. Their work demonstrated that cooperative action of load-resisting motors in a force-percolating structure integrates local contraction/buckling events into a global contractile state through an active coarsening process.

18.9 Percolation Transition in the Mutation-Propagating Module Cluster in Colorectal Tumorigenesis

Tumorigenesis, the initial formation of a tumor in the body, involves epistatic interactions between genes and multiple somatic mutations. Certain types of genes do not mutate randomly, but tend to co-occur in cancer patients. In particular, colorectal cancer develops through the sequential accumulation of driver mutations, including adenomatous polyposis coli (APC), Kirsten rat sarcoma viral oncogene homolog (KRAS; pronounced K-raz), phosphatidylinositol 3-kinase (PI3K), and tumor protein p53 (TP53), hence suggesting that an important factor for the development and spread of cancer may be the cooperativity of driver mutations. The mechanism underlying such cooperative multiple mutations during tumorigenesis is not well understood, however.

Shin et al. (2017) proposed a percolation model in order to understand the cooperative phenomenon of somatic mutations in tumorigenesis in colorectal cancer. They first obtained somatic mutation data for colorectal cancer patients. Then, in order to use only high-confidence pathogenic variants, Shin et al. filtered the mutations based on their predicted pathogenicity, derived from five functional mutation prediction tools, which were then assigned to tumor samples by abstracting binary event calls, such that a genomic event either occurred—represented by “1”—or did not occur - denoted by “0”—in a gene for a given sample. Only patients with less than 300 mutations from the 223 colorectal cancer patients in The Cancer Genome Atlas (TCGA) were used, reducing the total number of patients to 198, for whom the expression data were obtained from the Firehose website.

The protein–protein interaction (PPI) network had been previously constructed by Hofree et al. (2013), of which Shin et al. considered only the largest connected subnetwork, because the network propagation was not available between a pair of nodes that did not connect with each other, either directly or indirectly. Shin et al. then integrated the expression and mutation profiles of the cancer patients with the PPI network, in order to obtain the adjacency matrix A of the network, as well as a patient-by-gene matrix that displayed the mutation profiles of the aforementioned binary (1, 0) states on $N=10,968$ genes for the patients.

The entry $A_{ij}$ of A represents the probability of interaction between nodes i and j along which the mutation influence propagates, which, for a given sample, was assumed to be proportional to the product of expression values of both genes in that sample, i.e., $A_{ij}\rightarrow A_{ij}E_iE_j$, with $E_i$ indicating the expression value of gene i. Thus, patient-specific PPI networks were obtained, enabling the implementation of more realistic propagation of somatic mutations on the network (see below). Shin et al. used only samples that had both RNA-seq and somatic mutation data, obtained from the Dana-Farber Cancer Institute, and downloaded from the cBioPortal website.

To simulate the propagation of mutation effects in the PPI network, Shin et al. utilized an algorithm based on a random walk with restart on a network, developed by Vanunu et al. (2010). When a mutation occurred on a gene in the PPI network, a value of 1 was initially attributed to the mutated gene, which then propagated along the network neighborhood, such that higher values were assigned to non-mutated genes that were closer to the mutated one, according to, $\textbf{F}_{t+1}=\alpha \textbf{A}'{} \textbf{F}_t+(1-\alpha )\textbf{F}_0$ at step or time t, with $\textbf{F}_0$ being a matrix of binary (1,0) states on genes such that 1 or 0 indicated whether a gene was mutated or not in the corresponding patient; $\textbf{A}'$ a degree-normalized adjacency matrix of the PPI network, and $\alpha $ is the degree of diffusion of a mutation influence in the network. $\textbf{A}'$ was set to be, $\textbf{A}'=\textbf{D}^{-1/2}{} \textbf{A}\textbf{D}^{1/2}$, in which D is a diagonal matrix such that its entries $D_{ii}$ were the sum of row i of A

A mutation cannot successfully affect all the nodes in the PPI network, but covers at most a few layers of nearest neighbors of the nodes. In this way, the effective boundary of mutation influences centered on the mutated node is predicted, inside which nodes have influence scores above a given threshold, and form a sub-network called a mutation-propagating module (MPM). The MPMs of individual mutations are distributed throughout the network, forming connected modules, the largest one of which, referred to as the giant cluster (GC; see Chap. 17)—the analog of the largest percolation cluster on lattices—could cover the network on the order of network size N. The GC can, therefore, impact significantly the properties of the network.

The threshold of mutation influences was set high enough that each MPM approximately included only nearest neighbors. Simulation of Shin et al. indicated that, on average, the GC covered about 20–40% of the entire network, even though patients only had an average of 20–40 somatic mutations, hence indicating that multiple somatic mutations interact cooperatively to expand the effective boundary of mutation influences. To check whether the cooperative effects in the GC in the network for each individual patient was due to nonrandom mutation profile, Shin et al. compared the size $S_\textrm{GC}$ of the GC with the expected size $S_\textrm{GC}^r$ of the same number of mutations that occurred randomly on the network. The former turned out to be much larger than the latter, with the differences being more striking when the cancer-related or driver genes were considered, rather than all the mutations. Thus, the cooperative effect among somatic mutations in cancer cannot be due to a random selection of mutations, but is determined by topological— connectivity—properties of somatic mutations in the PPI network. That is, of course, the signature of the percolation process.

Other aspects were examined as well. For example, if a GC forms based on nonrandom mutation profiles of a cancer patient, genes within it must be correlated with relevant biological and clinical characteristics of the patient. To test this and identify the biological processes that are inherent in the cluster, Shin et al. (2017) used a set of 50 “hallmark genes,” and and identified hallmark gene sets enriched in the GC. The hallmark gene sets summarize and represent specific well-defined biological states or processes, display coherent expression, and were originally generated by a computational methodology based on identifying gene set overlaps and retaining those that display coordinate expression. They reduce noise and redundancy, providing a better-delineated biological space for GSEA (gene set enrichment analysis). It was hypothesized that the cooperation of multiple somatic mutations in a GC activates multiple cancer-related hallmarks that cooperatively enhance tumorigenesis.

Two tests were carried out to check the hypothesis. First, Shin et al. examined how closely the GC-based patient classification matches the previous consensus molecular subtypes (CMS) classification, proposed by Guinney (2015), which is considered a robust molecular classification for colorectal cancer. They also investigated which key biological features the subtypes have. The result was that the analysis of the hallmark gene set based on the giant percolating cluster (GPC) extracts all the features of the CMS groups, even though factor analysis was used, which is a type of unsupervised learning for feature selection or for data reduction. Moreover, the patient classification with the selected features exhibited very strong correlations with the CMS groups.

Second, Shin et al. showed that the hallmark gene sets enriched in cancer patients correlate with clinical tumor stages. To do so, they carried out statistical clustering analysis of the enrichment score of the hallmark gene sets for cancer patients according to the significant hallmark gene sets, of which 12 were selected to distinguish the tumor stages of cancer patients. The result indicated that the patient population may be divided into five clusters, of which three are strongly correlated with various tumor stages.

According to the analysis of Shin et al., scattered MPMs formed a GC in a PPI network and, more importantly, genes within the GC represented the phenotypic properties of cancer, including cell proliferation and metastasis. Therefore, Shin et al. proposed that the GC that results from somatic mutations of a cancer patient should be viewed as a GPC, i.e., the largest percolating cluster that integrates the influences of scattered somatic mutations, so that it confers phenotypic changes corresponding to cancer hallmarks. The question then is whether in practice the GC would also undergo a percolation transition on the accumulation of somatic mutations during cancer progression. The analysis of Shin et al. indicated that the answer is affirmative.

Thus, a GPC—cooperative mutation effects represented by a large connected cluster in a PPI network—can undergo a percolation transition during tumorigenesis. In the development of normal tissue, a number of somatic mutations can occur at various places in the PPI network, as a result of which the MPMs form a connected module, or a GC along with the accumulation of somatic mutations. Thus, tumorigenesis can be initiated by a certain driver mutation that connects scattered clusters into one, leading to the formation of a GPC that represents cancer hallmark. This is shown schematically in Fig. 18.6. Shin et al. also found that the most frequently observed sequence of driver mutations that characterize colorectal cancer development might have been optimized in order to maximize the size of the GPC, hence providing a intriguing insight into the relationship between cancer development and percolation transition, which can be useful for understanding the fundamental mechanism of tumorigenesis.

18.10 Essential Nodes in Brain: Optimal Percolation

A most important problem in neuroscience is how the brain integrates distributed and specialized networks into a coherent information processing system. It is generally believed (see, for example, Sporns 2013) that integrating specialized modules in the brain are facilitated by a set of what is referred to as the essential nodes. Integrated brain networks exhibit long-range correlations in their activity, with the correlation structure often referred to as functional connectivity. Any perturbation in the essential nodes results in large disturbances in functional connectivity that affect global integration. Therefore, a fundamental step toward understanding how information is processed in brain circuits entails identifying the essential nodes, which in turn may contribute to the design of targeted interventions to restore or compensate dysfunctional correlation patterns in disease states of the brain.

There have been many studies utilizing various network centrality measures in order to identify the essential nodes. Such measures include hubs, nodes that have many connections to other nodes; betweeness centrality, a measure of centrality in a connected graph since there exists at least one shortest path between every pair of nodes; closeness centrality, a measure of centrality of a node in a network and calculated as the sum of the length of the shortest paths between the node and all other nodes in the graph; eigenvector centrality that measures the influence of a node in a network by assigning relative scores to all nodes in the network based on the fact that connections to high-scoring nodes contribute more to the score of a given node than equal connections to low-scoring ones, and other measures.

Each of such measures may be used as a way of ranking the nodes in order to determine the most influential nodes, with the nodes having the highest ranking viewed as essential for integration in the brain. Although each measure provides a distinct aspect of influence, they all predict that when the essential nodes are inactivated in a targeted intervention, integration in the overall network is to a large extent prevented. In other words, inactivated nodes with the highest rank inflict the largest damage to the long-range correlations. Therefore, the optimal centrality measure would be one that prevents integration of the network by inactivating the fewest number of nodes.

Del Ferraro et al. (2018) suggested that the problem can be addressed using percolation theory. One must identify the smallest set of nodes that, when inactivated, destroy the integration of the network, which is an NP-hard (non-deterministic polynomial-time hardness) problem. Del Ferraro et al. suggested that the problem can be mapped onto the optimal percolation (OP), first described by Morone and Makse (2015). To understand the OP, first recall that the most influential nodes in a network are those that form a minimal set that guarantees global connection of the network, and are referred to as the optimal influencers. To be concrete and to see the relation with the percolation problem, consider a network of N nodes that are connected by M bonds or links, with an arbitrary degree distribution (i.e., arbitrary distribution of the nodal connectivity). As the reader knows, if the nodes are removed at random, then at a critical fraction $q_c$ of the removed nodes—the percolation threshold—the network undergoes the percolation transition and its structure collapses. Therefore, the optimal influence problem is equivalent to estimating the minimum $q_c$ of the influencers to fragment the network, so that the probability of having the GC would vanish.

This was the problem that was addressed by Morone and Makse (2015). They reformulated the problem as follows so that it can be addressed in a new framework. Consider a vector $\textbf{n}=(n_1,\cdots ,n_N)$ in a network of N nodes, such that $n_i=0$ for the removed influencers, and $n_i=1$, otherwise, with the fraction q of the removed influencers given by, $q=1-(N\sum _i n_i )^{-1}$. Suppose that $p_{i\rightarrow j}$ is the probability that a node i belongs to the GC in a network in which j has been removed. If there is no GC, then $p_{i\rightarrow j}=0$ for all $i\rightarrow j$. We construct a $2M\times 2M$ matrix, defined by

$$\begin{aligned} \left. \mathcal{M}_{k\rightarrow \ell ,i\rightarrow j}\equiv \frac{\partial p_{i\rightarrow j}}{\partial p_{k\rightarrow \ell }}\right| _{\{p_{i\rightarrow j}=0\}}\;. \end{aligned}$$

Morone and Makse (2015) showed that for locally tree-like random graphs one has, $\mathcal{M}_{k\rightarrow \ell ,i\rightarrow j}=n_i\mathcal{B}_{k\rightarrow \ell ,i\rightarrow j}$. Here, $\mathcal{B}_{k\rightarrow \ell ,i\rightarrow j}$ is the non-backtracking matrix of the network, a binary matricial representation of the topology of a network whose entries represent the presence of non-backtracking paths between pairs of different nodes, traversing a third intermediate one (Pastor-Satorras and Castellano, 2020). The entries of the matrix are non-zero, taking up a value of 1, only when $(k\rightarrow \ell ,i\rightarrow j)$ form a pair of consecutive non-backtracking directed edges; i.e., $(k\rightarrow \ell ,\ell \rightarrow j)$ with $k\ne j$.

One obtains a stable solution for the problem $\{p_{i\rightarrow j}=0\}$, if the eigenvalues of $\mathbf{\mathcal{M}}$, $\lambda (\textbf{n};q)\le 1$. Thus, for a given $q\ge q_c$, the optimal influence problem may be reformulated as one of identifying the optimal configuration n that minimizes the largest eigenvalue $\lambda (\textbf{n};q)$ such that the optimal set $\textbf{n}^*$ of $Nq_c$ influencers is obtained when the minimum of the largest eigenvalue is at its critical threshold,

$$\begin{aligned} \lambda (\textbf{n}^*;q_c)=1\;. \end{aligned}$$

(18.22)

Morone and Makse then devised a numerical scheme to determine the configuration for which the largest eigenvalue approaches one. To do so, they recast the problem as one of optimization, whereby one minimizes the “energy” (or the cost function) of a many-body system (i.e., many interacting nodes in the network), where the form of the interactions was fixed by the non-backtracking matrix of the network, defined above, and minimizing the energy is tantamount to determining the configuration $\textbf{n}^*$ with minimum $q_c$ at which the GC collapses. To do so, they defined a ball of radius $\ell $ around every node, considered the nodes that belonged to the frontier $\partial $Ball$(i,\ell )$ and assigned to i a collective influence (CI) strength at level $\ell $, defined by

$$\begin{aligned} \textrm{CI}_\ell (i)=(k_i-1)\sum _{j\in \partial \textrm{Ball}(i,\ell )}(k_j-1)\;, \end{aligned}$$

(18.23)

where $k_i$ is the degree of node i.

Thus, the computational algorithm based on Eq. (18.23) is as follows. Initially, all the nodes are present and, therefore, $n_i=1$ for all of them. One node $i^*$ with the largest CI$_\ell $ is removed, $n_{i^*}$ is set to zero, and the degree of neighbors of $i^*$ is reduced by one. The procedure is repeated to identify a new node with the largest CI$_\ell $ for removal. The computation is terminated when the GC does not exist anymore. The algorithm was implemented with three large networks, namely, the Erdös–Rényi graph (see Chap. 17), the web of Twitter users, and an immunization scheme on a personal contact network built by the phone calls by 14 million people in Mexico. In the last case, the method saved a large number of vaccines or, equivalently, identified the smallest possible set of people to quarantine.

Del Ferraro et al. (2018) used the OP algorithm of Morone and Makse, together with pharmacogenetic interventions in vivo, to predict and target nodes that are essential for global integration of a memory network in rodents. The computations predicted that integration in the memory network is mediated by a set of low-degree (low connectivity) nodes that are in the nucleus accumbens, which was confirmed with pharmacogenetic inactivation of the nucleus accumbens that eliminates the formation of the memory network, while inactivations of other brain areas leave the network intact.

18.11 Flexibility of Thought in High Creative Individuals

Flexibility of thought is considered a hallmark of creativity, because it indicates how an individual can process new data and ideas and expand his/her thought process. But, how can one quantify the flexibility of thought in a person? One approach is based on network science that represents semantic memory structure as a network, and was proposed by Borge–Holthoefer and Arenas (2010) and Baronchelli et al. (2013). Once the network representation is set up, percolation analysis can be used to study semantic memory and how it is affected by failures, i.e., cutting links in the network, and how it facilitates dynamical processes that act upon it. The former was used to study cognitive phenomena in patients with Alzheimer’s disease (Borge-Holthoefer et al. 2011), while the latter was investigated in the context of memory retrieval (Arenas et al. 2012).

Kenett et al. (2018) proposed a percolation approach to study high-level cognition, and to test the hypothesis that the semantic network of highly creative people is more flexible than that of low creative individuals and, therefore, more robust. Semantic networks of low semantic creative (LSC) and high semantic creative (HSC) persons have been reported by Kenett et al. (2014), which were utilized by Kenett et al. (2018). The two networks, each with about 4,000 links, consisted of 96 cue words that were divided into groups of 4 concrete words from 24 categories, including fruits, musical instruments, vehicles, and others, and representing a priori components of the two networks. Kenett et al. (2018) modified the two networks by taking into account the number of participants who generated the associative features that determine the strength of semantic similarity between the cue words. The more similar associative response generated and the larger number of participants who generated the association responses to a pair of cue words, the stronger the link between this pair of cue words is.

Kenett et al. (2018) normalized the links in each network according to the mean number of associations per cue word, so as to remove the effect of the HSC individuals who generate a higher number of associations per cue words, when compared with the LSC individuals. The links between all pairs of cue words define a symmetric correlation matrix whose (i, j) entry represents the semantic similarity between cue words i and j, which can be studied in terms of an adjacency matrix of a weighted, undirected network in which each cue word is represented by a node, and a link between two cue words represents the semantic weight between them.

As usual, percolation analysis of the two networks was carried out by removing their links according to a pre-set threshold $\mathcal{T}$, such that all links with a weight smaller than the threshold were removed, the size (number of nodes) of the giant component in the network was determined, and the detached components (clusters) were identified. $\mathcal{T}$ was varied from the smallest weight in the network, the initial threshold $\mathcal{T}_i$, to a weight strength in which the giant component was smaller than three nodes, the final threshold $\mathcal{T}_f$, with its resolution selected according to the smallest difference between the sorted weights. Thus, one has a variety of step sizes that are the number of nodes of a component that disconnects from the GC in each step, most of which were groups of several nodes.

An illuminating quantity that was calculated is the percolation integral $I_p$, defined by

$$\begin{aligned} I_p=\int _{\mathcal{T}_i}^{\mathcal{T}_f} C_G(x)dx=\sum _{\mathcal{T}_i}^{\mathcal{T}_f} C_G(\mathcal{T})\mathcal{T}_r\;, \end{aligned}$$

(18.24)

where $C_G$ is the size of the GC at a threshold and $\mathcal{T}_r$ is the threshold resolution. $I_p$ measures the speed by which the GC breaks. If the network breaks at a low value of the thresholds and has a steep percolation curve (see below), it will have a smaller $I_p$ than a network with the same size of the GC that breaks with large thresholds and a flat percolation curve.

The effect of noise is important and must be studied. If one adds noise to the weights of the links and, for example, makes weaker weights stronger, and as a result obtains a value of $I_p$ that is dramatically different from its pre-noise value, it implies that the percolation analysis has no physical significance. To test this, Kenett et al. (2018) computed $I_p$ over 500 realizations of the networks before carrying out the percolation analysis. The realizations were generated by adding a Gaussian noise to the network with a mean value of zero and a standard deviation $\sigma $ that varied between $10^{-4}$ and $10^{-2}$.

In addition to studying the effect of noise, analysis of the networks whose links were shuffled was also carried out in order to understand the effect of the structure of the network on the percolation analysis. In the shuffling, two links from two pairs of nodes are selected at random and are exchanged. For example, nodes $n_1$ and $n_2$ with link strength 0.5 and nodes $n_3$ and $n_4$ with link strength 0.7 are exchanged, so that the new network topology would be one in which $n_1$ and $n_3$ are connected with strength 0.5 and $n_2$ and $n_4$ with strength 0.7. This process was reiterated ten times the number of links in the network (i.e., about $\sim 40,000$ times) to ensure that most of the links in the network were replaced, and was also reiterated to generate 500 realizations, for each of which $I_p$ was computed.

The hypothesis is that the semantic network of HSC individuals is more flexible than the LSC network and, therefore, their network will be more robust to link removal, hence breaking apart slower than the LSC network. The analysis of Kenett et al. (2018) revealed that the GC of the LSC network breaks apart faster at lower thresholds than the HSC network. This is shown in Fig. 18.7A. For the same threshold values, the LSC giant cluster is mostly smaller than its counterpart in the HSC network, indicating that the latter is more robust and its components are better connected. The percolation integral $I_p$ turned out to be 14.76 for the LSC network, and 15.64 for the HSC one, computed for the type of curves shown in Fig. 18.7A.

The effect of noise is shown in Fig. 18.7B, indicating that both the LSC and HSC networks have a stable percolation process with significantly different integrals. The effect of link shuffling is demonstrated in Fig. 18.7C. Independent analysis of samples with the t test, carried out for the distribution of $I_p$ values of the LSC networks and its comparison with the HSC ones, indicated that the percolation shuffled integral of the former is even slightly larger than that of the latter, implying that the robustness of the HSC network is more influenced by the shuffling of its structure. This suggests strongly that the difference between the empirical $I_p$ of the two groups is driven by the structure of the networks, rather than by the link weights.

18.12 Cardiac Fibrosis and Arrhythmia

Cardiac contraction is controlled by abnormal propagation of nonlinear excitation waves in the heart, which can create life-threatening cardiac arrhythmias. Typically, such arrhythmias develop when the waves form turbulent rotational activity, or reentry. A well-known example of a pathology that leads to reentry formation in cardiac tissues is fibrosis (see, for example, Nguyen et al. 2014), which is the presence of a large number of nonexcitable cells that are scattered in the cardiac tissue, formed by myocytes that are excitable cells, inhomogeneously distributed, and are connected by gap junctions. The extracellular matrix, non-conducting gap junctions, blood vessels or fibroblasts that nonexcitable cells embedded in the tissue, can disturb the gap-junction coupling among the myocytes. Non-conducting connections in cardiac tissue strongly affect wave propagation and increases the possibility of arrhythmias. An important issue is the link between the pattern of fibrosis and the onset of arrhythmia, as it has practical implications for cardiology and development of arrhythmia treatment.

Alonso and Bär (2013) first suggested that heterogeneities in cardiac tissue may provide an anatomical mechanism for reentry. They considered a 2D cellular network with some link being non-conducting whose fraction was $\phi $. Thus, $\phi =0$ and 1 corresponded, respectively, to completely homogeneous tissue, and a set of disconnected myocytes. For small $\phi $ a wave propagates through the network, whereas for large $\phi $ the wave cannot propagate. Thus, as the reader knows, there is a percolation threshold$\phi _c$ where both behaviors interchange. Alonso and Bär (2013) carried out numerical simulations with the model in which they varied $\phi $. The governing equation for propagation of action potential at node i (or a cell with its center at i) is given by

$$\begin{aligned} \frac{\partial V_i}{\partial t}=-(I_{fi}+I_{so}+I_{si})+\sum _{j=1}^Zg_{ij} (V_j-V_i)\;, \end{aligned}$$

(18.25)

where $I_{fi}$, $I_{so}$ and $I_{si}$ are, respectively, fast inward sodium, slow outward potassium, and slow inward calcium ionic currents, which are all related to the potential $V_i$ (Fenton and Karma 1998), and $g_{ij}$ represents the conductance between i and j, taking on a constant value if i and j are connected, and zero if i and j are not linked.

The simulations indicated that a wave that propagates in homogeneous tissues and interacts with a heterogenous region—a damaged region on the cardiac muscle—can generate reexcitations through a small-scale reentry process, which may emerge as ectopic beats—extra heartbeats that are caused by a signal to the upper chambers of the heart (the atria) from an abnormal electrical focus—or as reentries. Generating such reexcitations requires the combination of two factors inside the damaged area, i.e., a fraction of non-conducting links close to but above the percolation threshold, and a change of the local properties toward weaker excitability. In addition, the probability of reentry was shown to be strongly correlated with the probability of formation of clusters with a size larger than a typical size that emerge near $\phi _c$.

In another paper, Alonso et al. (2016) carried out the same type of simulation in a 3D model, represented by a thin slab of cardiac tissue, in order to take into account the fact that atrial tissue is much thinner than ventricular one. Their simulations indicated the strong and non-trivial effect of the thickness, even for thin tissues, on the probability of micro-reentries and ectopic beat generation, and a strong correlation of the occurrence of micro-reentry with $\phi _c$. There was also a qualitative agreement between the 3D simulated electrograms in the fibrotic region with the experimentally observed complex fractional atrial electrograms (CFASs).

Vigmond et al. (2016) used an approach similar to Alonso and co-workers, except that they used 2D synthetic fibrosis patterns, discretized like a square network. The pattern was made heterogeneous by a percolation algorithm, namely, by generating a random number R for each elementary block and removing it if R was less than a threshold set for the element, representing the fibrotic density. Isotropic fibrosis was generated by using a uniform threshold for the removal of the elements, whereas anisotropic fibrosis was produced by setting the threshold of removal for even rows much higher than for the odd ones and, therefore, even rows represented disruptions to lateral coupling. The simulations indicated that reentry may be obtained with a two-phase behavior (conducting and non-conducting) and depends on fibrotic density. The CFAEs were recorded above fibrotic regions and, consistent with clinical data, electrogram duration and fractionation increased with more rapid pacing.

Since the early twentieth Century, two conditions have been widely recognized for the formation of reentry (Mines,^{Footnote 4} 1913), namely, the existence of a critical wavelength, and the so-called unidirectional blocks. In the case of the former, in order for a rotation to continue, its period must be longer than what is referred to as the refractory time $t_r$, which is the period of time in which an organ cannot repeat a particular action. In other words, $t_r$ is the amount of time that an excitable membrane requires to become ready for a second stimulus, once it returns to its resting state following an excitation. Thus, the existence of such a temporal condition suggests that the wavelength of the rotational activity must be larger than a critical value, namely, $vt_r$, where $v_r$ is the wave velocity. Structural heterogeneity, such as clusters of fibrosis, provides a natural morphology around which the rotation can occur.

The second condition, the existence of unidirectional blocks, is often necessary for the initiation of reentry (Kléber and Rudy 2004). The wave is locally blocked from propagating in one direction, but can do so in the opposite direction. Morphological heterogeneity created by fibrosis contributes to the existence of unidirectional block. As pointed out above, Alonso and Bär (2012) also showed that the probability of reentry in fibrotic tissue is strongly correlated with the probability of the formation of fibrotic clusters with a size larger than the typical size.

Pashakhanloo and Panfilov (2021) introduced a model for predicting the probability of reentry formation based on the size distribution of what they referred to as minimal functional clusters in the tissue, which is essentially a percolation model and represents a refinement of the earlier model by Alonso and co-workers. In their model, one begins with a $L\times L$ square lattice in which fibrosis is modeled by randomly removing elements from the lattice with probability p. For any given p, the simulations are carried out using the Beeler–Reuter–Drouhard–Roberge ionic model, which is similar to Eq. (18.25). The stimulus is applied from the lattice’s boundary in order to generate planar waves. Reentry was considered successful if the activation in the tissue lasted longer than a certain time, following application of the stimulus. The model successfully simulated the formation of reentry, when macroscopic rotational activity was initiated by a microscopic unidirectional block, such that the wave was blocked from propagating from the left to right direction, but could do so in the opposite direction.

A cluster is defined (Alonso and Bär 2013) as a set of connected inexcitable elements that form an obstacle, a purely topological characteristic of the texture that, without any conduction blocks, would solely determine the path of the propagating wave. Such blocks, in the presence of fibrosis-induced conduction blocks, form at the junction of, or close proximity to, two or more regions of fibrosis that connect the fibrotic clusters to form a larger cluster of unexcitable tissue, which is what Pashakhanloo and Panfilov referred to as the functional cluster, which is a more refined concept of fibrotic cluster, as it takes into account not only the morphology of the texture but also the properties of propagating waves, i.e., conduction blocks.

Pashakhanloo and Panfilov studied the size distribution of the functional clusters. A “correlation length” $\xi $ was defined in a manner similar to that for percolation clusters (see Chap. 2),

$$\begin{aligned} \xi ^2=\frac{2\sum _s R_s^2s^2n_s}{\sum _s s^2n_s}\;, \end{aligned}$$

(18.26)

where $n_s$ is the number of clusters of size s, and $R_s$ is the radius of gyration of the cluster. $\xi $ depends not only p but also on excitability e. The important point is that $\xi $ for the functional clusters is not equal to that of the fibrotic clusters, and is strongly affected by excitability. Moreover, the percolation threshold at which a spanning functional or fibrotic cluster forms was found to agree with the wave percolation threshold, i.e., the probability of formation of a total wave block in the lattice. The functional clusters manifest the dependence of percolation on the excitability, such that lower excitability resulted in a lower threshold, because a larger number of blocks existed.

Dynamically, after the occurrence of a unidirectional block, the tissue behind it can be re-excited from the opposite direction, if the propagation delay, caused by wave propagation around clusters to reach the conduction block from the opposite side, is longer than the refractory period $t_r$. But because there are multiple paths for the wavefront to reach there, the delay should be determined by the shortest path. This gives rise (Pashakhanloo and Panfilov 2021) to a minimal functional cluster associated with a unidirectional block, which is constructed by disconnecting the encompassing functional cluster at the location of the block, and identifying among the resulting subclusters the one with the smallest perimeter. The physical significance of the cluster with minimum perimeter is that the wave must—partially or fully—travel along its perimeter in order to reach the other side of the block, and possibly cause reexcitation.

Finally, we note that Rabinovitch et al. (2021) utilized a 2D cellular automata model to simulation transport, percolation properties, and tortuosity in heart, and Falkenberg et al. (2019) used a similar, but more complex, 3D model to study mechanisms of atrial fibrillation (AF). Their simulations produced spontaneous emergence of the AF that was in agreement with the clinically observed diversity of activation. Moreover, the simulations indicated that the difference in surface activation patterns is a direct consequence of the thickness of the discrete network of heart muscle cells in which electrical signals flow and percolate to reach the imaged surface.

18.13 Connectivity of Temporal Hierarchical Mobility Networks During COVID-19

The coronavirus disease 2019 (COVID-19), an unprecedented pandemic, affected more than 200 countries. As of the time of the final editing of this section (12 August 2022), the total number of infected people has been 585,086,861, with 6,422,914 having lost their lives (COVID-19, as reported by the World Health Organization). In order to understand how the infectious disease spreads, one must study its relation with the mobility pattern of a population. For example, some studies indicated (see, for example, Kraemer et al. 2020) that the early spatial patterns of COVID-19 infection in China correlated with the population mobility fluxes, and that the correlation decreased as some local control policies for limiting mobility were put in place.

Deng et al. (2021) demonstrated the relevance of percolation process and phase transition to human mobility networks on the county level. Their study indicated the possibility of developing effective strategies for controlling the mobility flow at critical bridges [connecting urban areas] and containing the transmission of COVID-19. To gain a deeper understanding of the percolation effect, one must construct mobility-based networks in order to accurately model the transmission of COVID-19. To do so, one must consider some key properties: (a) The networks must have proper resolution in order to account for both the spatial and temporal dynamics of mobility networks. For example, some studies aggregated human mobility on the census block group or county levels, whereas aggregation removes details and information on other levels, resulting in inaccurate prediction for the transmission patterns. (b) Some studies pointed to the effectiveness of inter-city travel restrictions for reducing the imported COVID-19 incidence rate, but one must gain a better understanding of how local measures and responses could change.

He et al. (2022) studied the percolation of the human mobility network that included over 175,000 of the census block groups in the United States. In their network model each bond between two sites, which represent the blocks, is given a weight. Two block groups have a connection (a bond between them) on any day in which people travel between them. If more people travel between them, the weight increases. To analyze the properties of the network, He et al. (2022) proposed a hierarchical structure of the mobility network to understand human dynamics during the pandemic, in which the network of inter-Metropolitan Statistical Area (MSA) was composed of 378 MSAs, with each comprising of many census block groups, which was at the intra-MSA level. To implement the percolation process in the network, He et al. (2022) deleted a fraction q of the weakest links (lowest weights) and determined the giant component of the network, since it is easier to break the weak links than the strong ones. They demonstrated the existence of a critical value $q_c$ of q, such that when $q>q_c$, the network breaks down abruptly because the large enough second largest component of the network becomes disconnected from the giant cluster. The value of $q_c$ is dynamic, and He et al. (2022) computed its daily value for each intra-MSA and inter-MSA.

The study of He et al. (2022) pointed to three crucial features of percolation in dynamic mobility networks. One is the hierarchical structure of the temporal mobility networks based on the concept of cluster formation in the classic percolation theory. The second feature is the universality of phase transition in mobility networks at $q_c$, a concept that by now is very familiar to the reader of this book, demonstrating that it is the percolation process that controls the universality of transitions in such networks on various levels, from block groups to inter-MSAs. The third feature is the association between critical threshold $q_c$ and the key characteristics of local MSAs, which improves the ability for assessing and predicting the vulnerability of mobility networks.

A related problem is vaccinating people and developing a strategy for maximum efficiency so that the society can eventually reach herd immunity. One such strategy is vaccinating those that have more exposure to others, which has been shown theoretically to be highly effective. The problem with this strategy is that identifying such individuals—sometimes referred to as “hot spotting”—is difficult, hence preventing its large-scale implementation. Penney et al. (2021) proposed an approach for using the technology underlying digital contact tracing in order to boost vaccine coverage. To do so, Penney et al. (2021) modeled the spread of disease using percolation in random graphs, described in Chap. 17. They defined a quantitative measure, the efficiency, which they defined as the percentage decrease in the reproduction number per percentage of the population vaccinated. Penney et al. (2021) demonstrated that optimal implementations of their percolation model can achieve herd immunity with as little as half as many vaccine doses as a non-targeted strategy.

18.14 Protein Structure

Theoretically, there are $10^{400}$ possible sequences for a protein of medium size. But, the sequences of only about $10^8$ proteins are currently known, whereas the number of extant sequences is estimated to be $10^{34}$, and it is believed that over the past four billion years of evolution up to $10^{43}$ different protein sequences may have been explored. Such statistics imply that the Darwinian protein evolution, based on mutation of the genotype and subsequent natural selection of the phenotype, excludes the possibility that extant sequences are randomly scattered in the theoretical sequence space. Instead, they are expected to form a connected network in which functional sequences and mutations form the nodes and links or edges (Smith 1970). Then, the fundamental questions (Smith 1970) are whether all existing proteins are part of a single network with a single starting point, and what fraction of the functional sequence space has been explored. One is also interested to learn how large is the space of functional, but never-born, proteins.

The $10^8$ known protein sequences are not equally distributed, but cluster into homologous families (Rappoport et al. 2012). Due to their sparsity (separation between the clusters), however, most clusters have neighboring nodes that differ by multiple mutations. The sparsity may be explained by the extinction of ancestor sequences (Thornton 2004). Alternatively, one may predict that as our knowledge of known sequences of proteins expands, the sparsity decreases.

Buchholz et al. (2017) used percolation theory to analyze the sequence space of known proteins. In particular, they utilized the power-law cluster-size distribution of percolation, $n_s\sim s^{-\tau }$ (see Chap. 2), to study the clusters in the network of the extant sequences. They began by updating the databases on $\alpha /\beta $ hydrolases (abH, 395,000 sequences), cytochrome P450 monooxygenases (CYP, 53,000 sequences), thiamine diphosphate-dependent decarboxylases (DC, 39,000 sequences), and $\beta $-hydroxyacid dehydrogenases/imine reductases (bHAD, 31,000 sequences) (see Buchholz et al. 2017 for references to the databases)). Representative sequences for each homologous family were selected as seed sequences, after which family databases for short-chain dehydrogenases/reductases (SDR, 141,000 sequences) and $\omega $-transaminases (oTA, 121,000 sequences) were established based on seed sequences derived from literature. Then, sequence identities of high-scoring sequence pairs were calculated for each protein database, and sequence pairs with a distinct sequence identity cutoff were clustered.

The cluster-size distributions $n_s$ for the six datasets were then constructed. To do so, $n_s$ was analyzed for cluster sizes, $s=1,2,\ldots , 1000$ and, because for large s the data become increasingly sparse, the corresponding histogram distributions were constructed by counting the number of clusters $n_{i,j}=\sum _{s=i}^j n_s$ for cluster sizes $i\le s\le j$. Three theoretical models were then tested against the data: (a) a Gaussian distribution; (b) the exponential distribution, $n_s\sim \exp (-bs)$, where b is a fitting parameter, and (c) the power-law distribution, as predicted by percolation theory.

When the data were plotted in log-log fashion, they differed considerably. $\log n_s$ according to the Gaussian distribution increased gradually with $\log s$ and decayed rapidly for s larger than the mean value of the distribution. The exponential distribution decays rapidly, of course, for all s, whereas for the power-law distribution $\log n_s$ depends linearly on $\log s$. Figure 18.8 presents the results for the six sets of the data that follow the power-law distribution. The exponent $\tau $ was found to be between 2.4 and 3.3. As described in Chap. 17, biological networks represent structures with high dimensionality. Therefore, the lower value of 2.4 is consistent with the mean-field value, $\tau =5/2$ (see Table 2.3). As described in Chap. 17, for scale-free networks described by the degree distribution (17.7), the exponent $\tau $ is given by, $\tau =(2\gamma -3)/(\gamma -2)$, where $\gamma $ represents the exponent of the power law (17.7). Thus, if we use $\tau \approx 3.3$, we obtain $\gamma \approx 2.77$. In other words, the protein sequences may be represented by a scale-free network with the exponent $\gamma \approx 2.77$.

For related work see Deb (2009), Weber and Pande (2015), and Peng et al. (2015).

18.15 Molecular Motors and Mechanical Properties of Cells and Active Gels

Molecular motors in the cytoskeletal filamentous actin (F-actin) network generate internal stresses that determine the mechanical properties of cells. Some aspects of this were already described in 18.8. It is due to motor activity that the cells have the ability to contract the surrounding extracellular matrix (ECM), which is essentially a biopolymer network whose elastic properties are strongly influenced by the active contractility. This happens in both reconstituted intracellular F-actin networks with myosin motors (see, for example, Mizuno et al. 2007; Gordon et al. 2012) and in the ECM with contractile cells (Lam et al. 2011). Active polymer networks are heterogeneous and contain soft or floppy modes of deformation (see Chap. 11). Motor-induced contractile stresses are coupled to the soft modes (Broedersz and MacKintosh 2011), which produce a nonlinear elastic response, which is different from the nonlinearities arising from single fiber elasticity. In addition, the coupling casts doubts on the presumed equivalence of internal—motor—and external stress. As a result, as has been emphasized throughout this book, their properties cannot be quantitatively described by the classical linear continuum theory of elasticity.

As described in Chap. 11, one can form elastic networks by interconnected straight fibers (bonds) with linear stretching and bending elasticity. If a fraction of the fibers is removed, one has an elastic percolation model. Sheinman et al. (2012; see also Sheinman et al. 2015) proposed such a model to study the effect of motor-generated stresses in percolating fiber networks. They used an FCC network that allowed them to study the effect of the fiber connectivity over an extended range, since the coordination number of the intact FCC network is $Z=12$. To introduce motor activity, they incorporated contractile, static, and strain-independent force dipoles between neighboring nodes of the network, inserted randomly with a probability q, while the fibers were modeled as linear elastic bonds with a stretching modulus $\alpha _1$ and bending rigidity $\alpha _2$ (see Chap. 11, where the two were denoted, respectively, by $\alpha $ and $\gamma $). Assuming that the equilibrium length of the fibers is unity, the total energy of the network is given by (see also Chap. 11)

$$\begin{aligned} \mathcal{H}=\frac{1}{2}\alpha _1\sum _{\langle ij\rangle }e_{ij}(|\textbf{r}_{ij}|-1)^2+ \frac{1}{2}\alpha _2\sum _{\langle ijk\rangle }e_{ij}e_{jk}\left( \frac{\textbf{r}_{ij}\times \textbf{r}_{jk}}{|\textbf{r}_{ij}||\textbf{r}_{jk}|}\right) +f\sum _{\langle ij\rangle }q_{ij}|\textbf{r}_{ij}|\;. \end{aligned}$$

(18.27)

Here, $\textbf{r}_{ij}=\textbf{u}_i-\textbf{u}_j$, with $\textbf{u}_i$ denoting displacement) of node i, $e_{ij}=1$ or 0 if the fiber is, respectively, present or cut, f is the dipole rigidity, and $q_{ij}=1$, if a motor acts between nodes i and j, or $q_{ij}=0$, otherwise.

The model was then studied by an effective-medium approximation (EMA)—poor man’s percolation—and numerical simulation. To develop the EMA, Sherman (2012) set $\alpha _2=0$, which reduced their network to one somewhat similar to a central-force spring network described in Chap. 11. The development of their EMA paralleled what was described in Chap. 11, but somewhat more complex, due to the additional term in the total energy $\mathcal{H}$, Eq. (18.27), representing the contribution of the motors. Their final EMA equation for the effective stretching modulus $\alpha _e$ is given by

$$\begin{aligned} \int _0^\infty \frac{\alpha _{ij}-\alpha _{1e}(\sigma _M)}{1/u_\textrm{EM}+\alpha _{ij}- \alpha _{1e}(\sigma _M)}h(\alpha _{ij})d\alpha _{ij}=0\;, \end{aligned}$$

(18.28)

where $u_\textrm{EM}$ is the displacement of a bond in the unperturbed effective medium due to a unit force acting on a bond, $\alpha _{ij}$ is the stretching modulus of fiber ij, $h(\alpha _{ij})$ represents the probability density function of $\alpha _{ij}$, and $\sigma _M=\sqrt{8}\;qf$ within the EMA. Thus, with a percolation distribution, $h(\alpha _{ij})=p\delta (\alpha _{ij}-1)+(1-p) \delta (\alpha _{ij})$, where $\alpha _{ij}=\alpha _1e_{ij}$, one obtains the following expression for the shear modulus of the network:

$$\begin{aligned} \mu _\textrm{EM}=\frac{5\sqrt{2}}{72}\alpha _{1e}+\frac{5}{6}\sigma _M\;. \end{aligned}$$

(18.29)

Recall from Chap. 11 that a central-force FCC network at the percolation threshold $p_c$ has a mean connectivity of $Z_\textrm{CF}=6$ [see Eq. (11.40)], and that for coordination numbers $Z<Z_\textrm{CF}$ the network is not rigid. In the present case, however, the presence of the motors means that even below $Z_\textrm{CF}$ the network is rigid with a finite effective shear modulus. Thus, one may view motor stress as an external field that stabilizes floppy networks. Numerical simulations of Sherman (2012) indicated that close to $Z_\textrm{CF}$ one has the anomalous scaling,

$$\begin{aligned} \mu _e\sim \alpha _e^{1-x}\sigma _M^x\;, \end{aligned}$$

(18.30)

with $x\approx 0.4$. The EMA predicts that, $x=1/2$. If the bending rigidity is not ignored, i.e., if $\alpha _2\ne 0$ in Eq. (27), then, Sherman (2012) showed that

$$\begin{aligned} \mu _e\sim \alpha _{2e}^{1-y}\sigma _M^y\;, \end{aligned}$$

(18.31)

with their numerical simulations yielding $y\approx 0.6$, where $\alpha _{2e}$ is the effective bending rigidity.

In another paper from the same group (Alvarado et al. 2013), extensive experimental data were analyzed to show that myosin motors contract crosslinks actin polymer networks with a scale-free size distribution (see Chap. 17), and that the effect occurs over a broad range of crosslink concentrations. Their analysis indicated that the motors reduce connectivity of the networks by unbinding the crosslinks, and that, in order to coordinate global contractions, motor activity should be low. Otherwise, they force initially well-connected networks to a critical state at which rupture forms across the entire network.

To understand the interplay between motor activity and the connectivity of active cytoskeletal network, Alvarado et al. (2013) varied the density of myosin motors and fascin crosslinks, given by the molar ratios $R_M$ and $R_C$, and prepared a series of active networks with constant myosin activity, $R_M= 0.01$, and gradually increased $R_C$. The experiments indicated that even at low $R_C$, the motors contract actin networks, albeit on small length scales. As $R_C$ was increased, contraction occurred on larger length scales. Eventually, the motors break up the network into multiple disjoint clusters. At still higher $R_C$, motor activity contracted the entire network into a single dense cluster. The transition from local to macroscopic contraction represented a percolation transition: Below the transition point, the network is only locally correlated without long-range connectivity, whereas above a well-defined critical connectivity the system establish global correlations and corrections.

To demonstrate the connection to the percolation transition, Alvarado et al. (2013) noted that, according to percolation theory, the largest cluster of size $\xi _1$ (which is the percolation correlation length $\xi _p$; see Chap. 2) increases monotonically with p, the probability of making connections between the nodes, whereas the second-largest cluster with a size $\xi _2$ should exhibit a peak right at the connectivity threshold, with $\xi _1$ and $\xi _2$ both approaching the system’s linear size L. The data of Alvarado et al. (2013) was in agreement with this picture; see Fig. 18.9.

Another test was carried out by Alvarado et al. (2013) to further establish the connection of the problem to percolation theory. It was based on the cluster-size distribution of percolation networks. As described in Chap. 2, close to the connectivity threshold, one has a power-law size distribution, with an exponent (in 2D) of, $\tau =187/91\approx 2.05$ [see Eq. (2.14) and Table 2.3]. To see whether this prediction is borne out by their data, Alvarado et al. (2013) looked for networks that satisfied the constraints, $\xi _1 \sim \xi _2\sim L$, and replotted all the data separately in the $\xi _1-\xi _2 $-space. This is shown in Fig. 18.10a. Since by definition, $\xi _2<\xi _1$, all samples are within a triangle in $\xi _1-\xi _2$-space. Figure 18.10a identifies clearly the samples at the triangle’s peak, where $\xi _1\sim \xi _2 \sim L$, and represents the critically connected network. Samples to the left of the peak are for low $R_C$, which is identified as the local contraction regime, whereas those to the right represent samples with high $R_C$ and, therefore, the global contraction regime.

Next, the entire distribution of cluster sizes was plotted; see Fig. 18.10b. The data are again consistent with percolation theory, since the critically connected regime exhibits a power-law cluster-size distribution over more than two orders of magnitude in measured area. The power-law exponent is $\tau \approx 1.9$. Although the estimate is only about 5% smaller than the predicted value of 2.05 in 2D, $\tau $ for the percolation clusters cannot be less than 2.

If one reverses the roles of occupied and vacant sites or bonds in a percolation lattice, the problem can be thought of as finding the distribution of hole sizes within a single large cluster (Hu et al. 2016). Thus, the holes studied are not simple percolation holes, but are within the surrounding backbone. Using scaling analysis, Hu et al. (2016) suggested that, $\tau =1+D_\textrm{BB}/2$ for the backbone holes, where $D_\textrm{BB}\simeq 1.64$ is the fractal dimension of percolation backbone in 2D (see Chap. 2). Thus, one obtains, $\tau \simeq 1.82$. If, however, one considers a model of simple holes within a percolation cluster, one obtains (Hu et al. 2016), $\tau =1+D_f/2$, where $D_f=91/48\simeq 1.895$ is the fractal dimension of 2D sample-spanning percolation cluster. Therefore, one obtains $\tau \simeq 1.94$, consistent with the experimental estimate (see Fig. 18.10).

Note also that the distributions of the other two regimes also agree with percolation theory, because the local contraction regime has a short-tail distribution with a sharp cutoff. On the other hand, the global contraction regime exhibits a bimodal distribution with two well-separated length scales, one with the percolating cluster with a size $\xi _1\sim L$, and a second one for small disjointed clusters with a typical size, $\xi _2\ll L$L.

In addition, Alvarado et al. (2013) demonstrated that restructuring the network breaks it into clusters because motor stresses unbind crosslinks. Thus, one may predict that increased motor activity—increased $R_M$ or motor force—should generate smaller clusters. The simulations and experiments of Alvarado et al. (2013) confirmed the predictions. If, however, motor forces become larger than the unbinding threshold of the crosslinks that maintain network stability, crosslink constraints may fail across the system, similar to fracture of heterogeneous networks (Sahimi and Goddard 1986; de Arcangelis et al. 1989; see Sahimi 2003b for a comprehensive review) and motors are free to slide actin filaments. Actin–myosin gels that are well-connected respond to motor sliding by undergoing a global contraction in which the meshwork collapses to a tightly packed cluster with a density higher than that of the the surrounding fluid.

For discussions of active stress in fiber networks in the context of biological functions see also Dasanayake et al. (2011), Lenz (2014), and Ronceray et al. (2016). A comprehensive review of these fascinating problems and their relationship with percolation theory is given by Alvarado et al. (2017). The percolation model of Sheinman et al. (2012, 2015) has been studied by several groups; see, for example, Böttcher et al. (2016), and Lee and Pruessner (2016) described earlier.

18.16 Percolation and Evolution

One may describe evolution as a percolation problem in a multidimensional fitness landscape (Gavrilets and Gravner 1997; Gravner et al. 2007; Mustonen and Lassig 2009; Catalan et al. 2017; see also Katsnelson et al. 2018). In this model selectable traits or genes correspond to “dimensions,” and their number is in the thousands, the model corresponds to a high-dimensional percolation. Thus, one may imagine a landscape with mountain ridges and plateaus that are partially covered with water. If the paths that go under water never continue, in the sense of being sample-spanning, whereas there are accessible paths (in the percolation sense) above the water level whose existence depends on the evolutionary temperature, or effective population size, then, one has a percolation model.

The percolation model becomes clearer if we consider the critical (sample-spanning) percolation cluster consisting of paths that follow the “shores”—the lines of minimal fitness necessary for survival. There is considerable evidence that, due to the cost of selection, actual evolution does not deviate much from such paths (Barton 1995; Bell 2013). Since the percolation model of evolution is one in a high-dimensional space, which is well represented by a Bethe lattice or Cayley tree (see Chap. 2), the implication is that the SSC—or the GC since the model is in high-dimensional space (see Chap. 17)—has a tree structure without double paths or dead ends and, as such, is an “optimal” simple line. Moreover, recall from Chap. 2 that the tree–like structure of the critical percolation cluster (in a high-dimensional space) has only one route between any two points on such a landscape.

Using such concepts, Gravner et al. (2007) studied two such percolation models of evolution. One represented a phenotype space in which individuals were characterized by a large number n of continuously varying traits. If the fitnesses are randomly assigned, viable phenotypes form a giant connected cluster that percolates throughout the phenotype space, if the viability probability is larger than $1/2^n$. In the second model, genotype-to-phenotype and phenotype-to-fitness maps were explicitly accounted for, allowing for neutrality at both phenotype and fitness levels. This generated a fitness landscape with tunable correlation length. Phenotypic neutrality and correlation between fitnesses may reduce the percolation threshold (see Chap. 3), and correlations at the phase transition point between local and global are most conducive to the formation of the giant cluster.

18.17 Microbial Communications

It was assumed for a long time that bacteria are “loners” that act on their own as organisms without any interactions with others. Accumulated experimental evidence over the years forced researchers to rethink this view. Not only do bacteria communicate with each other through signal transmission but also do so over length scales that are much larger than their own individual cell sizes. Such communications result in the formation of biolfilms and producing antibiotics, as well as pathogenesis.

Communication between bacteria is based on chemical signaling molecules that are called autoinducers. Such molecules regulate bacterial gene expression in a process that is referred to as quorum sensing. Similar to various languages, the signals also vary between species, with some bacteria being able to interpret many different signals, whereas others could respond to only a select few. It is through quorum sensing that individual bacteria within colonies coordinate and carry out colony-wide functions, such as sporulation, bioluminescence, and biofilm formation.

The way quorum sensing works is as follows (American Society for Microbiology, 2020). Individual bacterium synthesizes autoinducers. Acyl-homoserine lactone autoinducers, produced by Gram-negative bacteria, diffuse passively (no reaction) through their thin cell wall, whereas Gram-positive bacterial autoinducers that consist of peptide must be actively transported through their peptidoglycan cell wall using a transporting system called the adenosine triphosphate (ATP)-binding cassette. As the bacteria keep reproducing, more individual cells also produce autoinducers that move out of the cells, hence increasing their concentration in the extracellular matrix, eventually reaching a critical threshold for their mass that makes it energetically unfavorable for intracellular autoinducers to keep leaving the cell, thus increasing their intracellular concentration. As their concentration in intracellular increases, autoinducers bind to their receptors and trigger signaling cascades that alter transcription factor activity and, therefore, gene expression, which for many bacteria, includes downregulating autoinducer synthesis in a negative feedback loop.

For example, the Gram-negative pathogenic bacteria Vibrio cholerae (V. cholerae) uses quorum sensing for virulence during a cholera infection by building biofilms to help transport nutrients between colonies, while protecting them at the same time. It is due to their ability of forming biofilms within the host that bacteria’s reproduction cycle and eventual secretion of cholera toxin are successful, one of two virulence factors that contribute to between 21,000 and 143,000 cholera deaths worldwide each year.

To transmit signals through, for example, electrical signals in the nervous system, and coordinate cellular actions at long distances, bacterial communities must overcome two major challenges. One is how to achieve reliable signal transmission to desired target sites. The second challenge is due to the fact that bacterial communities manifest significant cell-to-cell heterogeneity, which may act as obstacles for long-range signal propagation (Raj and van Oudenaarden 2008; Symmons and Raj 2016). Thus, if, for example, only a fraction of the cells contribute to signal transmission, the cell-to-cell heterogeneity could result in the propagating signal to dieing out before reaching its target (Cao et al. 1999; Alonso et al. 2016). Therefore, the question of the effect of heterogeneity on bacterial communities in long-range signal transmission is important. Then, once one talks about heterogeneity and only a fraction of the cells contributing to signal propagation, the relevance of percolation theory becomes immediately clear.

Larkin et al. (2018) developed a percolation model to describe the connectivity and clustering statistics of firing cells in the biofilm. They used a triangular lattice, since cells in their experimental biofilms had a modal value of six neighbors. Thus, the cells had six nearest neighbors to which they connect, forming triangles. If a fraction $\phi $ of the cells are firing, then, at the percolation threshold $\phi _c$ a SSC of such cells forms. This corresponds to site percolation. Larking et al. utilized a $L\times W$ lattice with $L=30$ and $W=200$ cells, roughly corresponding to the size of their experimental systems, with the signal propagating in the L direction. Due to the asymmetric dimensions of the lattice, they estimated that $\phi _c \approx 0.45$, about 10% less than the exact site percolation threshold of the triangular lattice, $\phi _c=1/2$ (see Chap. 2). To describe the dynamics of signal propagation between firing cells, Larkin et al. used the FitzHugh–Nagumo model,^{Footnote 5} which is a slightly more complex model than what is described by Eq. (18.25).

Larking et al. also carried out experiments to compare with their simulations. The bacteria used was Bacillus (B.) subtilis NCIB 3610. B. subtilis, commonly found in soil, is a gram-positive model organism for basic research and an industrial workhorse, while B. subtilis NCIB 3610 is an undomesticated strain that exhibits phenotypes lost from the more common domesticated laboratory strains. In order to measure membrane potential of individual bacteria within the biofilms during electrochemical signaling, Larkin et al. (2018) used the fluorescent reporter thioflavin-T (ThT), which acts as a Nernstian membrane potential indicator, such that the higher the membrane potential of the cell, the larger is the amplitude of the fluorescent ThT signal. Measurements at the scale of a single cell in the biofilm indicated that only some cells exhibited pulses in electrical activity, with the rest not participating in signaling. This is shown in Fig. 18.12a and b. The distribution of membrane potential amplitudes developed during signal propagation was bimodal, with the fraction of signaling cells being, $\phi =0.43\pm 0.02$, which is in agreement with the percolation threshold obtained from their simulation.

Also measured was the spatial distribution of signaling, or firing, cells that formed clusters of various sizes; see Fig. 18.12c. The distribution of cluster sizes followed a power law with an exponent, $\tau \simeq 2.0$, in excellent agreement with 2D percolation, $\tau =187/91\simeq 2.05$; see Fig. 18.12d. Thus, the spatial organization of firing cells within the bacterial community may be organized near the percolation threshold.

Silva et al. (2019) took advantage of the percolation effect in order to study how to interrupt bacterial communications, both experimentally and with computer simulation. Robustness of the communication networks to interference was studied by analyzing spatial patterns of activation in a model bacterial community that consisted of a sender cell, which produced and responded to autoinducer (AI) signals, and a degrader cell that produced an enzyme capable of destroying the signal. Silva et al. used B. subtilis that secretes the degradative enzyme AiiA. B. subtilis interacts with many soil strains that use acyl-homoserine lactone (AHL)-based quorum sensing. Thus, by varying the ratio of the two strains one can study how they influence the ability of the senders to communicate over long distances.

Silva et al. (2019) mixed sender and degrader cells with a specified ratio and applied it to the top of an LB-agar plate. Cells emit signals that diffuse and spread out via cell division, but are not motile. Sender cells are a strain of Escherichia (E.) coli that contain the LuxRI quorum sensing system and a green fluorescent protein (GFP) reporter gene for quorum sensing activation. As a result of increasing the concentration of the AHL, the signaling molecule increased the expression of the GFP reporter gene, with the response function being nonlinear, exhibiting a sharp increase in cellular fluorescence when the amount of signal exceeded a threshold concentration. Degrader cells, on the other hand, are a strain of E. coli that express constitutively the lactonase enzyme AiiA, which remained inside the degrader cells and degraded AHL quorum sensing signals that diffuse into the cell. The complete details of the experiments are given by Silva et al. (2019).

Experiments by Silva et al. (2019) indicated the existence of a percolation transition that leads to a sudden breakdown in long-range communication at a critical strain ratio. The scaling exponents at the percolation transition point were found to be independent of the specific properties of individual strains. Three percolation exponents, namely, $\beta $, $\gamma $, and $\nu $ that characterize, respectively, the power-law behavior of percolation probability, mean cluster size, and the correlation length near the percolation threshold—see Eqs. (2.8), (2.11), and (2.12)—were estimated based on the experimental data. The results were, $\beta \simeq 0.17\pm 0.06$, $\gamma \simeq 2.28\pm 0.36$, and $\nu \simeq 1.36\pm 0.12$, which are consistent with their exact values in 2D, namely, $\beta =5/36\simeq 0.139$, $\gamma =43/18\simeq 2.39$, and $\nu =4/3\simeq 1.33$.

Notes

1.
The original version of this section, published in the first edition of this book, was written by the late Professor Dietrich Stauffer. The present section represents an updated version of it.
2.
Niels Kaj Jerne (1911–1994) was a Danish immunologist who is known for three important contributions to biological processes. (1) Instead of the body producing antibodies in response to an antigen, Jerne postulated that the immune system already has the specific antibodies it needs to fight antigens. (2) It was known that the immune system learns to be tolerant to the individual’s own self, and Jerne postulated that the learning occurs in the thymus. (3) It was known that the T and B cells communicate with each other, and Jerne’s network theory proposed that the active sites of antibodies are attracted to both specific antigens and to other antibodies that bind to the same site. Jerne shared the Nobel Prize in Physiology or Medicine in 1984 with Georges Jean Franz Köhler (1946–1995) and César Milstein (1927–2002) for “theories concerning the specificity in development and control of the immune system and the discovery of the principle for production of monoclonal antibodies.”
3.
Named after mathematician Richard Wesley Hamming (1915–1998).
4.
George Ralph Mines (1886–1914) was a pioneering English cardiac electrophysiologist who made two fundamental contributions to cardiac electrophysiology. One was his proposal that an action potential wavefront could propagate in a circle, repeatedly activating the tissue, while his second contribution was the discovery of the vulnerable period of the heart—the time during the cardiac cycle when a single stimulus can induce ventricular fibrillation. On 7 November 1914, Mines was found unconscious in his laboratory, and died tragically later that day, cutting short a brilliant career.
5.
Named after Richard FitzHugh who suggested the model (FitzHugh 1961), and Jin-Ichi Nagumo et al. (1962) who developed its equivalent circuit, the model describes a prototype of an excitable system, such as a neuron. It is an example of a relaxation oscillator because, if the external stimulus $I_\textrm{ext}$ exceeds a threshold, the system will exhibit a characteristic excursion in phase space, before the variables relax back to their rest values.

Author information

Authors and Affiliations

Department of Chemical Engineering and Materials Science, University of Southern California, Los Angeles, CA, USA
Muhammad Sahimi

Authors

Muhammad Sahimi
View author publications
You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Muhammad Sahimi .

Rights and permissions

Reprints and permissions

Copyright information

About this chapter

Cite this chapter

Sahimi, M. (2023). Percolation in Biological Systems. In: Applications of Percolation Theory. Applied Mathematical Sciences, vol 213. Springer, Cham. https://doi.org/10.1007/978-3-031-20386-2_18

Download citation

DOI: https://doi.org/10.1007/978-3-031-20386-2_18
Published: 18 March 2023
Publisher Name: Springer, Cham
Print ISBN: 978-3-031-20385-5
Online ISBN: 978-3-031-20386-2
eBook Packages: Mathematics and StatisticsMathematics and Statistics (R0)

Publish with us

Policies and ethics