Dynamical phase separation on rhythmogenic neuronal networks
DDynamical phase separation on rhythmogenic neuronal networks
Mihai Bibireata, Valentin M. Slepukhin, and Alex J. Levine
1, 2, 3 Department of Physics and Astronomy, UCLA, Los Angeles California, 90095-1596, USA Department of Chemistry and Biochemistry, UCLA, Los Angeles California, 90095-1596, USA Department of Computational Medicine, UCLA, Los Angeles California, 90095-1596, USA (Dated: January 10, 2020)We explore the dynamics of the preB¨otzinger complex, the mammalian central pattern generatorwith N ∼ neurons, which produces a collective metronomic signal that times the inspiration.Our analysis is based on a simple firing-rate model of excitatory neurons with dendritic adaptation(the Feldman Del Negro model [Nat. Rev. Neurosci. 7, 232 (2006), Phys. Rev. E 2010 :051911])interacting on a fixed, directed Erd˝os-R´enyi network. In the all-to-all coupled variant of the model,there is spontaneous symmetry breaking in which some fraction of the neurons become stuck ina high firing-rate state, while others become quiescent. This separation into firing and non-firingclusters persists into more sparsely connected networks, and is partially determined by k -cores in thedirected graphs. The model has a number of features of the dynamical phase diagram that violatethe predictions of mean-field analysis. In particular, we observe in the simulated networks thatstable oscillations do not persist in the large-N limit, in contradiction to the predictions of mean-field theory. Moreover, we observe that the oscillations in these sparse networks are remarkablyrobust in response to killing neurons, surviving until only ≈
20% of the network remains. Thisrobustness is consistent with experiment.
I. INTRODUCTION
Networks of connected neurons, or neuronal microcir-cuits, play a variety of roles [1]. Their collective dynami-cal properties depend upon both their network structure, i.e. , the pattern of which neurons are synaptically cou-pled to each other, and the neurons’ individual proper-ties, i.e. , the relationship between input and output atthe level of a single neuron. The preB¨otzinger complex(preB¨otC) is a particular neural microcircuit consisting of ∼ neurons that acts as a central pattern generator,establishing the inspiratory rhythm in mammals [2–4].It is particularly interesting for theoretical study sinceit may be modeled as a network of essentially identicalneurons that produces well-characterized collective dy-namics – a metronomic oscillation between states of highfiring rate and low firing rate, which sets the inspiratoryrhythm.Feldman and Del Negro [4] suggested a simple model ofthe PreB¨otC neurons that includes a combination of ex-citatory interactions and dendritic adaption, which wasstudied theoretically by Schwab et al. [5] using a sim-ple leaky integrate-and-fire model for the neurons [1, 6].The principal idea enabling an oscillatory phase of thenetwork is the introduction of a slow internal variable(identified by experiment to be dendritic calcium con-centration [6]) in each neuron controlling its dendriticsensitivity to excitatory postsynaptic potentials (EPSP).The calcium concentration increases with each EPSP. Inthe proposed model, when the calcium concentration isbelow a threshold, the neurons are sensitive to externalvoltage signals (EPSPs) and, though mutual excitatoryinteractions, they collectively increase their firing rate.This collective high-firing-rate period represents the in-spiratory signal. When many neurons have a high firingrate, however, the total input voltage to a typical neu- ron becomes sufficiently high that its dendritic calciumconcentration rises above the threshold rendering thatneuron insensitive to further input. As a result the neu-rons’ firing rate rapidly decreases and remains small untilthe next period of mutual excitation that occurs once thedendritic calcium concentration of a typical neuron hasfallen below the threshold restoring those neurons’ sensi-tivity to input.In addition to the stably oscillating phase of the net-work, the model was shown to admit two other phasescharacterized by either steady-state low-firing rate orhigh firing-rate throughout the network. We refer tothese as the quiescent and high-activity phases respec-tively. A dynamical phase diagram of this system wasobtained as a function of network size and basal (low cal-cium) neuronal excitability both in a mean-field analysisand numerically on a set of Erd˝os -R´enyi (ER) graphs.Intriguingly, the numerically obtained phase boundarybetween the stably oscillating and high activity phasesdemonstrates significant deviations from the mean-fieldtheory predictions, with discontinuous “jumps” whoseposition on the phase diagram corresponds to numbersof neurons at which the highest k k -core [7] of the net-work vanishes [5].In this manuscript, we further explore the Feldman-DelNegro model, showing that the dynamical phase space ismuch richer than previously thought. In particular, wefind that region of phase space consistent with stable os-cillations is bounded in both network size and neuronalbasal excitability. This is inconsistent with the mean-field predictions. Moreover, most of the phase bound-aries are jagged (the interfaces are rough), which is alsoinconsistent with the mean-field predictions for this dy-namical system. The main reason for the failure of themean-field approximation is due to the spontaneous sep-aration of the neurons into groups with disparate firing a r X i v : . [ q - b i o . N C ] J a n rates. In fact, even in an all-to-all coupled network, thepermutation symmetry of the neurons is spontaneouslybroken into these high and low firing-rate groups. Wecall this spontaneous activity phase separation. We an-alyze this separation process on random networks bothnumerically and analytically, showing that the connectiv-ity disorder of the random networks guides the separationprocess.We organize the rest of the manuscript as follows. Insection II, we demonstrate spontaneous activity separa-tion into high and low firing-rate groups on a small net-work. We show how that separation leads to the observedroughness of the phase boundaries. For the special caseof an all-to-all coupled network, one can analytically de-rive activity phase separation. We do so and comparethese results to numerical simulations on all-to-all cou-pled networks in section III. From that analysis we learnthat the steepness of the neuronal firing rate function(as a function of somatic potential) controls this sponta-neous symmetry breaking on the network. In section IV,we move to the case of more sparsely connected net-works, chosen from the ensemble of ER networks, wherewe prove that the activity separated solution, if it exists,is stable for sufficiently sharp neuronal firing-rate func-tions. The cases where such activity separated solutionsdo not exist is reminiscent of converse symmetry break-ing [8], where symmetric solutions can be paradoxicallystabilized by system asymmetry. Finally, in section Vwe consider the role of k -cores in determining which neu-rons fall into the high-activity state in sparsely connectednetworks. We prove that, when setting the low somaticvoltage firing rate to zero, activity phase separation isexactly controlled by the k -cores. We suggest that k -cores remain relevant in controlling the phase boundarybetween the quiescent and high activity phases of thedisordered system, but these topological features cannotalone account the roughness of the high-activity/stableoscillation phase boundary. We conclude with a com-parison of these results to current experiments and makesuggestions for future work.The reader interested primarily in our experimentalpredictions for in vitro PreB¨otC system may turn to sec-tion VI. Readers interested in the simulations will findthe reference to the software and appropriate parametersin Appendix E.
II. THE FELDMAN DEL NEGRO MODEL
Following Ref. [5], we describe the two-compartmentneuron model of preB¨otC neurons. The i th neuron ischaracterized by two dynamical variables, its somaticpotential V i and its dendritic calcium concentration C i . Their dynamics are controlled by the equations dV i dt = 1 τ V ( V eq − V i ) + ∆ V ( C i ) (cid:88) j M ij r ( V j ) (1) dC i dt = 1 τ C ( C eq − C i ) + ∆ C (cid:88) j M ij r ( V j ) , (2)where ∆ V ( C ) and r ( V ) are defined by∆ V ( C ) = ∆ V max σ (cid:18) C ∗ − Cg C (cid:19) (3)and r ( V ) = ( r m − r b ) σ (cid:18) V − V ∗ g V (cid:19) + r b . (4)In Eqs. 3, 4 we have introduced the standard sigmoid(Fermi) function σ ( x ) = 11 + e − x . (5)Here and throughout the manuscript we work in dimen-sionless calcium concentration units obtained by setting C eq = 0 and C ∗ = 5.Eq. 1 is typical of a leaky integrate and fire model foran excitatory neuron. The principal addition in the two-compartment model is dendritic adaptation which is builtinto ∆ V ( C i ) defined in Eqs. 3 and 5. An incoming EPSPproduces both an increase in dendritic calcium concen-tration ∆ C and somatic potential ∆ V ( C ). But above athreshold concentration C ∗ , ∆ V ( C ) becomes small, ren-dering the neuron insensitive to EPSPs. In the absenceof incoming EPSPs, the dendritic calcium concentrationreturns to C eq on a time scale of τ C at which point theneuron is once again sensitive to EPSPs.The parameter space of the neuron model is controlledby a small set of physiological constants. There arethe steady-state dendritic calcium concentration and so-matic potential C eq and V eq respectively. The voltage-dependent firing rate is determined by the basal andmaximal firing rates r b and r m as well as g V , which con-trols the steepness of the transitions around the thresh-old voltage V ∗ . Dendritic adaptation is parametrizedby the maximum voltage increment associated with anEPSP ∆ V max , a calcium concentration threshold C ∗ anda steepness parameter g C , analogous to g V discussedabove. In addition to the two time scales τ V < τ C forthe relaxation of somatic potential and dendritic calcium,there is a fixed calcium concentration increment ∆ C as-sociated with the response to an EPSP. Table I providesthe currently available values of the model parameters.The model also depends on the size and connectivity ofthe underlying network of synaptic connections betweenthe neurons. The network’s structure can be encoded byan adjacency matrix M whose matrix elements M ij =1 if neuron i synapses on neuron j , and equal to zerootherwise. In this manuscript, we consider only networks TABLE I. Model parameters known from experiment
Parameter Approximate value References V eq -65 mV [1] V ∗ -50 mV [1] τ V
20 ms [9],[10] r m
40 Hz [9],[10] r b V max p N [9], [11] built from uncorrelated stochastic connections – Erd˝osR´enyi directed graphs [12]. An ensemble of such networksis determined by a single probability p that any non-diagonal matrix element is equal to one. We exclude thepossibility of a neuron synapsing on itself. The all-to-all network is simply the case of such an ER graph with p = 1. A. Dynamical phase diagram
For a given set of parameters and given network of N neurons, the dynamical system evolves deterministi-cally from a set of 2 N initial conditions leading to eithera fixed point, limit cycle, or chaotic dynamics at longtimes. We find multiple fixed points, which can be fur-ther distinguished as quiescent (Q) where the somaticpotential averaged over the network of neurons (cid:104) V (cid:105) liesbelow the transition to the high-firing state V ∗ , or highactivity (HA), where (cid:104) V (cid:105) > V ∗ . Similarly, we can distin-guish three classes of stable limit cycle oscillations: belowthreshold oscillations (BTO) where the oscillatory aver-age voltage remains below the threshold for high firingrate, above threshold oscillations (ATO), where the os-cillatory average voltage remains above the threshold forhigh firing rate, and true metronomic activity (TMA),where the stable limit cycle oscillations carry the systembetween high and low firing rates, producing the physio-logically observed inspiratory rhythm.To examine the dynamical phase behavior of the sys-tem, we vary the basal excitability of the neurons andthe size of their network – ∆ V max and N – while fix-ing the rest of the parameters. Typical results for net-works are shown in Fig. 1 for three different choices ofthe fixed variables. We observe in Fig. 1A the numeri-cally determined phase diagram for all-to-all coupled net-works, which agrees with the mean field solution of themodel shown in Appendix C. In general, we find thatthe numerically determined phase diagram agrees withthe mean field approximation (for arbitrary initial con-ditions) as long as the transition in dendritic sensitivityis sufficiently smooth, i.e. g C (cid:38) i.e. , decrease g C (cid:46) V in the mean-field prediction. Inthis regime, we do not see chaotic dynamics unlike in thecases where g C is small. Changing other parameters ofthe model changes the shape of the phase boundaries, butdoes not introduce new dynamical phases. Both routes tothe breakdown of mean field theory (small g C and moresparsely connected networks) are related to an inherentinstability toward activity phase separation. We discussthis in more detail below.Before discussing the phase separation, we note thatthe roughness of the phase diagram in the non-mean-fieldregime, as shown in Fig. 1C, implies that the physiolog-ically desirable stably oscillating phase (TMA) admits atype of reentrant behavior in which one can remove neu-rons (decrease N ) from a network in the high activitystate to render it in the stably oscillating TMA phase.In Fig. 2 we see examples of such possible transitions at∆ V = 18mV where by decreasing the number of neuronsfrom N = 90 to N = 65, one encounters two transitionsHA to TMA, before remaining in the HA phase below N = 72. This suggests a specific experimental test of thefundamental model that can be made by looking for thesereentrant dynamical transitions upon killing neurons inthe network. In this figure we also show with black hor-izontal lines the values of N at which various k -cores ofthe network vanish. The positions of the tongues of ex-tra stability of the oscillating TMA phase appear to bebounded by these k -core transitions, suggesting that thedisappearance high- k k -cores changes the stability of theoscillatory (TMA) phase. We return to this point in sec-tion V, where we show that k -cores play a dominant role FIG. 1. Dynamical phase diagram of the model as a functionof the size of the network N and basal neuronal excitability∆ V . (A) An all-to-all coupled network with large g C = 3produces phase behavior consistent with mean-field predic-tions, but (B) sharp sigmoids (small g C = 0 .
01) produce adisordered diagram in which all dynamical phases are stronglymixed and the network dynamics is highly dependent on ini-tial conditions. Finally, in (C) randomly connected networks( p = 0 .
2) with large g C = 3, have initial-condition indepen-dent results with a modified dynamical phase diagram. Inall three panels the phases are: Q (blue), BTO (yellow), HA(dark green), ATO (light green), TMA (purple). All param-eter values are listed in appendix E. FIG. 2. Reentrant behavior along the of the TMA-HA phaseboundary. k -cores transition are shown as black lines, colorsare the same as at the previous figure. in the phase stability of a somewhat simplified version ofthe model.The fact that removing neurons from the network canenhance its ability to maintain stable oscillations seemsto be counter-intuitive. This reentrant behavior appearsat many phase boundaries in the system, including theone between the high-activity (HA) and quiescent (Q)phases. An example of such reentrant behavior at thisphase boundary is shown in Fig. 3. To understand howthis behavior emerges, it is simpler to study this casewhere the neurons’ dynamics reaches a fixed point ratherthan a limit cycle. Consider a fixed point of the sys-tem; setting the time derivatives on the left hand side ofEqs. 1, 2, we obtain V i = V eq + ∆ V ( C i ) τ V (cid:88) j M ij r ( V j ) (6) C i = C eq + ∆ Cτ C (cid:88) j M ij r ( V j ) (7)For neuron i to be rapidly firing, it must receive a num-ber of EPSPs consistent with both V i > V ∗ and C i < C ∗ .In this way, its somatic voltage is maintained above thethreshold and it remains sensitive to EPSPs. Too manyEPSPs will drive C i > C ∗ , resulting in the neuron’s so-matic potential falling below that threshold, while toofew EPSPs will allow V i < V ∗ even while maintainingdendritic sensitivity. As a result, the stable configura-tion of V i > V ∗ and C i < C ∗ can be destroyed by eitheradding or removing neurons that synapse on neuron i .We can see precisely how this works in an example of asmall network of seventeen neurons poised near the HA-Q boundary.In Fig. 4 we see that eliminating a low-firing rate neu-ron (number 16) from the network causes neuron 11 tochange from high to low firing rate. Removing an ex-citatory neuron has the expected behavior of reducingthe total activity of the network. But the subsequentremoval of another low-firing rate neuron (neuron 15) re- FIG. 3. Phase diagram showing reentrant behavior at the Q(blue) HA (dark green) phase boundary. There are also smallregions of the oscillatory phases: TMA (purple), ATO (lightgreen), BTO (yellow). sults in neurons 7 and 11 once again returning to highfiring rate. Finally, by removing the low-firing rate neu-ron 14, the entire network collapses into the quiescentstate.
III. SPONTANEOUS SYMMETRY BREAKINGON ALL-TO-ALL NETWORKS
In this section we explore phase separation on all-to-allcoupled networks, i.e., those having an adjacency matrixof the form M ij = 1 for all i (cid:54) = j and M ii = 0 for all i . The steady-state of the system spontaneously breaksthe permutation symmetry of the neurons. To explorethis symmetry breaking, we first investigate the symme-try preserving solution, that is obtained from the pair ofdifferential equations dVdt = 1 τ V ( V eq − V ) + ( N − V ( C ) r ( V ) (8) dCdt = 1 τ C ( C eq − C ) + ( N − Cr ( V ) , (9)which results from setting C i = C ( t ) and V i = V ( t ) forall i in Eqs. 1, 2 and using the all-to-all adjacency matrix.We demonstrated numerically that the dynamics ofthe full system Eqs. 1, 2 evolves towards this permu-tation symmetric solution for arbitrary initial conditionsif the sigmoidal functions P ( V ) and ∆ V ( C ) are smoothenough. If, on the other hand, these sigmoids are sharper,the system is unstable towards activity phase separation(breaking the original permutation symmetry of the un-derlying network) into time-independent subnetworks ofhigh and low firing rate neurons when the initial con-ditions are not themselves identical across the network.This symmetry broken state is, of course, not capturedby the above mean field analysis. To explore it, we needto analyze the full system of equations. Defining the sum of firing rates over the entire networkas R = (cid:80) i r ( V i ), we rewrite the dynamical system as dV i dt = 1 τ V ( V eq − V i ) + ∆ V ( C i ) [ R − r ( V i )] (10) dC i dt = 1 τ C ( C eq − C i ) + ∆ C [ R − r ( V i )] . (11)Now we can look for a self-consistent solution of this sys-tem, i.e. , we find V i ( R ) such that R = (cid:80) i r ( V i ). Study-ing the nullclines of Eqs. 10, 11, we see that it can haveone to three fixed points. The cases of one and threefixed points are shown in Fig. 5, where we see the inter-sections of the nullclines of Eqs. 10 and 11 in orange andblue respectively.Since we are looking for a self-consistent solution, wecan not analyze its stability directly from the graph;however, we can find the fixed points and later analyzetheir stability. For neuronal parameters consistent withsmooth sigmoids, there is only one fixed point V f , C f for afixed value of R . This self-consistent solution is both per-mutation symmetric and consistent with our mean fieldprediction. In contrast, for sharp sigmoids, there is morethan one fixed point, so it is possible to find some frac-tion of the network neurons at a high-voltage fixed point,while the remainder is at a low-voltage fixed point. Thenumber of neurons in these two categories is determinedby the condition R = (cid:80) i r ( V i ). There is also a rangeof parameters with g V > , g C = 0, such that the self-consistent solution does not exist.Therefore, there is nofixed point, and only oscillations are allowed. See Ap-pendix D for the details of the analytical calculation.For the small values of g C (sharp sigmoid), the phaseseparation into firing and quiescent neurons is the onlystable state. For intermediate values of g C we still ob-serve this stable state in Fig. 6A. We also show activityseparation into oscillating and quiescent subnetworks inFig. 6B. Continuing to increase g C we obtain synchronousoscillation of the whole network in Fig. 6C and, finally,fixed point with constant uniform activity in Fig. 6D. A. Step function limit: All-to-all networks
To better understand activity phase separation on thenetwork, it is useful to consider a non-physiological limitof the model in which the sigmoidal functions describ-ing both the firing rate and the dendritic adaptationare taken to be infinitely sharp, i.e. , step functions: g V = g C = 0. In this case, neurons with above-thresholdvoltage V ∗ fire at the maximal rate r m , while neuronsbelow that threshold voltage fire at the basal rate r b . Ifthe number of high and low firing rate neurons are n h and n l respectively ( n h + n l = N ) we find that V h = V eq + ∆ V ( C h ) τ V [( n h − r m + n l r b ] (12) C h = C eq + ∆ Cτ C [( n h − r m + n l r b ] . (13) FIG. 4. Example of reentrant high activity on a small network. Red neurons have
V > V ∗ and yellow V < V ∗ . The neuronsare numbered and the last neuron in each network is removed when going A → B → C → D . With the removal of neuron16 (from A to B) , the somatic potential of the neuron 11 drops below the threshold, as it has insufficient voltage input, andthe average network voltage falls below V ∗ , too. Going from B to C neuron 15, which synapses to neuron 4, is removed, whatlowers its calcium concentration. As a consequence, the somatic potential of the neuron 4 increases as well as its firing rate,resulting in increasing of the firing rate and voltage input to the neuron 7. The somatic potential of the neuron 7 then goesabove the threshold, too. The increasing firing rate of neuron 4 also raises the somatic potential of neuron 0, which raisessomatic potential of neuron 13, which in turn raises it for the neuron 11. Although somatic potentials of neurons 0 and 13does not exceed V ∗ , for neuron 11 it does. As a net result, the average voltage of the network raises above V ∗ . Finally, whenneuron 14 is removed (from C to D), all neurons are deactivated and ∆ V must increase to restore high activity. Similarly, the low firing-rate neurons have V l = V eq + ∆ V ( C l ) τ V [( n l − r b + n h r m ] (14) C l = C eq + ∆ Cτ C [( n l − r b + n h r m ] . (15)One can check that rate of spikes received by a highfiring-rate neuron R high = [( n h − r m + n l r b ] is lessthan that received by a low firing-rate neuron R low =[( n l − r b + n h r m ]. However, the condition for beingat a high firing rate is V ∗ < V h and a low firing rateis V ∗ > V l . For these inequalities to hold simultane-ously with the result that R high < R low , one needs thehigh firing-rate neurons to be more sensitive to incomingspikes than the low firing-rate ones. Thus we concludethat this state requires C l > C ∗ > C h . From this conclu-sion, we find n l , the number of low-firing rate neurons, to be n l = (cid:22) ( N r m − r b )∆ Cτ C + C eq − C ∗ ∆ Cτ C ( r m − r b ) (cid:23) , (16)where we have introduced the floor function: (cid:98) x (cid:99) = theinteger part of the real number x .To observe the phase separated state, we require thatthe high-firing rate neurons remain sufficiently sensitiveto incoming spikes. The lower bound of their sensitivity∆ V ( C h ) is given by V ∗ − V eq < ∆ V ( C h ) τ V [( n h − r m + n l r b ] . (17)While the number of neurons at low firing-rate neu-rons n l is fixed by Eq. 16, the identity of these neurons FIG. 5. Nullclines of the all-to-all N = 10 network describedby Eqs. 10 (orange) and 11 (blue) in the text. There are eitherthree fixed points or one fixed point depending on parameters.Assuming R is constant (and not fixed self-consistently) twoof the fixed points annihilate in a standard pitchfork bifurca-tion [13]. (A) g V = 0 . g C = 0 .
3, three fixed points. (B) g V = 5 mV, g C = 3, one fixed point. is determined solely by the initial conditions on the all-to-all network. There are a large number n ! n l ! n h ! of other-wise identical fixed points that are related by permuta-tion symmetry of the network. If, however, the networkis more sparsely connected and thus does not have thispermutation symmetry, there are fewer fixed points, asis discussed in the following section. IV. SYMMETRY BREAKING ON SPARSENETWORKS
If we randomly remove edges from the all-to-all net-work, we break the permutation symmetry of the neu-rons, and produce an instance of a network selected fromensemble of ER networks with probability p < p , as shown in Fig. 7. Below p ≈ . V vs. C graph, shown in Fig. 8.While for the smooth sigmoids the most common casewhen oscillations occur is an unstable fixed point, it isnot possible when sigmoids are very sharp. Indeed, anyfixed point that not exactly on the threshold in this caseis stable, as shown in Appendix C. Therefore, the onlyopportunity for the oscillatory or approximately chaoticbehavior is the absence of fixed point. A. Oscillations on star networks
In order to understand how all fixed points can vanishfor sparser networks, we can consider the special case ofa star network, in which one central neuron is bidirec-tionally coupled to N − V (0) τ V r b < V ∗ < ∆ V (0) τ V r m . (18)Furthermore, we require that all the peripheral neuronsfiring together at their basal rate are able to excite thecentral neuron over threshold. If, however, the centralneuron’s dendritic calcium is above threshold, then allthe peripheral neurons firing at their maximal rate arecollectively insufficient to excite the central neuron:( N − V ( C > C ∗ ) τ V r m < V ∗ < ( N − V (0) τ V r b . (19)We also demand two conditions on the calcium threshold.First, a single neuron cannot fire rapidly enough to pushanother neuron’s dendritic calcium over threshold: C ∗ > ∆ Cτ C r m . (20) N − N − Cτ C r b < C ∗ < ( N − Cτ C r m , (21)but they cannot do so when they are all firing at theirbasal rate.By obeying all of the above inequalities, the systemcannot reach a fixed point. Instead the network withthese step function neurons oscillates. The central neu-ron excites the peripheral ones and then those neuronsdrive the central neuron’s calcium concentration above FIG. 6. Activity phase separation on all-to-all connected network of N = 10 neurons. The traces show somatic potential ofindividual neurons as a function of time. Corresponding values of g C and g V are shown on the panels. (A) One neuron is athigh voltage , nine are quiescent. (B) Two neurons oscillate, eight are quiescent. (C) Synchronous oscillations of all neurons.(D) All neurons at high voltage. threshold rendering it insensitive. As a result, the cen-tral neuron returns to its low firing state and then so dothe peripheral ones. At this point the cycle begins again.Recall, however, that the step-function neurons on theall-to-all coupled network do not oscillate, reach one ofmany fixed points characterized by activity phase sepa-ration. By breaking the permutation symmetry of thenetwork, the star network admits a new synchronous os-cillatory phase. This is reminiscent of an effect calledconverse symmetry breaking [8], where the necessary con-dition for synchronous activity of a coupled network ofoscillators is an asymmetry of this system. We observea similar stability of globally synchronous oscillations inrandom networks that break the permutation symmetrysuch as the ER graphs discussed above. V. THE EFFECT OF NETWORKHETEROGENEITY ON PHASE SEPARATION
We have established that the neuron model leadsgenerically to dynamical phase separation. On permu- tation symmetric all-to-all networks, this phase separa-tion is a form of a spontaneously broken symmetry, butit exists on more sparse networks too. This poses thequestion: how does the network topology modify phaseseparation? To address this, we consider another simpli-fied limit of the model by setting r b = 0 (i.e. a neuronis either firing or not), using a step function firing rate( g V → C ∗ → ∞ . In this limit, we find phase separationinto groups of firing and nonfiring neurons at the fixedpoint. We demonstrate that the cluster of firing neuronsis determined exactly by the k -cores of the network foran integer k determined by the remaining neuron param-eters.Consider n i actively firing input neurons synapsing onneuron i . From the fixed point condition ˙ V i = 0 we find V i τ V = n i ∆ V r m . (22)For the i th neuron to be part of the group of activelyfiring ones, V i > V ∗ , which implies that number n i of FIG. 7. Number of stable fixed points as a function of thenetwork connectivity probability p for N = 100 neurons. For p = 1 this number coincides with n ! n l ! n h ! , and rapidly falls toone or zero when p (cid:46) . firing inputs exceeds a lower bound: n i ≥ V ∗ τ V ∆ V r max . (23)We intend to relate the actively firing group with atopological feature of the network: a k -core. This struc-ture is defined to be the maximal subnetwork, such thatwithin it each neuron has k or more inputs from the otherneurons in that subnetwork. k -cores have been discussedin a variety of applications in neuroscience, bioinformat-ics, ecology, and the study of social networks [14–17]. Thecondition for a neuron to be in the actively firing group,Eq. 23, is equivalent to membership within a k -core withthe integer k given by k = (cid:22) V ∗ τ V ∆ V r max (cid:23) . (24)If a k -core with some k given by Eq. 24 is absent fromthe network, the dynamical system on that network re-laxes to the quiescent fixed point, but if such a k -core ispresent, the neurons making up the k -core become fixedin the HA phase. We note that for typical values of k ,most of the network will be part of that k -core [7] be-cause for k > k -core has a discontinuous jump from zero to a signif-icant value as a function of the density of synaptic con-nections (in the thermodynamic limit of large networks).We compare our predicted phase boundary of the sys-tem and k -cores in Fig. 10. The predicted k value for aparticular set of n and ∆ V parameters corresponds ex-actly to the point in the phase space where the HA phasegives way to the Q phase. For this restricted version ofthe Feldman Del Negro model, at least, there is a precisecorrespondence between the neuronal network’s dynami-cal phase behavior and the prediction made purely fromthe topology of the underlying network. k -cores com-pletely determine the dynamical phase transition of the FIG. 8. Phase trajectories of the network with step func-tions in the mean V − C space. (A) Almost chaotic behav-ior. True chaos is not observed since the number of possiblestates is finite, but the voltage varies wildly. (B) Limit cy-cle with self-intersections, indicating asynchronous firing. (C)Standard limit cycle with synchronous firing, correspondingto true metronomic activity (TMA), rarely observed for thecase of step function limit. neurons interacting on them, and how the network phaseseparates into groups of high and low activity neurons. VI. DISCUSSION
We have explored the FDN model of oscillations in thepreB¨otC and found a form of dynamical phase separa-0
FIG. 9. A star network with N = 9 neurons. The peripheralneurons are bidirectionally coupled to the central neuron, butnot to each other.FIG. 10. The phase diagram for the simplified model dis-cussed in section V. There is no oscillatory phase, only quies-cent (Q, blue) and high activity (HA, green). Black horizontallines correspond to k -core transitions. We see almost exactcorrespondence between k -cores transtions and steps on thephase boundary. Small deviations are due to the fact that theaverage voltage of the whole network can be below V ∗ evenin the presence of the active k -core due to the averaging overall neurons including quiescent ones. tion on the network in which groups of neurons separateinto high and low firing-rate fixed points. This firing-quiescent phase separation plays the crucial role in thetermination of the TMA phase. One feature that emergesfrom this work is that the permutation-symmetric system(the all-to-all coupled network) admits a type of sponta-neous symmetry breaking into these high and low activ-ity phases. In more sparse, and physiologically relevantnetworks, the permutation symmetry is broken. In thiscase, the details of the network connectivity modify theinherent instability of the system toward phase separa- tion into groups of high and low firing-rate neurons. Inone particular limit of the model, we found that this in-teraction of neuronal dynamics and network topology isparticularly simple. By examining only the k -core struc-ture of the network, one can precisely predict both thedynamical phase diagram and which neurons will end upin the high and low firing rate groups.In the full model, the effect of the k -cores in deter-mining the phase boundaries of the dynamical systemremains, but no longer does it completely control thesedynamics. The incomplete influence of k -cores was ob-served earlier [5]. Here we believe we have better eluci-dated the underlying mechanism and explained why theircontrol of the dynamics is not complete.To assess the importance of these observations for thephysiological preB¨otC, we first present the numericallycomputed phase diagram for 1000 neurons using param-eters consistent with physiological measurements. Thisis shown in Fig. 12. Please see appendix A for a discus-sion of how the neuronal and network parameters wereselected.As discussed in the appendix, there is a remaining un-certainty in determining the value of ∆ C . Moreover,the full preB¨otC has somewhere between two to threetimes as many neurons as used in the simulation. Wenote, however, a scaling argument, based on the mean-field analysis of the model, that allows us to shift ∆ C as a way of effectively changing the network’s size. Inthe mean-field theory, three parameters ∆ C , ∆ V , and pN appear in only two combinations pN ∆ V and pN ∆ C .As a result, if we change ∆ C → ∆ Cλ , ∆ V → ∆ Vλ , and FIG. 11. Phase diagram of large networks with N up to 1000.All five phases are present: true metronomic activity (TMA)is purple, below threshold oscillations (BTO, yellow), abovethreshold oscillations (ATO, light green), high activity (HA,dark green), and quiescent (Q, blue). The pattern is approx-imately the same as in figure 1C, supporting the scaling ar-gument. pN → pN λ the mean-field solutions are invariant. Wecan test this scaling hypothesis in the full model by com-paring the phase diagram of the N=1000 network with∆ C = 2 . × − , shown in Fig. 11, to a much smaller1network of N = 100 and ∆ C = 0 .
1, shown in Fig. 1C.Their correspondence supports out exploration of largernetworks using calcium scaling.The scaling hypothesis suggests that, if we were able toexpand the network size used in Fig. 12 to the preB¨otC’strue physiological size, we would find that the region ofstable oscillations is bounded from above as well as forhigh and low neuronal excitability. We see for Fig. 11that the large N network has a rough phase boundarybetween the TMA and Q on the right side of the boundedTMA domain, which is incompatible with the mean-fieldanalysis and reflects the role of dynamical phase sepa-ration on the network. This result makes an interestingprediction in that there are regions of the phase diagramwhere increasing neuronal excitability can actually pro-duce globally quiescent networks, through calcium inhi-bition. This feature cannot be reproduced by the mean-field model.Another consequence of this large N phase diagram isthat we can predict the robustness of the network to dam-age. From Fig. 11, we see that under optimal conditions,one can destroy about eighty percent of the network be-fore causing the collapse of the oscillating phase. Thisagrees with experimental observations [11]. FIG. 12. Phase diagram of the network with physiologicallyrelevant parameters. It shows three stable dynamical phases:true metronomic activity (TMA, purple) consistent with thepreB¨oC’s physiological dynamics, as well as a high activity(HA, dark green), and a quiescent (Q blue) regime. Thereis a narrow band of above threshold oscillations (ATO, lightgreen).
We propose three types of experimental tests of theabove analysis. The first of these, alluded to above, isthat the network should be able to be silenced by in-creasing neuronal excitability. Secondly we predict thatthe roughness of the phase boundaries, particularly when N is large, suggests presence of multiple reentrant tran-sitions in which the network goes from being oscillatoryto quiescent, and back to oscillatory as neurons are re-moved from it. Third, one should be able to directlyobserve dynamical phase separation in the system. Ineither the high activity or quiescent phase, one should be able to find neurons trapped at the other fixed point sothat the globally quiescent state of the network shouldharbor some fixed fraction of high firing rate neurons.Conversely, the network in its globally highly active stateshould contain a subpopulation of neurons trapped intheir low firing-rate state. ACKNOWLEDGMENTS
AJL and VMS acknowledge partial support from NSF-DMR-1709785. VMS acknowledges support from Peccei-Holmes Graduate Research Fellowship and the BhaumikInstitute for Theoretical Physics. AJL, VMS and MB arethankful to Jack Feldman, Robijn Bruinsma, and SufyanAshhad for fruitful discussions.
Appendix A: Determining physiological parametersfor the model
To prepare Fig. 12 we must address the current under-standing of the physiological parameters of the neuronsas well as the network connectivity parameters. Mostof them can be fixed using experimental data shown inTable I. The exception is the set of parameters relatedto the dendritic calcium concentration, that are not cur-rently as well known. We chose them to reproduce theobserved dynamics of the system.2
FIG. 13. (A) Phase diagram for the all-to-all coupled network,using arbitrary initial conditions and smooth sigmoids. It isidentical to (B), the mean-field phase diagram.
Specifically, we set τ C to reproduce the observed pe-riod of stable oscillation. Taking into account that, forus, the units for the calcium concentration are arbitraryas is the choice of the zero for that concentration, weare left with only two independent parameters: g C and∆ C . g C is chosen to be large enough that we havereproducible phase behavior, avoiding highly heteroge-neous and initial-condition dependent results as shown inFig. 1B. We also require it to be small enough to producea true threshold for calcium inactivation of the dendrite, i.e. , C eq − C ∗ g C >
1. The last parameter to be fixed is ∆ C .The proper choice of ∆ C is facilitated by scaling behaviorobserved in the mean-field approximation of the model.Indeed, consider Fig. 1C. To quantitatively fit the invitro data, the network must support stable oscillationswhen ∆ V ≈ . C = 0 . N ≈ neurons –see Figure 12. It is computationally difficult to studynetworks with N > N . If wechoose ∆ C = 0 . , p = 0 . Appendix B: Mean-field solution
The mean field solution is obtained by assuming thatsomatic potentials and dendritic calcium concentrationsof all the neurons are the same: V i = V, C i = C . An-other assumption is that the network is on average ho-mogeneously connected, i.e. each neuron on average has pN inputs and outputs, with N the total number of neu-rons and p the connection probability. Then Eqs. 1, 2are rewritten as pair of equations for V and CdVdt = 1 τ V ( V eq − V ) + ∆ V ( C ) p ( N − r ( V ) (B1) dCdt = 1 τ C ( C eq − C ) + ∆ Cp ( N − r ( V ) , (B2)The phase diagram for such system is shown inFig. 15B. One may see that it is identical to the phasediagram for the all-to-all coupled network in Fig. 15Ain case of smooth transition in dendritic sensitivity. Themean-field approximation is valid for large g V and g C butbreaks when they become smaller. More detailed resultsare shown in Figure 14. FIG. 14. The results of the simulations for different g V and g C . Blue points corresponds to the case when these results(phase diagrams) fit mean-field approximation, and red onesto the cases when they do not (we disstinguish it by visu-ally comparing corresponding phase diagrams). Transparentregions are guides for the eye. Appendix C: Stability of fixed points on sparsenetworks
For the smooth sigmoid neurons, we found in sec-tion III that there was only a single fixed point – see3
TABLE II. Parameters used in the simulation
Fig. g C g V , mV C ∗ V ∗ − V eq ,mV τ C ,ms τ V ,ms r m ,Hz r b ,Hz ∆ C ∆ V max ,mV p N
1A 3 5 5 15 500 10 75 5 0.1 0-100 1 2-201B 0.01 5 5 15 500 10 75 5 0.1 0-100 1 2-201C 3 5 5 15 500 10 75 5 0.1 0-100 0.2 2-1002 5 5 15 15 500 10 75 5 0.035 0-25 0.75 2-1003 0 0 15 15 500 10 75 5 0.1 5-30 0.5 2-204 0 0 15 15 500 10 75 5 0.1 10 0.5 6-95A 0.3 0.5 5 15 500 10 75 5 0.1 50 1 105B 3 5 5 15 500 10 75 5 0.1 50 1 106A 0.3 0.5 5 15 500 10 75 5 0.1 50 1 106B 0.3 0.5 5 15 500 10 75 5 0.1 50 1 106C 1.1 0.1 5 15 500 10 75 5 0.1 50 1 106D 10.8 1.8 5 15 500 10 75 5 0.1 50 1 107 0 0 20 15 500 10 70 5 0.015 7.3 0.5-1 10010 0 0 ∞
15 500 10 70 5 0.1 1-5 0.5 1-5011 3 5 5 15 500 10 75 5 0.025 1-50 0.083 1-100012 3 5 5 15 500 20 40 0.1 0.015 1-10 0.065 1-100015A 3 5 5 15 500 10 75 5 0.1 0-100 1 2-2014A 0-3 0-20 5 15 500 10 75 5 0.1 0-100 1 2-20
Fig. 5. When that fixed point is unstable, the system ex-ecutes limit cycle oscillations. To determine the param-eter range of these oscillations, here we investigate thestability of that fixed point. Expanding the equations ofmotion near the fixed point { V fi , C fi } in the 2 N dimen-sional space of V i , C i , we define v i = V i − V fi , c i = C i − C fi and obtain dv i dt = − v i τ V + ∆ V (cid:48) ( C fi ) c i (cid:88) j M ij r ( V fj ) ++∆ V ( C fi ) (cid:88) j M ij r (cid:48) ( V fj ) v j (C1) dc i dt = − c i τ C + ∆ C (cid:88) j M ij r (cid:48) ( V fj ) v j (C2)Using r ( V ) and ∆ V ( C ) from Eqs. 3, 4, we find dv i dt = − v i τ V − g C ∆ V max σ (cid:48) ( C ∗ − C fi g C ) c i (cid:88) j M ij r ( V fj ) ++ 1 g V ( r m − r b )∆ V ( C fi ) (cid:88) j M ij σ (cid:48) ( V fj − V ∗ g V ) v j , (C3)and dc i dt = − c i τ C +∆ C ( r m − r b ) (cid:88) j M ij σ (cid:48) ( V fj − V ∗ g V ) v j . (C4)Dynamical phase separation requires that neither C ∗ − C fi nor V ∗ − V fi vanish. Therefore, if g V and g C are small,the terms with sigmoids are exponentially suppressed. Neglecting these we see that only the first term on theright hand side of Eqs. C3 and C4 remains. These implystability of the fixed point. For large g C and g V , however,we cannot ignore the terms proportional to σ (cid:48) , whichdestabilize the phase separated fixed point. Appendix D: Quasi-periodical phase diagrams
We consider the case g C = 0 , g V >
0, where quasi-periodical phase diagram can emerge. Following subsec-tion III A assuming the activity phase separation to n h neurons with V h , C h and n l neurons with V l , C l we writeequations for the fixed point V h − V eq = ∆ V ( C h ) τ V [( n h − r ( V h ) + n l r ( V l )](D1) C h − C eq = ∆ Cτ C [( n h − r ( V h ) + n l r ( V l )] . (D2)and V l − V eq = ∆ V ( C l ) τ V [( n l − r ( V l ) + n h r ( V h )](D3) C l − C eq = ∆ Cτ C [( n l − r ( V l ) + n h r ( V h )] . (D4)In case g C = 0 we have ∆ V ( C l ) = 0 , ∆ V ( C h ) =∆ V max (for C h < C ∗ , C l > C ∗ . Then in the same way asin the subsection III A we obtain the number of neuronsfiring at low rate n l = (cid:22) ( N r ( V h ) − r ( V l ))∆ Cτ C + C eq − C ∗ ∆ Cτ C ( r ( V h ) − r ( V l )) (cid:23) , (D5)4 FIG. 15. (A) Phase diagram for the all-to-all coupled network, g C = 0 , g V >
0. All five previously mentioned phases arepresent. There is a quasi-periodical pattern on the BTO-Q boundary (yellow and blue). This phase diagram fits thetheoretical prediction (B) where blue corresponds to the casethat Eq. D6 has a solution and yellow to the case that it doesnot. and simplify equations for the somatic potentials to get V l = V eq and V h − V eq ∆ V max τ V = N r ( V l ) − r ( V h ) + (cid:38) C ∗ ∆ Cτ C − ( N − r ( V l ) r ( V h ) − r ( V l ) (cid:39) (D6)where the ceil function (cid:100) x (cid:101) is the smallest integer that islarger or equal to x .Eq. D6 does not have a solution for V h for some pa-rameters. Due to the ceil function on the right hand side,the changes happen when its argument is incremented byone, what causes quasi-periodical structure of the phasediagram. Appendix E: Simulation details
Source code is available at https://github.com/mbibireata/Networks. For the simulation, we use theparameters from Table II. [1] C. Koch,
Biophysics of Computation: Information Pro-cessing in Single Neurons (Computational NeuroscienceSeries) (Oxford University Press, Inc., New York, NY,USA, 2004).[2] J. C. Smith, H. H. Ellenberger, K. Ballanyi, D. W.Richter, and J. L. Feldman, Science (New York, N.Y.) , 726 (1991), 1683005[pmid].[3] S. Coombes, P. Bressloff, and eds,
Bursting: The Genesisof Rhythm in the Nervous System (2005).[4] C. A. D. Negro, C. Morgado-Valle, and J. L. Feldman,Neuron , 821 (2002).[5] D. J. Schwab, R. F. Bruinsma, J. L. Feldman, and A. J.Levine, Phys. Rev. E , 051911 (2010).[6] C. Morgado-Valle, L. Beltran-Parrazal, M. DiFranco,J. L. Vergara, and J. L. Feldman, The Journal of physi-ology , 4531 (2008), 18635649[pmid].[7] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes,Phys. Rev. Lett. , 040601 (2006).[8] T. Nishikawa and A. E. Motter, Phys. Rev. Lett. ,114101 (2016). [9] J. C. Rekling, X. M. Shao, and J. L. Feldman, The Jour-nal of neuroscience : the official journal of the Society forNeuroscience , RC113 (2000), 11090613[pmid].[10] S. Ashhad and J. L. Feldman, bioRxiv (2019),10.1101/664946.[11] P. A. Gray, W. A. Janczewski, N. Mellen, D. R. Mc-Crimmon, and J. L. Feldman, Nature neuroscience ,927 (2001), 11528424[pmid].[12] P. Erd¨os and A. R´enyi, Publicationes Mathematicae De-brecen , 290 (1959).[13] S. H. Strogatz, Nonlinear Dynamics and Chaos: WithApplications to Physics, Biology, Chemistry and Engi-neering (Westview Press, 2000).[14] F. Morone, G. Del Ferraro, and H. A. Makse, NaturePhysics , 95 (2019).[15] G. D. Bader and C. W. Hogue, BMC Bioinformatics ,2 (2003).[16] S. B. Seidman, Social Networks , 269 (1983).[17] N. Lahav, B. Ksherim, E. Ben-Simon, A. Maron-Katz,R. Cohen, and S. Havlin, New Journal of Physics18