## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Stochastic cycle selection in active flow networks

Edited by David A. Weitz, Harvard University, Cambridge, MA, and approved June 1, 2016 (received for review March 1, 2016)

## Significance

Nature often uses interlinked networks to transport matter or information. In many cases, the physical or virtual flows between network nodes are actively driven, consuming energy to achieve transport along different links. However, there are currently very few elementary principles known to predict the behavior of this broad class of nonequilibrium systems. Our work develops a generic foundational understanding of mass-conserving active flow networks. By merging previously disparate concepts from mathematical graph theory and physicochemical reaction rate theory, we relate the self-organizing behavior of actively driven flows to the fundamental topological symmetries of the underlying network, resulting in a class of predictive models with applicability across many biological scales.

## Abstract

Active biological flow networks pervade nature and span a wide range of scales, from arterial blood vessels and bronchial mucus transport in humans to bacterial flow through porous media or plasmodial shuttle streaming in slime molds. Despite their ubiquity, little is known about the self-organization principles that govern flow statistics in such nonequilibrium networks. Here we connect concepts from lattice field theory, graph theory, and transition rate theory to understand how topology controls dynamics in a generic model for actively driven flow on a network. Our combined theoretical and numerical analysis identifies symmetry-based rules that make it possible to classify and predict the selection statistics of complex flow cycles from the network topology. The conceptual framework developed here is applicable to a broad class of biological and nonbiological far-from-equilibrium networks, including actively controlled information flows, and establishes a correspondence between active flow networks and generalized ice-type models.

Biological flow networks, such as capillaries (1), leaf veins (2) and slime molds (3), use an evolved topology or active remodeling to achieve near-optimal transport when diffusion is ineffectual or inappropriate (2, 4⇓⇓–7). Even in the absence of explicit matter flux, living systems often involve flow of information currents along physical or virtual links between interacting nodes, as in neural networks (8), biochemical interactions (9), epidemics (10), and traffic flow (11). The ability to vary the flow topology gives network-based dynamics a rich phenomenology distinct from that of equivalent continuum models (12). Identical local rules can invoke dramatically different global dynamical behaviors when node connectivities change from nearest-neighbor interactions to the broad distributions seen in many networks (13⇓⇓–16). Certain classes of interacting networks are now sufficiently well understood to be able to exploit their topology for the control of input–output relations (17, 18), as exemplified by microfluidic logic gates (19, 20). However, when matter or information flow through a noisy network is not merely passive but actively driven by nonequilibrium constituents (3), as in maze-solving slime molds (6), there are no overarching dynamical self-organization principles known. In such an active network, noise and flow may conspire to produce behavior radically different from that of a classical forced network. This raises the general question of how path selection and flow statistics in an active flow network depend on its interaction topology.

Flow networks can be viewed as approximations of a complex physical environment, using nodes and links to model intricate geometric constraints (21⇓–23). These constraints can profoundly affect matter transport (24⇓⇓–27), particularly for active systems (28⇓–30) where geometric confinement can enforce highly ordered collective dynamics (20, 31⇓⇓⇓⇓⇓⇓⇓⇓–40). In symmetric geometries like discs and channels, active flows can often be effectively captured by a single variable

## Model

### Lattice ϕ 6 Field Theory for Active Flow Networks.

Our network is a set of vertices *e*, where *e* in *SI Local Potential*).

Incompressibility, appropriate to dense bacterial suspensions or active liquid crystals, is imposed as follows. The net flux into vertex *v* is *e* is directed out of *v*, *e* is directed into *v*, and 0 if *e* is not incident to *v* (45). Exact incompressibility corresponds to the constraint *λ* and *μ*. This energy is comparable to that of a lattice spin field theory, but with interactions given by higher-dimensional quadratic forms akin to a spin theory on the vertices of a hypergraph (*SI Model*).

### Network Dynamics.

Appealing to recent results showing that bacterial vortex lattices obey equilibrium-like physics (41), we impose that *β* the inverse effective temperature. This stochastic dynamical system has a Boltzmann stationary distribution **2** are**3** arises in an otherwise equivalent fashion to how an elastic energy *V* (*SI Model*).

We now characterize the behavior of this model on a variety of forms of underlying graph *SI Model*), we will focus on connected, simple graphs

## Results

### Stochastic Cycle Selection.

The combination of energy minimization and noise leads to stochastic cycle selection (Movies S1 and S2). A local energy minimum comprises a maximal edge-disjoint union of unit flux cycles: edge fluxes seek to be at *A*–*F* depicts flow on the four-vertex complete graph *A*–*C*) and the generalized Petersen graph *D*–*F*)—the tetrahedron and triangular prism, respectively—where we have integrated Eq. **2** to yield flux–time traces of each edge (*SI Numerical Methods*). The coordinated switching of edges between states of mean flux at *B* and *C*).

A graph possessing an Eulerian cycle—a nonrepeating tour of all edges starting and ending at one vertex—has global energy minima with all edges flowing, which exists if and only if all vertices are of even degree. By contrast, a graph possessing many vertices of odd degree will have minimum energy states with nonflowing edges, because edges flowing into and out of such a vertex pair up to leave an odd number of zero-flow edges. Such “odd” networks are particularly interesting dynamically, as they are more susceptible to noise-induced state switches than graphs with even degree vertices (*SI Model*). For this reason, we specialize from now on to cubic (or 3-regular) graphs where all vertices have degree three.

### Waiting Times and Graph Symmetries.

The cycle-swapping behavior can be quantified by the distribution of the waiting time for an edge to transition between states in **1** are purely topological, with no reference to an embedding of *SI Automorphic Equivalence*); this determines an equivalence relation on *SI Numerical Methods*). The resultant survival function *T* of any edge in *λ* (Fig. 1*G*), and is well approximated by a two-part mixture of exponential distributions.

### Transition Rate Estimation.

Reaction-rate theory explains the form of the waiting time distribution (46). In a system obeying damped noisy Hamiltonian dynamics such as ours, a transition from one local energy minimum to another, respectively *T* for the system to change its state between either minimum is a mixture of two exponential distributions weighted by the equilibrium probabilities of the system to be found in each state. For *H* and *I* shows *λ*, as determined by maximum likelihood estimation on simulation data. Our nonquadratic potential means these rates are not precisely determined by Eq. **4**, but it does suggest an Arrhenius-type dependence *SI Energy Barriers*) and fitting the proportionality constant for each of *H* and *I*).

The complete graph *D*) is vertex transitive, in that any vertex can be permuted to any other by its automorphism group *J*, *Inset*), one containing the two triangles and the other containing the three edges between them. The waiting times then cluster into two distinct distributions according to these two classes. However, when more than two inequivalent minima exist, as they do for *λ* (Fig. 1*J*), consistent with transitions obeying Eq. **4**. But why does one set of edges transition more slowly, on average, than the other? We shall now explore this question for both highly symmetric and totally asymmetric graphs.

### Edge Girth Determines Rate Band Structure.

Global symmetries and local graph structure play distinct roles when determining the transition rates. Fig. 2*A* shows the edge state transition rates for the first eight generalized Petersen graphs *B*), at a representative choice of parameter values that we fix henceforth to focus on the effects of network topology. Inspection of Fig. 2*A* reveals that there are some graphs, such as

When *m* in *ρ* of the edges in *s* running from 0 to 1, as follows. Suppose that only edges in *SI Energy Barriers*). Using the symmetry of *V*, the energy *s* of the transition is given by *H* is maximized precisely when *m*. Therefore, for fixed *ρ*, *m*. This argument suggests that, because the transition rate *e*-girth *e*, so that the usual graph girth is the minimum *e*-girth. Categorizing edge classes in Fig. 2 by *e*-girths.

### Asymmetric Networks.

Even for graphs with no symmetry, the behavior of each edge can still be predicted by a simple local heuristic. For our purposes, a graph with no symmetry is one possessing only the identity automorphism, in which case we say it is asymmetric (45). In this case, edges can have transition rates entirely distinct from one another. Fig. 3*A* depicts the mean transition rates for the edges of 20 asymmetric bridgeless cubic graphs on 21 edges (*SI Numerical Methods*), exemplified in Movie S2. As in Fig. 2, categorizing edges by their *e*-girths (illustrated in Fig. 3*B* for the starred graph; see Fig. S4 for all 20 graphs) splits the rates into near-distinct bands, despite the total absence of symmetry. However, the bands are not perfectly distinct, and high-girth edges, in particular, display a range of transition rates both within and across graphs. A large portion of this variation is accounted for by considering the sizes of all cycles containing an edge. Although the full dependence is highly complex, we can obtain a good transition rate estimate by considering just two cycles. Let *e*. (It may be that *m*-cycle, suppose that flips of these two cycles occur independently with waiting times *A* yields an exponent *C*): the different *e*-girth categories now spread out along the fit line, showing that Eq. **5** yields an easily computed heuristic to estimate the transition rates of edges in a given graph better than

## Discussion

### Incompressible Limit.

Thus far, we have been considering approximate incompressibility with *SI Incompressible Limit*). The most physically intuitive decomposition uses a cycle basis comprising a nonorthogonal set of unit flux cycles, so that each basis element corresponds to adding or removing a unit of flux around a single cycle. Finite planar graphs, in particular, possess a highly intuitive cycle basis. Fix a planar embedding for *e* borders face *α*, and is then

### Example.

Fig. 4 shows an integration of Eq. **6** for an embedding of the graph *SI Incompressible Limit* and Fig. S2). [In fact, the dual of a polyhedral graph is unique (48).] Note that the *difference* *α* and *β* must be near *d* to the external face, which is constrained to zero flux, can be metastable at values up to

### Low Temperature Limit and Ice-Type Models.

Similar to how a lattice

### Complex Networks.

We have focused on small regular graphs, but the dynamical principles presented here will still apply to active flow on complex networks. The edges in a large random graph typically exhibit a wide distribution of *e*-girths, where topologically protected edges, whose *e*-girth is large enough to prevent them ever changing state within a realistic observation window, coexist with frequently switching edges of small *e*-girth. In fact, graphs drawn from distributions modeling real-life network phenomena (13, 14) seem to have far more small *e*-girth edges than their fixed degree or uniformly random counterparts (*SI Complex Networks* and Fig. S3). Furthermore, although large random graphs are almost always asymmetric (45), many real-life complex networks have very large automorphism groups (51), meaning that, as in Fig. 2, there will be large sets of edges in such a network with identical transition rates. Active flow on complex networks can therefore be expected to display a rich phenomenology of local and global state transitions.

## Conclusions

Our analysis shows that the state transition statistics of actively driven quasi-incompressible flow networks can be understood by combining reaction rate theory with graph-theoretic symmetry considerations. Furthermore, our results suggest that nonequilibrium flow networks may offer new insights into ice-type models and vice versa. The framework developed here offers ample opportunity for future generalizations from both a biophysical and a transport optimization perspective. For example, an interesting open biological question concerns how plasmodial organisms such as *Physarum* (3, 6, 7) adapt and optimize their network structure in response to external stimuli, such as light or nutrient sources or geometric constraints (52, 53). Our investigation suggests that a combined experimental and mathematical analysis of cycle structure may help explain the decentralized computation strategies used by these organisms. More generally, it will be interesting to explore whether similar symmetry-based statistical approaches can guide the topological optimization of other classes of nonequilibrium networks, including neuronal and man-made information flow networks that typically operate far from equilibrium.

## Materials and Methods

Eqs. **2** and **6** were integrated by the Euler–Maruyama method (54) with time step *SI Numerical Methods*. All data are available on request.

In this document, we cover details of statements and derivations in the main text. We first elaborate on the need for a local potential of order higher than quartic, before discussing details of our model: analogies to spin field theories and pure diffusive models, decoupling of loops, and the more active dynamics of a graph of odd degree versus one of even degree. We next show that edges that are equivalent under automorphism obey identical waiting time distributions, we compute the energy barriers for the 3-cycle/4-cycle transitions in *e*-girth distributions to predict the behavior of active flows on complex networks. We also provide details on numerical methods used throughout this work.

In addition to figures referenced in the following text, we also provide two further figures: Fig. S4, which supplements Fig. 3 by illustrating all 20 asymmetric cubic graphs considered; and Fig. S5, which shows an example integration of totally incompressible flow on a

## SI Local Potential

The typical symmetric, bistable potential is the quartic

Consider the elementary (although not simple!) two-vertex, three-edge graph in Fig. S1*A*. This has energy*λ*.

Consider the case

## SI Model

### Relationship to Lattice Spin Field Theory.

Our energy in Eq. **1** can be seen as a generalization of a lattice spin field theory. Suppose we switch to a typical vertex-based picture, where fluxes *e* in *i* in an interaction graph *i* and *j* is antiferromagnetic or ferromagnetic, respectively. In our theory, however, multiple spins are permitted inside each interaction term according to the degree of each vertex in **1** is equivalent to**S1** being the special case where **S1** has two types of interaction edge—antiferromagnetic and ferromagnetic—between two spins, the general theory has *n* spins for all

### Diffusive Dynamics When λ = 0 .

In the absence of a local edge potential (i.e., **2** reduces to noisy scalar diffusion on the edges of

As in the derivation of the incompressible limit (see *SI Incompressible Limit* and *Incompressible Limit*), by analogy with a spectral decomposition for the diffusion equation, we decompose *μ* whereas modes with

### Decoupling of Loops.

We show here that if *w*, because flow in along

### Odd- Versus Even-Degree Vertices.

Graphs with odd-degree vertices exhibit more active stochastic cycle selection dynamics than those with even-degree vertices. This is exemplified by the small graphs in Fig. S1, where adding an extra edge markedly slows transition rates. For the graph in Fig. S1*A* to change state while conserving flux, one edge changes from *B*, one edge changes from *A*.

## SI Automorphic Equivalence

In this section, we show that edges within an automorphism equivalence class obey identical waiting time distributions. As in the main text, we assume that

To permit multiple edges, we define an automorphism **2**, whose components read*σ*, so that *e* with **S2** and substituting this definition implies*σ* is a permutation, we can reorder the sums as*σ* preserves incidence but not necessarily orientation, *e* with respect to *v*. Therefore, Eq. **S3** becomes**S4** by *H*. Therefore, any edges

Note that, in the incompressible limit *σ* exists, even in connected simple graphs.

## SI Energy Barriers

We describe here the process of computing the transition energy barriers for *A*, let

Note that the transition is not exactly of the form

## SI Incompressible Limit

### Dimensional Reduction.

Let *i*. As *Incompressible Limit*.

### Cycle Basis for P 4,1 .

We detail here the derivation of the planar cycle basis representation in the incompressible limit for the cube *A*, orient and number the edges as indicated. Next, construct the dual of the (undirected) plane graph, with vertices numbered as in Fig. S2*B*, and assign an orientation to each edge of the dual such that *A* and *B* is *α* (Fig. S2*B*). This has incidence matrix*B*—reads**6** reading

## SI Complex Networks

Finally, we elaborate on discussions in the main text extrapolating our results to expected behavior on large complex networks.

Unlike the cubic graphs we have focused on, a complex network possesses a broad vertex degree distribution. This will certainly affect transition rates, because the presence of even-degree vertices deepens energy minima (Fig. S1). However, because the effect of vertex degree is broadly independent of cycle structure, we predict that the distribution of transition rates will still qualitatively match the *e*-girth distribution. In Fig. S3, we plot *e*-girth distributions derived from 10 instances each of four random graph distributions on 1,000 vertices: fixed degree 3, uniform, Barabási–Albert (“scale free”), and Watts–Strogatz (“small world”). Of the four, random cubic graphs display by far the highest *e*-girths, whereas the Barabási–Albert and Watts–Strogatz graphs, commonly used as prototypes of certain forms of real-life complex networks, both retain many edges with low *e*-girth despite their size. Therefore, by this measure, complex networks may exhibit a far greater proportion of fast-switching edges than random cubic graphs on the same number of vertices.

## SI Numerical Methods

### Numerical Integration and Waiting Times.

Eqs. **2** and **6** were integrated by the Euler–Maruyama method with time step *λ* in Fig. 1, and over 24 (Fig. 2) or 8 (Fig. 3) integrations to

### Graph Generation and Properties.

Mathematica (Wolfram Research, Inc.) was used to generate graphs and their incidence matrices, and to determine all graph-theoretic properties, including cycle lengths and edge equivalence classes. The graphs in Fig. 3 (see also Fig. S4) were chosen uniformly at random from the database of all nonisomorphic bridgeless connected cubic graphs on 14 vertices accessible in Mathematica after filtering to discard those with nontrivial automorphism group.

## Acknowledgments

We thank R. Goldstein, I. Pivotto, and G. Royle for discussions. F.G.W. was supported by Trinity College, Cambridge; A.F. and J.D. were supported by NSF Award CBET-1510768; J.B.F. was supported by ARC Discovery Grant DP130100106; and J.D. was supported by an MIT Solomon Buchsbaum Fund Award and an Alfred P. Sloan Research Fellowship.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: F.G.Woodhouse{at}damtp.cam.ac.uk.

Author contributions: F.G.W., A.F., J.B.F., and J.D. designed research, performed research, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1603351113/-/DCSupplemental.

## References

- ↵
- ↵
- ↵.
- Alim K,
- Amselem G,
- Peaudecerf F,
- Brenner MP,
- Pringle A

*Physarum polycephalum*organizes fluid flows across an individual. Proc Natl Acad Sci USA 110(33):13306–13311. - ↵
- ↵
- ↵
- ↵.
- Tero A, et al.

- ↵
- ↵
- ↵
- ↵.
- Garavello M,
- Piccoli B

- ↵
- ↵
- ↵.
- Barabási AL,
- Albert R

- ↵
- ↵.
- Bianconi G,
- Pin P,
- Marsili M

- ↵
- ↵
- ↵.
- Prakash M,
- Gershenfeld N

- ↵.
- Pearce DJ,
- Turner MS

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Lushi E,
- Wioland H,
- Goldstein RE

- ↵
- ↵
- ↵.
- Buhl J, et al.

- ↵.
- Yates CA, et al.

- ↵
- ↵
- ↵.
- Tjhung E,
- Marenduzzo D,
- Cates ME

- ↵
- ↵.
- Baxter RJ

- ↵
- ↵
- ↵.
- Godsil C,
- Royle GF

- ↵
- ↵
- ↵
- ↵.
- Grimmett G

- ↵
- ↵
- ↵.
- Reid CR,
- Latty T,
- Dussutour A,
- Beekman M

- ↵
- ↵

## Citation Manager Formats

### More Articles of This Classification

### Biological Sciences

### Biophysics and Computational Biology

### Physical Sciences

### Related Content

- No related articles found.

### Cited by...

- No citing articles found.