Mode selection in compressible active flow networks
MMode selection in compressible active flow networks
Aden Forrow, Francis G. Woodhouse, and J¨orn Dunkel ∗ Department of Mathematics, Massachusetts Institute of Technology,77 Massachusetts Avenue, Cambridge MA 02139-4307, U.S.A. Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, U.K. (Dated: September 12, 2018)Coherent, large scale dynamics in many nonequilibrium physical, biological, or information trans-port networks are driven by small-scale local energy input. Here, we introduce and explore an ana-lytically tractable nonlinear model for compressible active flow networks. In contrast to thermally-driven systems, we find that active friction selects discrete states with a limited number of oscillationmodes activated at distinct fixed amplitudes. Using perturbation theory, we systematically predictthe stationary states of noisy networks and find good agreement with a Bayesian state estimationbased on a hidden Markov model applied to simulated time series data. Our results suggest that themacroscopic response of active network structures, from actomyosin force networks to cytoplasmicflows, can be dominated by a significantly reduced number of modes, in contrast to energy equipar-tition in thermal equilibrium. The model is also well-suited to study topological sound modes andspectral band gaps in active matter.
PACS numbers: 47.63.-b, 05.70.Ln, 05.65.+b 05.40.-a
Active networks constitute an important class ofnonequilibrium systems spanning a wide range of scales,from the intracellular cytoskeleton [1, 2] and amoeboidorganisms [3–8] to macroscopic transport networks [9–12]. Identifying generic self-organization principles [13–15] that control the dynamics of these biological or ar-tificial far-from-equilibrium systems remains one of theforemost challenges of modern statistical physics. De-spite promising experimental [4–6, 16–18] and theoret-ical [1, 7, 19–21] advances over the past decade, it isnot well understood how the interactions between localenergy input, dissipation and network topology deter-mine the coordinated global behaviors of cells [16], plas-modia [4–6] or tissues [22]. Further progress requires ana-lytically tractable models that help clarify the underlyingnonequilibrium mode selection principles [23, 24].We inroduce here a generic model for active flows ona network, motivated by recent experimental studies ofbacterial fluids [20, 25] and ATP-driven microtubule sus-pensions [26] in microfluidic channel systems. Buildingon Rayleigh’s work [27] on driven vibrations and theToner-Tu model of flocking [28], the theory accountsfor network activity through a nonlinear friction [28–31]. We work in a fully compressible framework allow-ing accumulated matter at vertices to affect flow throughnetwork pressure gradients, generalizing previous workon incompressible pseudo-equilibrium active flow net-works [32, 33], as suited to the many biological systemsexhibiting flexible network geometry [4–6] or variationsin the density of active components [15]. Although in-herently nonlinear, the model can be systematically ana-lyzed through perturbation theory. Such analysis showshow slow global dynamics emerge naturally from the fastlocal dynamics, enabling prediction of the typical statesin large noisy networks; these states have significantly fewer active modes than for energy equipartition [34] inthermal equilibrium. More broadly, our model providesan accessible framework for investigating generic physi-cal phenomena in active systems, including topologically-protected sound modes [15] and the influence of spectralband gaps (SM [35]).We consider activity-driven mass flow on anarbitrarily-oriented graph G = ( V , E ) with V = | V | ver-tices and E = | E | edges. The elements of the V × E gradient (incidence) matrix ∇ are ∇ ve = − e is oriented outwards from vertex v , ∇ ve = +1 if e isoriented inwards into v , and ∇ ve = 0 otherwise. The dy-namical state variables are the deviations from the meanmass ¯ (cid:37) = M/V on the nodes, ( (cid:37) ( t ) , . . . , (cid:37) V ( t )), and themass fluxes on the edges, ( φ ( t ) , . . . , φ E ( t )), governed bythe non-dimensionalized (SM [35]) transport equations˙ (cid:37) v = (cid:88) e ∇ ve φ e , (1a)˙ φ e = − (cid:88) v ∇ (cid:62) ev (cid:37) v + (cid:15) µ − φ e φ e φ e + √ Dξ e ( t ) , (1b)where ξ e ( t ) is standard Gaussian white noise. Equa-tion (4a) ensures mass conservation. The first termon the r.h.s. of Eq. (4b) represents the gradient of anideal gas-type node pressure p v ∝ (cid:37) v , corresponding tothe leading term in a virial expansion; the second termis a Toner-Tu type (SM [35]) active friction force de-rived from a depot model [29, 36] with coupling (cid:15) > µ , which drivesthe edge fluxes φ e towards preferred values ±√ µ when µ >
0. Many networks have non-uniform edge and ver-tex weights, which can be incorporated into equations ofidentical form to Eqs. (4) with appropriate rescaling of (cid:37) , φ , and ∇ (SM [35]). a r X i v : . [ phy s i c s . b i o - ph ] J u l FIG. 1. Activity can select a single dominant oscillationmode on hierarchically weighted networks. (a) The edges inthe graph simulated in (b) and (c) are given weights decreas-ing exponentially with their distance from the central redpath. (b) Oscillations in pressure and flux develop primar-ily along the central high-weight path (Movie 1). (c) Edgefluxes φ e settle into steady synchronized oscillations as exem-plified for two edges indicated in (b), one on ( φ ) and one off( φ ) the path. (d) Plotting the time-dependent amplitudeof each analytically-determined flow eigenmode confirms se-lection of a single oscillatory mode. The ten modes with thehighest average amplitude in this simulation run are pictured;the marked top two rows are oscillatory modes, while the re-maining rows are cyclic modes. See Fig. S6 for all modes.Simulation parameters were (cid:15) = 0 . µ = 1, and D = 10 − . Active flow networks described by Eqs. (4) exhibit richoscillatory transport behavior, including the mode selec-tion illustrated in Movie 1 and Fig. 1 for a hierarchically-weighted network with vertex degrees at most 3 as is typ-ical of
Physarum polycephalum [37]. When this networkis initialized with zero pressure variation and flux, it typ-ically settles into a quasi-steady state with a single dom-inant oscillation frequency on the highest-weight path.This is a manifestation of the fact that single-frequencyselection is the norm on actively driven path graphs, aswe shall show analytically below.Generally, the features of the steady-state attractorwill be determined by the topology of the subgraph ofhigh-weight edges, which may be much sparser than theoriginal network. For this reason, as well as for easeof analysis and illustration, we will henceforth assume G to be a tree, as realized in certain peripheral sensoryneurons [38], though in general the full model in Eqs. (4)is not restricted to any particular class of graph. Thebehaviors observed on trees can be extended to denser graphs by choosing appropriate edge weights.The complex active flow dynamics encoded by Eqs. (4)can be understood analytically by considering the basisof oscillation modes of the network, as we illustrate nowin the fully deterministic case ( D = 0). To progress,we adopt a Rayleigh [27] approximation (cid:15) ( µ − φ e ) φ e forthe active friction (SM [35]). Now, expand the pressure (cid:37) v = (cid:80) En =1 r n ( t ) (cid:37) vn and flux φ e = (cid:80) En =1 f n ( t ) φ en in theright and left singular vectors (cid:37) n = ( (cid:37) vn ) and φ n = ( φ en )of ∇ (cid:62) corresponding to the E = V − λ n . (On a tree, there is a single zero eigenvalueof ∇∇ (cid:62) yielding an additional right singular vector forthe pressure, but this corresponds to a constant massshift and so can be safely neglected.) Defining mode am-plitudes A n = r n + f n , the network energy then takesthe simple form H = (cid:80) n λ n A n (SM [35]). When (cid:15) is small there are two distinct timescales, namely thefast oscillation timescale t and the slow friction timescale τ = (cid:15)t , which we separate in the perturbation ansatz r n = (cid:80) ∞ σ =0 (cid:15) σ r σn and f n = (cid:80) ∞ σ =0 (cid:15) σ f σn [39]. Activefriction does not contribute at lowest order, so the O (1)contribution to each mode ( r n , f n ) is an uncoupled har-monic oscillator r n ( t ) = A n ( τ ) cos[ λ n t − δ n ( τ )] and f n ( t ) = − A n ( τ ) sin[ λ n t − δ n ( τ )] with t -independent am-plitude A n and phase δ n (SM [35]).The influence of activity becomes apparent at first or-der in (cid:15) , introducing couplings between mode amplitudeswhose dynamics encode the state selection behavior ofthe active network. Requiring that the O ( (cid:15) ) amplitudes r n and f n remain small relative to the leading termsimplies that the secular (unbounded) terms in the firstorder equations must vanish [39]. Assuming negligiblemode degeneracies, the slow dynamics of the O (1) modeamplitudes A n ( τ ) are found to obey (SM [35]) d ( A n ) dτ = (cid:32) µ − E (cid:88) k =1 P nk A k (cid:33) A n , (2)where the overlap matrix P nk = (1 − δ nk ) (cid:80) e φ en φ ek encodes the network topology. Fixed points of Eq. (16)can then be found by choosing a subset of the A n tobe zero and solving (cid:80) Ek =1 P nk A k = µ for A n over theremaining non-zero modes. If all the non-zero solutionsfor A n are positive, then there is a stationary point withthose modes activated (SM [35]).Activity-driven fixed points with exactly one mode ac-tive always exist. If only mode p is active at leading or-der, then A n = (cid:112) µ/P pp δ np is a fixed point of Eq. (16).These amplitudes, which closely match both those calcu-lated with the full unapproximated active friction forceand those from averages computed over fully nonlinearsimulations (SM [35]), show that as µ crosses 0 there isa supercritical Hopf bifurcation with A n ∼ √ µ . How-ever, the stability of such a single-mode state depends ontopology: our simulations suggest that activity always se-lects exactly one oscillation mode in simple path graphs, Mode amplitude A M od e a m p lit ud e A Mode amplitude A M od e a m p lit ud e A Mode amplitude A M od e a m p lit ud e A T i m e A A A A E dg e f l ux V e r t e x p r e ss u r e M od e a m p lit ud e
000 200 200 200 2000Time Time Time (a) (b)(c)
FIG. 2. First order perturbation theory accurately predicts the stable states on small trees. (a) A five vertex tree possessingfour nontrivial modes, as illustrated. (b) On the tree in (a), mode amplitudes settle into one of two stable stationary states,as seen in simulations for three different initial conditions. Modes are ordered by frequency from high (top) to low (bottom).(c) Simulated mode trajectories (rainbow) in (b) match analytic predictions (blue streamlines) in the subspaces of activatedmodes. There are three possible arrangements of nonzero critical points in each 2D subspace: a saddle point on one axis anda stable node on the other axis (left), a stable node on each axis and a saddle point in the middle (center), or a saddle pointon each axis and a stable node in the middle (right; Movie 2). Higher order effects cause both the convergence to a point with A > (cid:15) = 0 . µ = 1, D = 0. whereas single-mode states are typically unstable in net-works with complex topologies. We can use this observa-tion to model more complex active networks with singlemode selection by appropriately weighting the edges: ifthe edge weights for a path are large enough compared tothe weights elsewhere in the network, the path behaviordominates (Fig. 1).Insight into stability is provided by the case with upto two modes active. Writing A n = A p δ np + A q δ nq ,Eq. (16) yields d ( A p ) /dτ = ( µ − P pp A p − P pq A q ) A p , (3)and symmetrically for A q . Depending on the topology-encoding overlap coefficients P nk , this gives up tofour fixed points: the zero state A p = A q = 0,which is always linearly unstable; the single-mode state( A p , A q ) = ( (cid:112) µ/P pp , P pq > P pp and a saddle if not, plus analogously for (0 , (cid:112) µ/P qq );and, potentially, a mixed state ( A ∗ p , A ∗ q ) where A ∗ p = (cid:113) µ ( P qq − P pq )/ (cid:0) P pp P qq − P pq (cid:1) with A ∗ q defined sym-metrically. When it exists, the mixed state is either sta-ble (if P pq < P pp P qq ) or a saddle (if P pq > P pp P qq ), butif one of the single-mode states is stable and one is un-stable, then one of A ∗ p and A ∗ q is imaginary and thereis no mixed state. Hence, we have three possible sce-narios (Fig. 2): one stable single mode and the othera saddle with no mixed state (Fig. 2b,c; left); two sta-ble single-mode states with a mixed saddle in-between(Fig. 2b,c; center); and two single-mode saddles with astable mixed state in-between (Fig. 2b,c; right). Thesepredictions match simulations quantitatively even for rel-atively large (cid:15) beyond the small- (cid:15) perturbation regime (Fig. 2). In fact, simulations show the same qualitativebehavior for (cid:15) = 2, suggesting perturbation analysis re-mains predictive at high activity.This two-mode analysis yields a simple topologicalheuristic for the stability of single-mode states. Since | φ p | = 1, P pp is small when φ p is spread over many edgesand large when φ p is localized to a few edges. If φ q islocalized to the same edges as φ p , P pq will also be largeand mode p will be stable to perturbations in mode q .However, if φ q is localized to a disjoint set of edges, P pq will be a scaled inner product of near-orthogonal vectors( φ ep ) and ( φ eq ) and will be small. Thus localized modeswill be unstable to modes in other regions, while con-versely if a mode is to be stable alone then it will bespread out across the entire network. Therefore, a stablecombination of modes will possess significant flows on alledges of the network.Biological systems exhibit vastly different macroscopicand microscopic time scales [40–43]. This phenomenonis present in our compressible active flow network, wherehigher-order nonlinear effects induce slow global timescales from faster small-scale dynamics. When thezeroth-order amplitudes A n are at a fixed point, thefirst-order corrections r n and f n are harmonic oscilla-tors with natural frequency λ n driven at linear combina-tions of the frequencies active at zeroth order (SM [35]).For instance, if two modes p and q are active at zerothorder, the driving frequencies are 3 λ p − k ( λ p ± λ q ) for k = 0 , . . . ,
3. This introduces new, slower timescalesinto the dynamics, including oscillations in the energy H = (cid:80) n λ n ( r n + f n ) with frequency λ p − λ q . Theirmagnitude depends on the difference in frequency: slower FIG. 3. States on larger trees possess surprisingly few active modes, which can be inferred from time series with non-zero noise.(a) The mean number of stationary states of Eq. (16) grows exponentially with edges E as 1 . E ≈ (2 E ) / (solid orange line),close to the upper bound of 2 E states (dashed black line), while the mean number of stable states grows as 1 . E ≈ (2 E ) / (solid blue line). We counted states on all nonisomorphic trees with E ≤
14 edges (filled circles) and on a random sample of ∼
175 trees per point for 15 ≤ E ≤
24 (open circles). Averages are over trees with a fixed number of edges. (b) As E increases,both the mean and the variance of the distribution of trees with each number of stable states increase rapidly. (c) Distributionof the average number of modes active in a stable state. The mean over trees scales like 0 . E ≈ E/ E/ × ) always only activate one mode; complex trees (+) have more modes active. (e) Noisynetworks ( D >
0) transition stochastically between stable states, exemplified by an amplitude-time trace for the tree shown.Modes are ordered by frequency from high (top) to low (bottom). Simulation parameters are (cid:15) = 0 . µ = 1, D = 5 × − .(f) States found by vbFRET from simulations on the tree in (e) (SM [35]). The second, first, and fifth columns are states seenin (e), indicated by the colored bars above. (g) States predicted by Eq. (16) for the tree in (e). The first five states in (f) matchthose in (g); the sixth column in (f) is likely a transient combination of analytically stable states. oscillations, driven by modes with similar frequencies λ p ≈ λ q , have higher amplitudes (SM [35], Fig. S7).The number of activated modes in an arbitrary com-pressible active network depends on intricate interac-tions between local activity and global flow configura-tions. The total number of available modes is equal tothe number of edges E , meaning that, were each com-bination of modes to be a fixed point, a tree could haveup to 2 E stationary states. To see how the true num-ber of stationary and stable states depends on tree size,we performed an exhaustive numerical fixed point searchof Eq. (16) over a large sample of trees with E ≤ E suggests expo-nential growth of the mean number of steady states withedges E ; this is indeed what we see, going as ∼ (2 E ) / .However, though still exponential in E , the mean numberof stable states is much smaller at ∼ (2 E ) / (Fig. 3a).Remarkably, these stable states have only ∼ E/ E modes under thermal equipartition [34].Path-like topologies lead to even more dramatic reduc-tions in the number of modes active (Fig. 3c), suggestingthat a biological system can further reduce the number of active modes through an optimal choice of topology;moreover, hierarchically tuned edge capacities as realizedin Physarum [5, 6, 37] can further enhance mode selectioneven in non-tree topologies (Fig. 1).Real active transport networks will have some nonzerolevel of thermal or athermal noise [44–46]. Provided thenoise is not too large, it will render previously stablestates now only metastable, with flow patterns exhibitingsmall fluctuations around these metastable states punc-tuated by noise-driven stochastic transitions betweenthem [32, 46]. Long-time simulations of Eqs. (4) with
D > ∗ [email protected][1] C. Broedersz and F. MacKintosh, Rev. Mod. Phys. ,995 (2014).[2] V. Ruprecht, S. Wieser, A. Callan-Jones, M. Smutny,H. Morita, K. Sako, V. Barone, M. Ritsch-Marte, M. Sixt,R. Voituriez, and C.-P. Heisenberg, Cell , 673 (2015).[3] A. Takamatsu, R. Tanaka, H. Yamada, T. Nakagaki,T. Fujii, and I. Endo, Phys. Rev. Lett. , 7 (2001).[4] A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber,M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Naka-gaki, Science , 439 (2010).[5] K. Alim, G. Amselem, M. P. Brenner, and A. Pringle,Proc. Natl. Acad. Sci. U.S.A. , 13306 (2013).[6] K. Alim, N. Andrew, A. Pringle, and M. P. Brenner,Proc. Natl. Acad. Sci. USA , 5136 (2017).[7] V. Bonifaci, K. Mehlhorn, and G. Varma, J. Theor. Biol. , 121 (2012).[8] C. R. Reid, H. Macdonald, R. P. Mann, J. A. R. Marshall,T. Latty, and S. Garnier, J. R. Soc. Interface , 44(2016).[9] G. Coclite, M. Garavello, and B. Piccoli, SIAM J. Math.Anal. , 1862 (2005).[10] B. Piccoli and M. Garavello, Traffic Flow on Networks:Conservation Laws Models (AIMS, Springfield, MO,2006).[11] S. Hata, H. Nakao, and A. S. Mikhailov, Phys. Rev. E , 020801 (2014).[12] L. L. Heaton, E. L´opez, P. K. Maini, M. D. Fricker, andN. S. Jones, Phys. Rev. E , 021905 (2012).[13] H. Nakao and A. S. Mikhailov, Nat. Phys. , 544 (2010).[14] A. S. Mikhailov and K. Showalter, Chaos , 026101(2008).[15] A. Souslov, B. C. van Zuiden, D. Bartolo, and V. Vitelli, arXiv:1610.06873.[16] N. Fakhri, A. D. Wessel, C. Willms, M. Pasquali, D. R.Klopfenstein, F. C. MacKintosh, and C. F. Schmidt,Science , 1031 (2014).[17] P. Ronceray, C. P. Broedersz, and M. Lenz, Proc. Natl.Acad. Sci. U.S.A. , 2827 (2016).[18] S. Marbach, K. Alim, N. Andrew, A. Pringle, and M. P.Brenner, Phys. Rev. Lett. (2016).[19] C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C.MacKintosh, Nat. Phys. , 983 (2011).[20] M. Paoluzzi, R. Di Leonardo, and L. Angelani, Phys.Rev. Lett. , 188303 (2015).[21] A. Bressan, S. Cani´c, M. Garavello, M. Herty, andB. Piccoli, EMS Surv. Math. Sci. , 47 (2014).[22] C. G. Vasquez and A. C. Martin, Dev. Dynam. , 361(2016).[23] W. Ebeling, U. Erdmann, J. Dunkel, and M. Jenssen, J.Stat. Phys. , 443 (2000).[24] J. Dunkel, W. Ebeling, U. Erdmann, and V. A. Makarov,Int. J. Bifurcat. Chaos , 2359 (2002).[25] H. Wioland, E. Lushi, and R. E. Goldstein, New J. Phys. , 075002 (2016).[26] K.-T. Wu, J. B. Hishamunda, D. T. N. Chen, S. J. De-Camp, Y.-W. Chang, A. Fern´andez-Nieves, S. Fraden,and Z. Dogic, Science , eaal1979 (2017).[27] J. W. S. B. Rayleigh, The Theory of Sound vol. 1 , 2nded. (Macmillan, New York, 1894) p. 81.[28] J. Toner, Y. Tu, and S. Ramaswamy, Annals of Physics , 170 (2005).[29] F. Schweitzer, W. Ebeling, and B. Tilch, Phys. Rev.Lett. , 5044 (1998).[30] P. Romanczuk, M. B¨ar, W. Ebeling, B. Lindner, andL. Schimansky-Geier, Eur. Phys. J. Spec. Top. , 1(2012).[31] P. S. Burada and B. Lindner, Phys. Rev. E , 032102(2013).[32] F. G. Woodhouse, A. Forrow, J. B. Fawcett, andJ. Dunkel, Proc. Natl. Acad. Sci. U.S.A. , 8200(2016).[33] F. G. Woodhouse and J. Dunkel, Nat. Commun. , 15169(2017).[34] A. I. Khinchin, Mathematical Foundations of StatisticalMechanics (Dover, New York, 1949).[35] See Supplemental Material, which includes Refs. [48].[36] P. Romanczuk, W. Ebeling, U. Erdmann, andL. Schimansky-Geier, Chaos , 047517 (2011).[37] W. Baumgarten, T. Ueda, and M. J. B. Hauser, Phys.Rev. E , 046113 (2010).[38] J. Kromer, A. Khaledi-Nasab, L. Schimansky-Geier, andA. B. Neiman, ArXiv:1701.01693.[39] S. H. Strogatz, Nonlinear Dynamics and Chaos (West-view Press, Boulder, CO, 2015).[40] J. Halatek and E. Frey, Cell Rep. , 741 (2012).[41] R. A. Kerr, H. Levine, T. J. Sejnowski, and W. Rappel,Proc. Natl. Acad. Sci. U.S.A. , 347 (2006).[42] I. H. Riedel-Kruse, C. M¨uller, and A. C. Oates, Science , 1911 (2007).[43] A. Varma, K. C. Huang, and K. D. Young, J. Bacteriol. , 2106 (2008).[44] B. Lindner, J. Garcıa-Ojalvo, A. Neiman, andL. Schimansky-Geier, Phys. Rep. , 321 (2004).[45] L. Gammaitoni, P. H¨anggi, P. Jung, and F. Marchesoni,Rev. Mod. Phys. , 223 (1998).[46] P. H¨anggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. , 251 (1990).[47] J. E. Bronson, J. Fei, J. M. Hofman, R. L. Gonzalez, andC. H. Wiggins, Biophys. J. , 3196 (2009).[48] J. W. S. B. Rayleigh, Proc. Lond. Math. Soc. s1-10 , 4(1878); P. Misra, Physics of Condensed Matter (ElsevierScience, 2011).
Supplemental material: Mode selection in compressible active flow networks
Aden Forrow, Francis G. Woodhouse, and J¨orn Dunkel
Nondimensionalization of governing equations
We can define the model in terms of the dimensional quantities ˆ (cid:37) v , ˆ φ e , and ˆ t ; global dimensional parameters ˆ (cid:15) ,ˆ β , and ˆ D ; dimensionless edge conductances γ e and vertex volumes m v ; and a dimensionless global parameter µ andfunction g as d ˆ (cid:37) v d ˆ t = (cid:88) e ∇ ve ˆ φ e ,d ˆ φ e d ˆ t = − ˆ γ e (cid:88) v ∇ (cid:62) ev m − v ˆ (cid:37) v + ˆ (cid:15)g (cid:18) µ, ˆ φ e ˆ β ˆ γ e (cid:19) ˆ φ e + (cid:112) D ˆ ξ e (ˆ t ) . The scaling by conductance in the argument of g is chosen to match the phenomenology observed in dense bacterialsuspensions, where activity selects a characteristic velocity φ e /γ e and not a fixed flux φ e . If we choose a conductancescale ˆ γ and volume scale ˆ m and insert the rescaled, nondimensional parameters γ e = ˆ γ − ˆ γ e , m v = ˆ m − ˆ m v , (cid:15) = ˆ γ − ˆ (cid:15), D e = ˆ β − ˆ γ − γ − e ˆ D and variables (cid:37) v = m v ˆ γ ˆ β − ˆ (cid:37) v , φ e = γ − e ˆ β − ˆ φ e , t = ˆ γ ˆ t, ξ e ( t ) = ˆ γ − ˆ ξ e (ˆ t ) , we are left with d(cid:37) v dt = (cid:88) e m − / v ∇ ve γ / e φ e ,dφ e dt = − (cid:88) v γ / e ∇ (cid:62) ev m − / v (cid:37) v + (cid:15)g (cid:18) µ, φ e √ γ e (cid:19) φ e + (cid:112) D e ξ e ( t ) . With constant conductances γ e = 1 and volumes m v = 1, we recover the model introduced in the main text, namely d(cid:37) v dt = (cid:88) e ∇ ve φ e , (4a) dφ e dt = − (cid:88) v ∇ (cid:62) ev (cid:37) v + (cid:15)g ( µ, φ e ) φ e + √ Dξ e ( t ) , (4b)with nonzero entries of the gradient matrix equal to ±
1. All of our analysis applies equally well to the varying weightscase: the only substantive change is replacing ∇ ve with the weighted gradient ∇ ∗ ve = m − / v ∇ ve γ / e .We can combine Eqs. (4a) and (4b) into one second order equation for the pressure dynamics reading¨ (cid:37) v = (cid:88) e ∇ ve (cid:32) − (cid:88) u ∇ (cid:62) eu (cid:37) u + (cid:15)g ( µ, φ e ) φ e + √ Dξ e ( t ) (cid:33) . (5)In the absence of friction, when g ( µ, φ e ) = 0, the dynamics are Hamiltonian with energy H = 12 (cid:88) v,e,u (cid:37) v ∇ ve ∇ (cid:62) eu ρ u + 12 (cid:88) e,v,f φ e ∇ (cid:62) ev ∇ vf φ f . (6)The energy is particularly simple when written in the basis of singular vectors of ∇ (cid:62) with non-zero singular values,giving H = 12 (cid:88) n λ n (cid:0) r n + f n (cid:1) ≡ (cid:88) n H n . Relation to physical flow systems
We chose to explore a minimal model coupling local active energy input to network structure, rather than capturethe details of any particular model system. Nevertheless, the key features of our model, namely mass conservationand a polynomial expansion of the active term, are generic enough to be straightforwardly adapted to a range ofapplications.Mass conservation and pressure driven flow are likely to remain in any active flow model; the form of the activeterm may change in different contexts. In our case, staying close to examples of bacterial suspensions, we modelactivity as driving spontaneous flow on all edges. An alternative option, more closely related to shuttle streaming innetworks, would be to apply an active force f v that compresses or expands each vertex and drives flow in or out, withmodified dynamics d(cid:37) v dt = (cid:88) e ∇ ve φ e ,dφ e dt = − (cid:88) v ∇ (cid:62) ev ( (cid:37) v + (cid:15)f v ) + √ Dξ e ( t ) . The correct form of the active force depends on the microscopic details of the driving. Some generic features, however,will not depend on the exact form of f v and will be discoverable by choosing a simple function of local quantities( (cid:37) v , ˙ (cid:37) v , etc.) as an approximate driving force.The same method is used to derive the Toner-Tu equations for continuous active flows [28]; our model can beunderstood as a discrete version of a special case of these equations. If advective and diffusive terms are renderednegligible in favor of pressure-driven and activity-driven flow by geometric effects or otherwise, and we take only thelinear term in the virial expansion of the active pressure, the general Toner–Tu model simplifies to ∂(cid:126)v∂t = α(cid:126)v − β | (cid:126)v | (cid:126)v − σ (cid:126) ∇ ( ρ − ρ ) + (cid:126)f , ∂ρ∂t + (cid:126) ∇ · ( (cid:126)vρ ) = 0 . In a limit where deviations from the mean density are small, so ρ = ρ + η(cid:37) for some η (cid:28)
1, we can further reduce to ∂(cid:126)v∂t = α(cid:126)v − β | (cid:126)v | (cid:126)v − ησ (cid:126) ∇ (cid:37) + (cid:126)f , η ∂(cid:37)∂t + ( ρ + η(cid:37) ) (cid:126) ∇ · (cid:126)v + η(cid:126)v · (cid:126) ∇ (cid:37) = 0 . Then on short time scales τ = t/η , we have ∂(cid:126)v∂τ = ηα(cid:126)v − ηβ | (cid:126)v | (cid:126)v − η σ (cid:126) ∇ (cid:37) + η (cid:126)f , ∂(cid:37)∂τ ≈ − ρ (cid:126) ∇ · (cid:126)v, where we neglect terms that must be of order η : if the coefficients α , β , and σ are sufficiently large, their terms willremain relevant. The scaling of τ ensures that t is small when τ is order one or smaller. Discretizing the velocity anddensity fields as well as the noise (cid:126)f and replacing the continuous gradient with either ∇ (cid:62) ev or −∇ ve as appropriateyields Eqs. (4). Compressibility
Compressibility as included in our model is intended to describe changes in density or volume of the active com-ponent, not the underlying fluid. For example, variations in (cid:37) may be interpreted as variations in the density ofswimmers in a bacterial system or variations in the tube volume in
Physarum polycephalum . Such systems may beeffectively compressible even though the solvent fluid (e.g. water) is incompressible.In some cases, compressibility is the primary object of interest. For example, a recent preprint [15] discussessound in active fluids in a network using a continuous wave equation derived from the Toner-Tu model. On top ofa background flow taking the form of a lattice of counter-rotating cycles, they find modes confined to the edges of aLieb lattice, which we can reproduce in our discretized setting (Fig. S1 and Movie 3). In both their setting and ours,these edge modes decay over time without propagating into the bulk (cf. discussion in App. I.B of Ref. [15]).We can recover an incompressible limit of our model by first extending it to include damping on the vertices: d(cid:37) v dt = (cid:88) e ∇ ve φ e − η(cid:37) v , (11a) dφ e dt = − γ (cid:88) v ∇ (cid:62) ev (cid:37) v + (cid:15)g ( µ, φ e ) φ e + √ Dξ e ( t ) . (11b) Fig. S1. Our active network model exhibits behavior similar to the topological edge modes of [15]. (a) A discretized version ofthe Lieb lattice considered in [15]. Edges shared by adjacent 8-cycles have weight γ e = 2 to account for the additional width ofthe corresponding channels. The most stable flow on this network consists of a lattice of counter-rotating cycles, in which boththe active friction term g ( µ, φ e / √ γ e ) and the pressure variations (cid:37) v are everywhere zero. (b) This lattice has modes confinedto the edges of the domain, allowing sound waves to propagate and decay without scattering into the bulk (cf. discussion inApp. I.B of Ref. [15]); one such mode is pictured. Simulations started in this mode as a perturbation to the most stable flowpattern do not cause density changes in the center (Movie 3). The network model allows study of such phenomena withoutresorting to full scale simulation of the flow patterns. This paper examines the limit η → η → ∞ , where Eq. (11a) can only be balanced if (cid:37) v → (cid:37) v = 1 η (cid:88) e ∇ ve φ e . Substituting this into Eq. (11b) gives dφ e dt = − γη (cid:88) v ∇ (cid:62) ev ∇ va φ a + (cid:15)g ( µ, φ e ) φ e + √ Dξ e ( t ) . With g ( µ, φ e ) = φ e (1 − φ e ), this is equivalent to the model discussed in [32]. If γ → ∞ so that γ/η is constant, smalldeviations from incompressibility are allowed; if γ/η → ∞ , incompressibility is fully enforced. However, compressibilityis a necessary ingredient for sound waves [15] and density oscillations [20]. Rayleigh friction approximation
While choosing the friction function to be [29] g ( µ, φ e ) = µ − φ e φ e has convenient theoretical properties, namely that it gives a passive constant friction coefficient (cid:15) for µ = − φ → ∞ , it is analytically difficult. To simplify the analysis, we approximate this g ( µ, φ e ) with a symmetricquadratic [27] ˆ g ( µ, φ e ) = a − bφ e , (12)where a = µ and b = 1 are chosen so that ˆ g ( µ,
0) = g ( µ,
0) and ˆ g ( µ, φ e ) has the same zeros as g ( µ, φ e ). This ensuresthat the two functions approximately match when they are both negative, that is, when activity is putting energyinto the flow. The large difference between g ( µ, φ e ) and ˆ g ( µ, φ e ) when the flux is large is less important, as the flowwill be damped down in either case. The larger damping in ˆ g ( µ, φ e ) does result in slightly lower steady amplitudes,both analytically and in simulations.0 Perturbation expansion If (cid:15) is small, there will be two widely separated timescales: the fast oscillation timescale t and the slow frictiontimescale τ = (cid:15)t . After writing (cid:37) v and φ e in the mode basis, we can further expand in (cid:15) as r n ( t ) = ∞ (cid:88) k =0 (cid:15) k r kn ( t, τ ) , (13a) f n ( t ) = ∞ (cid:88) k =0 (cid:15) k f kn ( t, τ ) , (13b)where we explicitly separate the dependence on the two timescales. Then¨ r kn ( t, τ ) = ∂ t r kn + 2 (cid:15)∂ t ∂ τ r kn + (cid:15) ∂ τ r kn , ¨ f kn ( t, τ ) = ∂ t f kn + 2 (cid:15)∂ t ∂ τ f kn + (cid:15) ∂ τ f kn . At zeroth order in (cid:15) , with D = 0, Eq. (5) becomes V (cid:88) n =1 ∂ t r n (cid:37) vn = − V (cid:88) n =1 λ n r n (cid:37) vn . The modes (cid:37) vn are orthonormal, so the terms decouple into separate harmonic oscillators; f kn can be found from r kn using Eq. (4a). The leading order solution is then r n ( t ) = A n ( τ ) cos( λ n t − δ n ( τ )) ,f n ( t ) = − A n ( τ ) sin( λ n t − δ n ( τ )) . At first order in (cid:15) , with g ( µ, φ e ) = ( µ − φ e ), V (cid:88) n =1 ( ∂ t r n + 2 ∂ t ∂ τ r n ) (cid:37) vn = − V (cid:88) n =1 λ n r n (cid:37) vn + (cid:88) e ∇ ve µ − (cid:32) E (cid:88) n =1 f n φ en (cid:33) E (cid:88) l =1 f l φ el . Multiplying by (cid:37) vm and summing over v , we find ∂ t r m + 2 ∂ t ∂ τ r m = − λ m r m + λ m µf m − (cid:88) e φ em (cid:32) E (cid:88) n =1 f n φ en (cid:33) . (14) Leading order amplitude dynamics
In order for the expansion in Eqs. (13a) and (13b) to make sense, the magnitudes of the summands r kn and f kn must remain bounded. From Eq. (14), r m is a harmonic oscillator with natural frequency λ m driven by the zerothorder oscillations. It will have bounded oscillations only if the resonant terms in Eq. (14), those that drive r m atits natural frequency, are zero. Finding the resonant terms and setting them to zero will fix the leading order modeamplitudes A n ( τ ).Expanding the cube in Eq. (14) gives ∂ t r m + 2 ∂ t ∂ τ r m = − λ m r m + λ m (cid:20) µf m − (cid:88) e φ em E (cid:88) k,(cid:96),n =1 f k φ ek f n φ e(cid:96) f n φ en (cid:21) = − λ m r m + λ m (cid:20) µf m + E (cid:88) k,(cid:96),n =1 (cid:32)(cid:88) e φ em φ ek φ e(cid:96) φ en (cid:33) × A k A (cid:96) A n sin( λ k t − δ k ) sin( λ (cid:96) t − δ (cid:96) ) sin( λ n t − δ n ) (cid:21) . (15)1Now, the product of sines can be expanded intosin( λ k t − δ k ) sin( λ (cid:96) t − δ (cid:96) ) sin( λ n t − δ n ) = 14 (cid:104) sin( δ k − δ (cid:96) − δ n − λ k t + λ n t + λ (cid:96) t ) − sin( δ k − δ (cid:96) + δ n − λ k t − λ n t + λ (cid:96) t ) − sin( δ k + δ (cid:96) − δ n − λ k t + λ n t − λ (cid:96) t )+ sin( δ k + δ (cid:96) + δ n − λ k t − λ n t − λ (cid:96) t ) (cid:105) . We seek only resonant terms, which only occur when ± λ k , ± λ (cid:96) , and ± λ n sum to λ m . This happens most oftenin one of two ways. First, we might have k = (cid:96) and n = m or similar. Alternatively, we might have degeneratemodes, λ k = λ (cid:96) and λ n = λ m . However, we ignore the latter possibility because degeneracies add significant analyticcomplications, including nontrivial dynamics of their relative phases. We also ignore the rare possibility of resonantterms arising from interactions of modes with three or four distinct singular values. The results we get with theseassumptions closely match simulated time series (Fig. 3e-g), suggesting that the existence of degeneracies has littleimpact on the dynamics of nondegenerate modes.The remaining resonant terms in Eq. (15) must cancel so that r m is not an oscillator of frequency λ m driven atfrequency λ m . Thus,2 ∂ t ∂ τ r m = λ m µf m + 14 (cid:32)(cid:88) e φ em (cid:33) A m (3 sin( λ m t − δ m )) + 3 E (cid:88) k =1 ,k (cid:54) = m (cid:32)(cid:88) e φ em φ ek (cid:33) A k A m
14 (2 sin( λ m t − δ m )) . Substituting in r m and f m , − A (cid:48) m λ m sin( λ m t − δ m ) + 2 λ m cos( λ m t − δ m ) δ (cid:48) m = λ m (cid:34) − µA m sin( λ m t − δ m ) + 14 (cid:32)(cid:88) e φ em (cid:33) A m (3 sin( λ m t − δ m ))+ 3 E (cid:88) k =1 ,k (cid:54) = m (cid:32)(cid:88) e φ em φ ek (cid:33) A k A m
14 (2 sin( λ m t − δ m )) , where primes denote differentiation with respect to τ . For this to hold for all t we need the coefficients of the sineand cosine terms to separately cancel. From the cosine term, δ (cid:48) m = 0; from the sine term, A (cid:48) m = 12 A m µ − (cid:32)(cid:88) e φ em (cid:33) A m − E (cid:88) k =1 ,k (cid:54) = m (cid:32)(cid:88) e φ em φ ek (cid:33) A k ≡ A m (cid:32) µ − E (cid:88) k =1 P mk A k (cid:33) , where the matrix P has entries P mk = (1 − δ mk ) (cid:80) e φ em φ ek . Rewriting in terms of the squared amplitudes, ddτ ( A m ) = 2 A m A (cid:48) m = A m (cid:32) µ − E (cid:88) k =1 P mk A k (cid:33) . (16)As a matrix equation, with x m = A m , this reads x (cid:48) = x (cid:12) ( µ − P x ) , (17)where denotes the vector of ones and (cid:12) is the component-wise product.To find stationary points, we set x (cid:12) ( µ − P x ) = 0. The obvious way to solve Eq. (17) for all stationary points isto exhaustively search over combinations of active modes: on picking certain elements of x to be zero, the remainingnonzero entries ˆ x are found by solving ˆ P ˆ x = µ , where ˆ P is P restricted to those modes chosen to be nonzero.Stability of a fixed point x then follows by standard perturbation analysis: inserting a small perturbation x + δ x ( τ )into Eq. (17) gives δ x (cid:48) = δ x − x (cid:12) ( P δ x ) − ( P x ) (cid:12) δ x + O ( δ x ) ≡ M δ x + O ( δ x ) , where I denotes the identity matrix, and the eigenvalues of M then determine stability in the usual fashion.2 Fig. S2. Steady state amplitudes A i as a function of activity µ for the tree pictured undergo a Hopf bifurcation as µ crosses 0.Dots are long-time root-mean-square amplitudes from simulations started in each mode; lines are numerical solutions of Eq. (18).Mode A is too unstable to reliably observe in simulations, so it is omitted. For µ <
0, all amplitudes go to zero in simulations;the dot included in that region is at µ = − (cid:15) (cid:54) = 0. Parameters were (cid:15) = 0 . D = 0. Accuracy of Rayleigh friction approximation
To verify that the Rayleigh friction approximation does not significantly impact the results, we check the amplitudeand stability of single modes for the full model with g ( µ, φ e ) = ( µ − φ e ) / (1 + φ e ) on all edges. Here setting the firstorder secular terms to zero in a perturbation expansion with A n = A p δ np leads to A p = ( µ + 1) (cid:88) e (cid:32) − (cid:115)
11 + A p φ ep (cid:33) . (18)Numerically solving Eq. (18) for µ = 1 yields solutions within a few percent of the Rayleigh approximation solution1 / (cid:112) P pp which additionally match numerical simulations of the full model even for (cid:15) as large as 0 . µ .If we assume µ (cid:28) A p (cid:28)
1) and expand the square root to order A p , we find A p + O ( A p ) = µ/P pp , exactlymatching the Rayleigh friction result. The scaling A p ∼ √ µ is typical of a supercritical Hopf bifurcation. Attractor characteristics on tree networks
The mode interactions of Eq. (17) can lead to complex oscillation patterns dependent on global, not local, topology,as shown for a 127-vertex complete binary tree in Movie 4 and Fig. S3. After initializing with zero pressure variationand flux, the system settles into quasi-steady states with dramatically different dynamics in separate regions of thetree (Fig. 3a,b). Flux in edges near the leaves of the tree tends to oscillate rapidly, driving large pressure fluctuationsin nearby vertices, whereas flux oscillations near the root are comparatively slow with nearly constant pressure in thevertices (Fig. 3b,d). Since, apart from the root and leaves, each vertex has the same local topology, the different timescales emerge from the interaction of the local active friction with the global structure of the tree.A comprehensive and precise characterization of the relative lifetimes of different attractors in large active flownetworks remains out of reach with current numerical methods, in part because the range of noise levels low enoughto observe state selection and high enough to observe transitions is quite small. Such a fine-tuning between thermaland active transport processes is a characteristic feature of many, if not all, biological systems that function optimallyin a narrow temperature range: bacterial flagellar motors are designed to barely beat Brownian diffusion at roomtemperature, ATP-driven intracellular transport is tuned such that it improves moderately over thermal diffusion,and so-on. Another well-known example in this context is stochastic resonance in driven multistable systems [45].However, as all these systems typically exhibit exponential Arrhenius-type waiting times, it is practically impossible3to completely explore their attractor statistics in the moderate-to-weak noise regime, except for the simplest two-statesystems [46].Nevertheless, long simulation runs as shown in Fig. S4 offer some insight into the qualitative behavior of attractorsin active flow networks. Specifically, our simulations suggest that, while there is considerable variation in the relativeoccupancy of different attractors, stable states can be approximately divided in two classes: (1) states with one highenergy mode at high amplitude and a few low energy modes at low amplitude and (2) states with multiple low-energymodes active at moderate amplitude, some of them degenerate. States of type (2) tend to quickly transition to otherstates of type (2) (Fig. S4); states of type (1) have a wide range of lifetimes but no obvious transition patterns.
Fig. S3. Activity causes depth-dependent separation of time scales on a large tree. (a) Most pressure variation occurs nearthe leaves on large binary trees (Movie 4). (b) The tree in (a) develops an activity-driven steady state with slow oscillations inthe center and fast oscillations near the edges, as illustrated by the flux φ e on the three edges labelled in (a). (c) Unnormalizedcorrelations between the Fourier transforms of the flux through the edges of the tree in (a), with phases ignored. Colors indicatethe tree level of the tail vertex of the edge. There are strong correlations within each level and between neighboring levels,but low correlations for edges in widely-separated levels. (d) Frequency spectra of each tree level, computed by taking Fouriertransforms of the edge fluxes as in (c) and averaging the magnitudes across all edges at each level. A distinct primary oscillationfrequency for each level can be seen, which increases with distance from the tree center. Simulation parameters in all panelsare (cid:15) = 0 . µ = 1, and D = 10 − . (e-h) While adding edges in the center leads to steady flow on cycles there, frequency stillincreases with distance from the center in the outer, tree-like sections. Fig. S4. Lower energy modes transition more often for the graph in Fig. 3e of the Main Text. Modes are ordered by frequencyfrom high (top) to low (bottom). Simulation parameters are (cid:15) = 0 . µ = 1, D = 5 × − , identical to those in Fig. 3. Notethat rows 7 and 8, the two modes that switch on and off most, are degenerate. Networks with cycles
We focus on tree networks in this paper as they allow substantial analytical progress. However, Eqs. (4) can beapplied without modification to networks with cycles. Cycles correspond to right singular vectors φ n of ∇ (cid:62) withsingular value zero. As these are always degenerate, we expect the conclusions of Section to be most accurate whenthere are few or no cycles. Alternatively, on a weighted graph where the edges of high conductance form a tree, theattractor characteristics will be similar to the attractors on that tree (Fig. 1; all modes pictured in Fig. S6).Qualitatively, we find the same stochastic switching between states with subsets of modes active in simulations ofEqs. (4) on cyclic graphs even with equal weights, with the additional feature that cyclic modes are particularly stableand take longer to transition on average (Fig. S5). For further discussion of similar dynamics on cycles, see [32]. Fig. S5. States on graphs with cycles, like the one shown, tend to be more stable. Modes are ordered by frequency from high(top) to low (bottom). Note that the eight modes at the bottom, which are the only ones active in the lower half of the trace,are all cycles. Simulation parameters are (cid:15) = 0 . µ = 1, D = 5 × − . Fig. S6. Including all of the modes from the simulation in Fig. 1 of the Main Text shows clear single mode selection on thisweighted network. Edges a distance d from the central red path were given weight e − d . Modes are ordered by frequency fromhigh (top) to low (bottom); the last thirty modes, marked in red, are cycles. The modes pictured in Fig. 1 are marked in black. Higher order oscillations
Before, by setting resonant terms to zero, we found the slow dynamics of A n . Now we look at the non-resonantterms driving r m to find higher order effects. If we let S i n ...i nkk = (cid:88) e k (cid:89) j =1 φ n j ei j , Fig. S7. Slow global oscillations emerge from the fast active dynamics. (a) First order considerations fix a constant mean flowenergy; higher order effects cause significant slow oscillations about that mean. Simulation parameters were µ = 1, (cid:15) = 0 . D = 0; the tree used is inset. (b) The mode amplitudes A and A , like the energy, oscillate much more slowly than theharmonic oscillations of f and f . All other mode amplitudes (unlabelled traces) are close to zero. (c) Frequency spectra ofthe two active modes and the energy H for the simulation in (a) and (b). The energy oscillates due to higher-order interactionsbetween modes at frequencies that are linear combinations of active mode frequencies, not the harmonic frequencies alone(dashed lines). assume the resonant terms are zero, and assume A m = A p δ mp + A q δ mq , the remainder of Eq. (14) is ∂ t r m + λ m r m = 14 λ m (cid:110) S mp A p sin(3 λ p t ) + 3 S mq p A p A q (cid:2) sin((2 λ q − λ p ) t ) − sin((2 λ q + λ p ) t ) (cid:3) + 3 S mqp A p A q (cid:2) sin((2 λ p − λ q ) t ) − sin((2 λ p + λ q ) t ) (cid:3) + S mq A q sin(3 λ q t ) (cid:111) . (19)Setting m = p and only looking at the terms closest to resonance, we obtain ∂ t r p + λ p r p ≈ λ p (cid:8) S q p A p A q sin((2 λ q − λ p ) t ) + 3 S p q A p A q sin((2 λ p − λ q ) t ) (cid:9) . Thus r p ≈ c cos((2 λ q − λ p ) t − δ ) + c cos((2 λ p − λ q ) t − δ ) ,f p ≈ − c sin((2 λ q − λ p ) t − δ ) − c sin((2 λ p − λ q ) t − δ ) , where c = 34((2 λ q − λ p ) − λ p ) λ p S q p A p A q ,c = 34((2 λ p − λ q ) − λ p ) λ p S qp A p A q . The energy in this mode to first order in (cid:15) is H p = λ p (cid:0) ( r p + (cid:15)r p ) + ( f p + (cid:15)f p ) (cid:1) + O ( (cid:15) )= λ p (cid:110) A p + 2 (cid:15)A p (cid:2) c cos((2 λ q − λ p ) t ) + c cos(( λ p − λ q ) t ) (cid:3)(cid:111) + O ( (cid:15) ) , exhibiting an order (cid:15) time dependence. The coefficients c and c are small unless λ p ≈ λ q . If we kept the frequency3 λ p , 3 λ q , 2 λ p + λ q , and 2 λ q + λ p terms from Eq. (19), we would find energy oscillations with frequencies 2 λ p , 2 λ q ,3 λ q − λ p , and λ p + λ q (Fig. S7); those oscillations have smaller amplitudes as the driving is farther from resonance. Noise and thermalization
In Eqs. (4a) and (4b) we add Gaussian white noise only to the flux as a physically intuitive source of randomfluctuations that preserve mass conservation. However, even with purely passive friction, this does not lead toequipartition of energy as seen in thermal systems.7Written as stochastic differential equations with g ( µ, φ ) = −
1, Eqs. (4a) and (4b) become d(cid:37) v = (cid:88) e ∇ ve φ e dt, (20a) dφ e = − (cid:88) v ∇ (cid:62) ev (cid:37) v dt − (cid:15)φ e dt + √ Dd ˜ B e ( t ) , (20b)where each ˜ B e ( t ) is standard Brownian motion. The components of the E -dimensional Brownian motion B ( t ) =( ˜ B , . . . , ˜ B E ( t )) in any orthonormal basis are also standard Brownian motions, so we can rewrite the system in themode basis as dr n = λ n f n dt, (21a) df n = − λ n r n dt − (cid:15)f n dt + √ DdB n ( t ) . (21b)The associated Fokker-Planck equation for the probability distribution p ( r , f , t ) is ∂ t p = (cid:88) n (cid:20) − ∂∂r n ( λ n f n p ) + ∂∂f n ( λ n r n p ) + ∂∂f n ( (cid:15)f n p ) + D ∂ p∂f n (cid:21) with p → r n , f n → ∞ and p integrating to 1. Now, without friction or noise, the dynamics are governed by theHamiltonian H = 12 (cid:37) v ∇ ve ∇ eu (cid:37) u + 12 φ e ∇ ev ∇ va φ a = 12 (cid:88) n λ n ( r n + f n ) ≡ (cid:88) n H n . If p is a function of the H n alone, the Fokker-Planck equation in steady state reduces to0 = (cid:88) n ∂∂f n ( (cid:15)f n p ) + D ∂ p∂f n , which has solution p ( H , . . . , H M ) ∝ M (cid:89) n =1 e − HnkTn , where kT n = λ n D/(cid:15) .Loosely, adding noise this way couples each mode to a heat bath with a distinct temperature. The result isequipartition of amplitude, not energy: the long-time average (cid:104) A n (cid:105) is independent of n . Adding weak couplingbetween modes by making µ > − (cid:112) D n dB n ( t ) ≡ √ Dλ n dB n ( t ) . This is only possible for λ n (cid:54) = 0, which precludes cyclic modes. Equation (4b) becomes dφ e = − (cid:88) v ∇ (cid:62) ev (cid:37) v dt + (cid:15)g ( µ, φ e ) φ e dt + (cid:88) n λ n φ en √ DdB n ( t ) . The previous analysis goes through identically, leading to kT n = λ n D n /(cid:15) = D/(cid:15) . Differential growth rates
While the E/ E . There are, however, several straightforward generalizations of ourmodel that may lead to more strict mode selection. We discuss two possibilities in this and the subsequent section:variations in activity across the network and variations in weights of vertices or edges.For simplicity, we introduced Eqs. (4) with a uniform activity level µ across the entire network. This leads to equaldriving on all modes: if Eq. (16) is initialized near zero, it can be linearized to ddτ ( A m ) = µA m , µ in Eqs. (4) with edge-dependent parameters µ e . With thequadratic driving of Eq. (12), Eq. (4b) becomes dφ e dt = − (cid:88) v ∇ (cid:62) ev (cid:37) v + (cid:15) (cid:0) µ e − φ e (cid:1) φ e + √ Dξ e ( t ) . Following through the previous calculations with this change, Eq. (14) becomes ∂ t r m + 2 ∂ t ∂ τ r m = − λ m r m + λ m (cid:88) e,l φ em µ e f l φ el − (cid:88) e φ em (cid:32) E (cid:88) n =1 f n φ en (cid:33) . The first term inside the square brackets no longer simplifies, since the φ en are not orthonormal with the weighting µ e . However, if we again ignore degeneracies, the only resonant term is (cid:80) e φ em µ e f m from l = m . In this case,defining ν m = (cid:80) e φ em µ e , Eq. (16) then reads ddτ ( A m ) = A m (cid:32) ν m − E (cid:88) k =1 P mk A k (cid:33) , (22)where modes have distinct growth rates independent of their interactions. Alternatively, one could specify ν m arbi-trarily in Eq. (22), though this would require more complex changes in Eq. (4b) coupling activity across edges. Band gaps
In addition to distinct activity levels µ e across edges, we can also introduce edge weights γ e or vertex weights m v that vary across the network. Changing the conductances γ e and volumes m v changes our system in two ways:first, by changing the modes to the singular vectors of γ ∗ ve ; and second, by changing the coupling matrix to ˜ P mk = (1 − δ mk ) (cid:80) e γ − e φ em φ ek , which depends explicitly on the edge weights.Such changes are known to cause qualitative changes in the physics of classical spring-mass networks, including theintroduction of band gaps. In an infinite one-dimensional line of beads of equal mass m connected by springs withequal spring constant f , for example, the dispersion relation between frequency ω and wavenumber q is ω ( q ) = 2 (cid:114) fm (cid:12)(cid:12)(cid:12) sin (cid:16) qa (cid:17)(cid:12)(cid:12)(cid:12) , Fig. S8. The emergence of an activity-driven spectral band gap is exhibited by a simulation on a 14-vertex path with (a)all weights equal to 1 and (b) alternating vertex weights 1 and 5. Modes are ordered by frequency from high (top) to low(bottom). Note that in (b) the central n = 7 mode is always active and the low energy states on the right half of the plotare significantly more suppressed than they ever are in (a). The qualitative difference is due to the presence of vertices withunequal weights, not the overall scale of the vertex weights; changing vertex weights uniformly is equivalent to rescaling otherparameters. Parameters were µ = 1 . D = 5 × − , and (cid:15) = 0 .
5. Both simulations used the same random seed. a is the size of the unit cell, in this case equal to distance between adjacent beads [48]. If instead of equalmasses the beads alternate between a smaller mass m and larger mass m , the dispersion relation splits into twobranches, ω ( q ) ± = f (cid:18) m + 1 m (cid:19) ± f (cid:115)(cid:18) m + 1 m (cid:19) − m m sin (cid:16) qa (cid:17) . Here a unit cell has two beads, so the distance between beads is a/
2. At q = π/a , there is a gap between ω + = (cid:112) f /m and ω − = (cid:112) f /m2