Interacting Turing-Hopf Instabilities Drive Symmetry-Breaking Transitions in a Mean-Field Model of the Cortex: A Mechanism for the Slow Oscillation

Electrical recordings of brain activity during the transition from wake to anesthetic coma show temporal and spectral alterations that are correlated with gross changes in the underlying brain state. Entry into anesthetic unconsciousness is signposted by the emergence of large, slow oscillations of electrical activity ( & 1 Hz) similar to the slow waves observed in natural sleep. Here we present a twodimensional mean-field model of the cortex in which slow spatiotemporal oscillations arise spontaneously through a Turing (spatial) symmetry-breaking bifurcation that is modulated by a Hopf (temporal) instability. In our model, populations of neurons are densely interlinked by chemical synapses, and by interneuronal gap junctions represented as an inhibitory diffusive coupling. To demonstrate cortical behavior over a wide range of distinct brain states, we explore model dynamics in the vicinity of a general-anesthetic-induced transition from ‘‘wake’’ to ‘‘coma.’’ In this region, the system is poised at a codimension-2 point where competing Turing and Hopf instabilities coexist. We model anesthesia as a moderate reduction in inhibitory diffusion, paired with an increase in inhibitory postsynaptic response, producing a coma state that is characterized by emergent low-frequency oscillations whose dynamics is chaotic in time and space. The effect of long-range axonal white-matter connectivity is probed with the inclusion of a single idealized point-to-point connection. We find that the additional excitation from the long-range connection can provoke seizurelike bursts of cortical activity when inhibitory diffusion is weak, but has little impact on an active cortex. Our proposed dynamic mechanism for the origin of anesthetic slow waves complements—and contrasts with—conventional explanations that require cyclic modulation of ion-channel conductances. We postulate that a similar bifurcation mechanism might underpin the slow waves of natural sleep and comment on the possible consequences of chaotic dynamics for memory processing and learning.


I. INTRODUCTION
Over the past two decades, extensive imaging studies using functional magnetic resonance imaging (fMRI) and positron emission tomography (PET) have demonstrated that, during specific cognitive tasks, human subjects exhibit a high degree of spatial organization in neuronal activation [1].The underlying basis of this spatiotemporal patterning of neural activity has not yet been unambiguously established.While there is evidence linking imaging patterns retrieved from blood oxygen-level-dependent (BOLD) studies to hard-wired axonal pathways defined by specific anatomic projections [2,3], known connectivities cannot account for many of the correlated BOLD observations [4][5][6].This absence of a clear mapping between anatomical connectivity and cognitive function suggests that some kind of spatial self-organization could be occurring.It is our belief that aspects of the observed patterns of neural activation can be generated by intrinsic Turing and Hopf instabilities that emerge spontaneously under physiologically realistic modulation of cortical parameters.Indeed, we advance the idea that normal brain function requires a delicate balance between these instabilities.
Our cortical model is based on early work by Liley et al. [7,8], but with enhancements motivated by Rennie et al. [9] and Robinson et al. [10], and recently augmented with the inclusion of electrical gap-junction synapses [11,12] to supplement communication via standard chemical synapses.Our state parameter is the mean excitatory firing rate, presumed to provide a proxy for the scalp-measured EEG.The model equations are set out in Sec.II; in summary, the cortex is pictured as a continuum of excitable tissue containing populations of excitatory neurons and inhibitory neurons (interneurons) that communicate locally via chemical and electrical synapses (gap junctions), and over longer ranges via myelinated axons.The gap junctions that link neighboring interneurons form direct resistive connections that permit ionic currents to flow, tending to equalize neuron voltages.The bulk effect of electrical synapses is to produce diffusion terms similar in form to those found in standard reaction-diffusion models that support Turing structures [13].We find that increasing the strength of the interneuronal gap-junction coupling D 2 can precipitate a Turing bifurcation that generates stationary labyrinthine patterns of raised and lowered cortical activity [11].
Spatial self-organization associated with Turing structures has been identified in a variety of physical systems and studied extensively in the context of the so-called Brusselator equations, first expounded by Prigogine and coworkers [14] and expanded subsequently by many other authors (e.g., [15,16]).For a Turing instability to occur, the inhibitor must diffuse more rapidly than the activator.Since the cortex can be modeled as a network of inhibitory and excitatory neural populations coupled via nonlinear interactions, it has the necessary prerequisites for emergence of spatial instability.The prediction of such Turing structures in a mean-field cortex has been made by several investigators [17][18][19][20][21][22], with most work being confined to the consideration of a one-dimensional (1D) cortex.Typically, the Turing mechanism in these cases has relied on a presumed ''Mexican-hat'' kernel featuring long-range inhibitory and short-range excitatory axonal connections.In contrast, in the present paper we present an alternative mechanism for Turing instability in a two-dimensional (2D) model of the cortex in which the longer-range inhibitory dominance arises from gap-junction currents between adjoining interneurons that may form a diffusive syncytium scaffolding the cortex [23].If inhibitory gap-junction diffusive coupling D 2 is sufficiently strong, a spatial instability is predicted.
Our model also predicts that a sufficient reduction in i , the inhibitory rate constant for postsynaptic response at chemical synapses, can promote a temporal (Hopf) bifurcation that manifests as global, whole-of-cortex oscillations of soma voltage; this form of oscillatory dynamics has been interpreted as grand mal seizure [24][25][26].
Significantly, a simultaneous tuning of i and D 2 can locate a codimension-2 Turing-Hopf (c2TH) point [27] at which Turing and Hopf instabilities interact, giving rise to complex mixtures of traveling waves and oscillating Turing structures [28,29].Such traveling patterns allow separated patches of cortex to exhibit phase-coherent oscillations, providing a plausible dynamic substrate for the ''binding'' and ''action at a distance'' behaviors supposed to be prerequisites for cognition and consciousness [30].
We have argued that the small-amplitude ultralowfrequency spatiotemporal oscillations generated near the Turing-Hopf point might represent some features of the default noncognitive idling background state of the conscious brain since these oscillations can provide a functional connection between brain areas that lack an identifiable anatomical connection; these coherent oscillations persist over time scales ranging from hundreds of milliseconds to seconds [29].Empirical and modeling studies of the resting default state by Honey et al. [31,32] show that structurally connected regions of cortex do exhibit stronger functional connectivity, but only over time scales of minutes.This suggests that, while brain dynamics is guided by anatomical connectivity, it is not necessarily constrained by it, particularly at the faster time scales (of the order of seconds) where c2TH activity patterns can emerge.
An important difference between this and previous iterations of our model [11,12,28,29] is the inclusion of a long-range axonal fiber that connects two points in the 2D neural field.While we recognize that other workers [31,33,34] have incorporated into their models anatomically detailed maps of the cortico-cortical links between brain areas, our intention here is simply to identify any new dynamics arising from a single heterogeneous connection.Our approach is a generalization of that described by Jirsa and Kelso for a 1D neural topology [21].
Locating the cortex close to the c2TH point gives it access to a rich variety of diverse spatiotemporal behaviors.To explore the alterations in dynamics across a range of cognitive states, we adopt the anesthetic phase transition as a paradigm.This framework, well studied by the authors, allows us to move with relative ease between the putative cortical states with minimal tuning of parameters.The system is characterized by a distribution of equilibrium states which can be used to identify specific phases of consciousness.
We select two reference states, ''wake'' and ''anesthetic coma,'' and carefully investigate the dynamic impact of alterations in Turing-Hopf balance by varying the gapjunction connection strength D 2 .For coma, the dynamics can range from a wakelike state that might correspond to delirium [35] (large D 2 ) to seizurelike bursting and eventually suppression (small D 2 ), while for wake, the progression is from wake to ''seizure.''In both cases, at intermediate values of diffusion, we see spontaneous emergence of a large-amplitude slow chaotic oscillation between up (activated) and down (quiescent) states.
In our previous study [12]-in which we examined the role of gap-junction modulation of Turing-Hopf dynamics with respect to transition from wake to seizure-we interpreted the emergent chaotic oscillations as an incoherent precursor signaling the onset of the pathologically coherent seizure state (see Fig. 2 of Ref. [12]).But when the anesthetic effect was boosted (Fig. 3 of Ref. [12]), downscaling of inhibitory diffusion gave a chaotic regime, closely followed by cortical collapse into a completely silent ''suppressed'' state.The addition of the long-range two-point connection described here (see Fig. 1) allows a new regime of ''anesthetized'' dynamics to emerge: The chaotic slow oscillations now give way to a seizurelike ''bursting'' state as diffusion is reduced (see Figs. 8 and 10), followed eventually by full suppression.Although highly idealized, the fact that the two-point connection evokes a more realistic anesthetic dynamics motivated us to identify the chaotic state as a potential source of the cortical slow oscillation, and to examine in detail the impact of the longrange connection on anesthetic cortical dynamics (Fig. 12).The interesting finding is that, in a mean-field sense, the slow oscillation can arise without recourse to the classical mechanisms involving ion channels in neurons.
The paper is structured as follows.In Sec.II, we define the cortical model, locate its equilibrium states, and explain how cortical dynamics can be quantified in terms of a phase coherence measure.In Sec.III, we report on a comprehensive series of numerical simulations that survey the effect of inhibitory diffusion on wake and anesthetic-coma states, and assess the impact of the two-point connections on system behavior.A significant conclusion of this theoretical study is the prediction that the slow-wave up-down oscillations of anesthetic coma are chaotic in nature, arising spontaneously from a competitive interaction between Turing (spatial) and Hopf (temporal) instabilities.This mechanism for the slow oscillation is quite distinct from the conventional view of a cyclic alternation in ionic concentrations that tends to suppress firing in the active state, then activates firing in the quiescent state.Our prediction of chaotic dynamics is supported by clinical studies of phase synchronization changes in EEG during induction of propofol anesthesia [36].The Turing-Hopf interaction may also underlie the slow waves of natural sleep, and, if so, it has interesting implications for memory processing during sleep.

A. Model equations
Our model envisions the cortex as a 2D network of excitatory and inhibitory neurons that are interconnected locally via resistive gap junctions and neurotransmittermediated chemical synapses, and over distances via longrange myelinated axons.Since scalp and cortical electrodes can only sense average voltages generated by populations of neurons (and not single-unit activity), we represent the cortex as a continuum of excitable tissue whose cortical parameters have been coarse grained over a spatial extent of order 1 mm 2 , corresponding to the area of a cortical macrocolumn [37].In our numerical simulations, we map the cortex to a 120 Â 120 grid (see Fig. 1).
The spatially averaged excitatory and inhibitory soma potentials V e and V i at grid location r ¼ ðx; yÞ obey the following differential equations: where we have adopted a left-to-right convention for labeling presynaptic and postsynaptic connectivity so that, for example, eb is given as ''e !b,'' with the postsynaptic subscript b standing for either e or i.The terms ½. . . in square brackets are chemical-synaptic voltage inputs.The nabla-square symbol denotes the 2D Laplacian operator The final term on the right of Eq. ( 1) is the voltage contribution arising from gap-junction currents diffusing between adjacent neurons.Gap junctions are connexinprotein channels that form spontaneously at the points of contact between the dendrites of adjoining neurons, and which allow the passive diffusion of ions driven by intercellular potential differences [38].We write D bb as the diffusive-coupling strength between electrically adjoined i $ i or e $ e neuron pairs, so D bb = b can be regarded as an effective diffusion coefficient (see [11] for a detailed derivation and magnitude estimates).To simplify notation, we now write ðD 1 ; D 2 Þ ðD ee ; D ii Þ.Because gap-junction coupling between inhibitory interneurons is substantially more abundant than that between excitatory neurons [38], we set D 1 to be a small fraction of D 2 (D 1 ¼ D 2 =100), with D 2 being an adjustable parameter.Layer 1 of the cortex is predominantly (over 90%) comprised of inhibitory neurons [23], so it is not implausible to posit a continuous syncytium of interneuron-to-interneuron diffusive scaffolding spanning the cortex.The physiological significance of dominant inhibitory diffusion is that it allows for the spontaneous formation of distributed Turing patterns that might modulate cortical activity.
The b in Eq. ( 1) are soma time constants for the neuron populations, and V rest b are their respective resting voltages.The b are the chemical-synaptic strengths given by the area under the excitatory or inhibitory postsynaptic potential (PSP), with e > 0 (EPSP) and i < 0 (IPSP).These strengths are scaled by dimensionless reversal-potential functions c ab , The diagram shows direct axonal connections from coordinate ðx 1 ; y 1 Þ ¼ ð20; 60Þ to ðx 2 ; y 2 Þ ¼ ð40; 60Þ (black arrow), and back (gray arrow).Communication between connected points occurs after a time delay Át ¼ Ár=v, where v is the axonal conduction speed and Ár is the physical separation ðL=N x Þjx 1 À x 2 j, with L ¼ 25 cm and N x ¼ 120.We apply periodic (toroidal) boundaries.
that are normalized to unity when the neuron is at its resting voltage and are zero when the membrane voltage reaches the relevant reversal potential, taken to be V rev e ¼ 0 mV for excitatory (AMPA) receptors and V rev i ¼ À70 mV for inhibitory (GABA) receptors.For simplicity, we assume that the V rest b term appearing in the denominator of c ab is a constant, independent of the ÁV rest b offset of Eq. ( 1).The four È ab functions in Eq. ( 1) are postsynaptic fluxes.Assuming that the impulse response at the postsynaptic dendrite obeys an alpha-function form 2 The and superscripts distinguish longer-range (cortico-cortical) versus local chemical-synaptic inputs; N eb , N eb are the number of such input connections, and eb , Q e;i the corresponding long-range and local spike-rate fluxes.sc eb is the unstructured stimulatory tone arriving from subcortical sources.This stimulus is modeled as a low-intensity spatiotemporal white-noise variation eb about a constant tonic background h sc eb i, with a being a dimensionless amplitude scale factor.The eb ð r; tÞ is Gaussian-distributed delta-correlated noise with statistics hð r; tÞi ¼ 0 and h m ðt 1 Þ n ðt 2 Þi ¼ m;n ðt 1 À t 2 Þ.In numerical simulation with fixed time increment Át, we approximate continuous noise ðtÞ with discrete noise samples f k g constructed using the MATLAB randn random number generator: k ¼ randn= ffiffiffiffiffi ffi Át p .Inclusion of a noisy stimulus encourages the cortical model to explore its repertoire of dynamic behaviors accessible from a resting, equilibrium state, and represents the reality that any living, biological system is unavoidably exposed to a continuous background wash of random perturbations.The fourth flux term on the right of Eq. (2), ;het eb , is the symmetrybreaking heterogeneity arising from the direct two-point connection illustrated in Fig. 1  where Ã eb is the inverse-length scale for e !b axonal connections, and v is the axonal conduction speed.Jirsa and colleagues [21,[39][40][41] have investigated the dynamical changes in a 1D neural-field model wrought by breaking isotropic symmetry with a direct point-to-point connection.Here, we generalize their approach to two spatial dimensions as a modest first step towards representing realistic anatomical heterogeneity.Note that these earlier works are primarily concerned with the effects of point-to-point transmission delays on cortical stability, showing that higher axonal speeds-such as those that occur during developmental myelination-tend to stabilize the homogeneous steady state against oscillatory instabilities [40].In contrast, here we keep the two-point axonal delay fixed and vary the inhibitory diffusion strength as a control parameter to alter the relative balance between spatial and oscillatory instabilities that can emerge in a 2D model cortex.
The heterogeneous flux input ;het eb in Eq. ( 2) is expressed in terms of a lossless two-point connectivity kernel [21], where Q e ð R; TÞ is the excitatory neural activity at location R at earlier time T ¼ t À jr À Rj=v.This timing offset accounts for the propagation delay jr À Rj=v from source point R to field point r with axonal speed v. Here, 'm is the connection strength from point ' to point m.Applying the properties of delta functions, Eq. ( 5) becomes with offset times T n ¼ t À jr À rn j=v for n ¼ 1, 2. This means that, in numerical simulations, there can be no flux contribution from either distant connection until a delay time jr 1 À r2 j=v has elapsed.Figure 1 illustrates the geometry of the direct two-point connections used in the present work.For simplicity, we place both connected points along the y ¼ 60 column of our 120 Â 120 cortical grid and impose balanced twoway connections of equal strength in both directions: To investigate the effect of variable separation distance Áx ¼ jx 1 À x 2 j, we fix x 1 ¼ 20 and allow x 2 to vary.For most of the work reported here, we set x 2 ¼ 60 to give a physical separation of Ár ¼ 8:33 cm (i.e., one-third of the side length of our cortical grid).Note that the values for our model parameters have been chosen to be physiologically plausible; see Table I for parameter values.(See Ref. [42].)

B. Modeling anesthesia
At the molecular level, inductive anesthetic agents suppress neural activity by prolonging the opening of gamma-aminobutyric acid (GABA) channels on the postsynaptic neuron [43], hyperpolarizing the neuron by increasing the net influx of chloride ions.For example, the anesthetic propofol prolongs the duration of the inhibitory postsynaptic potential (IPSP) without altering its peak amplitude [44].We model propofol's effect by scaling both the inhibitory rate constant i and synaptic strength i by a dimensionless scale factor that is set to unity in the absence of propofol, and which grows proportionately to propofol concentration, where 0 i and i are the anesthetic-free default values.This scaling ensures that the area of the IPSP response (representing the total charge transfer) increases linearly with drug concentration while the IPSP peak remains unchanged.We note that at very high propofol concentrations-well above the clinically relevant range-the charge-transfer versus drug-concentration curve shows saturation effects [44], but the assumption of linearity is accurate at low concentrations and has been used by Hutt and Longtin [45] in their anesthesia modeling.

C. Equilibrium states of the cortex
We view cortical dynamics as a spatiotemporal variation about, or relaxation towards, the family of underlying homogeneous equilibria that define one or more resting reference states of cortical activity.By linearizing about such a reference state, we can apply linear stability analysis to make predictions about the conditions under which spatial or temporal instabilities might arise [12].If the homogeneous equilibrium is stable, the cortex will tend to relax to this steady state when input stimuli are removed; if the equilibrium is unstable, it acts to organize the resulting spatiotemporal dynamics.And if the cortex has access to multiple steady states, very rich dynamical behaviors can emerge, particularly when more than one of these states is unstable.
We locate the homogeneous equilibrium states by setting the heterogeneous flux term in Eq. (2) to zero ( ;het eb ¼ 0) and eliminating all space and time derivatives in differential equations ( 1)-(4) (r 2 ¼ 0; @=@t ¼ @ 2 =@t 2 ¼ 0), then solving (numerically) the resulting set of nonlinear coupled algebraic equations for the steadystate firing rates ðQ e ; Q i Þ of the excitatory and inhibitory neural populations.For the present work (see Fig. 2), we have visualized the distribution of steady-state solutions as a (multivalued) function across a domain defined by the anesthetic effect and resting potential offset ÁV rest e .We note that, over part of the ð; ÁV rest e Þ domain, the steady-state manifold folds back on itself to form a reverse S-shaped reentrance whose top and bottom surfaces correspond to activated and suppressed neural states, respectively.Thus, we identify the upper surface with awake and the lower surface with anesthetic coma.In Fig. 3, we have plotted a vertical slice through the Fig. 2 manifold in the Q e À plane at resting voltage offset ÁV rest e ¼ 1:5 mV.In the work that follows, we investigate the impact of two-point heterogeneity on the dynamical behavior of the awake and comatose cortexes.We acknowledge that, as soon as spatial heterogeneity enters, the equilibrium states will also become spatially dependent, meaning that each spatial coordinate on the 2D cortical grid could own a distinct equilibrium manifold.Nevertheless, it is conceptually useful to regard the homogeneous distribution of steady states as our background reference, which is to be perturbed by the ;het eb heterogeneity.This is not unreasonable, given that the nonlocal term has a strength ( eb ¼ 200) equivalent to a quarter of the local connectivity (N eb ¼ 800) and a tenth of the longer-range connectivity (N eb ¼ 2000; see Table I), so its effect is to raise destination firing rates by only a few percent above the homogeneous steady-state values.

D. Linear stability predictions
Equations ( 1)-( 4) define the cortical model in terms of two first-order (V e;i soma voltages) and six second-order (È ee;ei ; È ie;ii ; ee;ei firing-rate fluxes) partial differential equations (DEs), equivalent to 14 first-order DEs.In Sec.IV, we examine model response to changes in both inhibitory diffusion (D 2 ) and inhibitory response rate ( i ) by time-stepping these equations to obtain numerical solutions.Since the long-range axonal strength has been set deliberately weak, we can make reasonable predictions as to which instability modes might emerge by looking at the linear stability of the homogeneous stationary states.
We do this by disabling the subcortical noise and disconnecting two-point axonal linkages ( ¼ 0).Noting the parameter symmetries evident in Table I ), the homogeneous system reduces to a set of two first-order and three second-order DEs, equivalent to eight first-order equations.We define the eight-variable state vector X ¼ ½V e ; V i ; È eb ; _ È eb ; È ib ; _ È ib ; eb ; _ eb T with an equilibrium value Xð0Þ .We then express the cortical state as its equilibrium value plus a small spatiotemporal disturbance, Xðt; rÞ ¼ Xð0Þ þ Xðt; rÞ; with Xðt; rÞ being a plane-wave perturbation Xðt; rÞ ¼ XðtÞe i qÁr ¼ Xð0Þe Ãt e i qÁ r; where q is the wave vector with wave number j qj ¼ q, and Ã is an eigenvalue whose real part gives the growth rate of the Xð0Þ initial perturbation: If ReðÃÞ > 0, an instability is predicted.Substituting X ¼ Xð0Þ þ X into Eqs.( 1)-( 4) and retaining only linear terms results in the matrix equation where J is an 8 Â 8 Jacobian matrix in which the r 2 Laplacians for excitatory and inhibitory diffusion [Eq.( 1)] and wave propagation [Eq.( 4)] appear as Àq 2 terms.The eight eigenvalues owned by J describe the linearized dynamics of the homogeneous cortex.For each wave number FIG. 3. Vertical slice through the homogeneous equilibrium manifold of Fig. 2 showing the distribution of steady-state firing rates as a function of the anesthetic effect for constant excitation ÁV rest ¼ 1:5 mV.We associate the upper, high-firing branch with the awake state, and the lower, low-firing branch with the ''asleep'' state corresponding to anesthetic-induced coma.Circled points indicate reference states at i ¼ 1:0 and 1.018 on awake and coma branches, respectively.q, we extract and plot the dominant eigenvalue-i.e., that eigenvalue whose real part is the most positive (or the least negative)-since this eigenvalue describes the most strongly growing (or most long-lived) mode at a given spatial frequency.
Figure 4 plots the predicted wave number q dependence of the real and imaginary parts of the dominant eigenvalue for three values of inhibitory diffusion (D 2 ¼ 0:7, 0.4, 0:1 cm 2 ) for the top-and bottom-branch equilibria at i ¼ 1:0. Figure 4(a) shows that, for all three diffusion values, the high-firing up-state is expected to destabilize in favor of a whole-of-cortex Hopf oscillation (peak instability occurs at q ¼ 0) of temporal frequency of approximately 3 Hz.In contrast, the low-firing down-state [Fig.4(b)] exhibits instability properties that are strongly D 2 dependent, starting from a Turing-dominated spatial mode of spatial frequency q % 0:4 waves/cm, wavelength around 2.5 cm, for strong inhibitory diffusion (D 2 ¼ 0:7 cm 2 ), evolving to interacting but damped Hopf and Turing instabilities at moderate diffusion (D 2 ¼ 0:4 cm 2 ), and giving way to a damped Hopf instability when diffusion is weak (D 2 ¼ 0:1 cm 2 ).
As we will see later in Sec.III, these linear stability predictions provide a reasonable prediction of spatial and temporal instabilities, but the range of complex and chaotic nonlinear dynamics that can emerge from the modal interplay of Hopf and Turing instabilities only becomes apparent when the full cortical equations are simulated numerically.

E. Phase coherence as a measure of cortical state
Although D 2 (inhibitory diffusion) and i (inhibitory synaptic rate constant) have no influence on the location and distribution of equilibrium states of our model cortex, these two inhibitory parameters are key determinants of the nonequilibrium dynamical states that can spontaneously emerge from the equilibrium manifold.These far-fromequilibrium states cover the gamut from static Turing patterns (large D 2 ) at one extreme, to global whole-of-cortex Hopf oscillations (small i ) at the other, and encompass a rich diversity of complex spatiotemporal dynamics in between.The intermediate settings for which neither instability dominates are particularly interesting, because minor alterations in the Turing-Hopf balance can produce very dramatic changes in the emergent cortical dynamics.
This dynamic complexity presents a challenge: How can we encapsulate the state of a noise-stimulated cortex whose oscillatory patterns of activity are fluctuating in time and space?We choose to compute the average phase coherence, R, as a measure of the degree to which oscillations at separated points on the cortical grid have a consistent phase relationship.Following Mormann et al. [46], we define the time-averaged phase coherence between reference point r0 and some other point r1 on the cortical grid as where the instantaneous phase angles ðt; r0 Þ, ðt; r1 Þ are computed from the Hilbert transforms of the excitatory firing rates Q e ðt; r0 Þ, Q e ðt; r1 Þ recorded at locations r0 , r1 at time t (see [12] for further details and for sample MATLAB code).The time average h. ..i is taken across a 2-s interval captured when the cortex has settled into a fully developed (and usually nonequilibrium) dynamical state, well after initial startup transients have decayed away.These startup transients typically die out within about 0.5 s for top-branch awake simulations [see strip charts of Fig. 5(b)] but can persist for up to 5 s for the bottom-branch coma state [Fig.8(b)].
If neural activity at r0 and r1 is tightly phase coupled, then mapping their sampled phase differences to the complex unit circle via Eq.( 7) will give a tightly clustered angular distribution of phasors leading to an enhanced mean coherence ðR % 1Þ.However, if Q e ð r0 Þ and Q e ðr 1 Þ have phase angles that are unrelated, their phase differences will tend to be randomly distributed around the unit circle, giving a degraded coherence statistic (R % 0).
Our numerical simulations have a grid resolution of 120 Â 120 elements, so, in principle, there are 120 2  2 % 10 8 pixel pairs for which phase coherence could be computed.However, to keep the calculation manageable, we choose instead to focus on the 120 elements making up the middle vertical column at y ¼ 60 (see Fig. 1), giving distinct pairs of coherence values that we visualize as Rðx 0 ; xÞ coherence maps (see Figs. 6 and 9).For a given simulation run, we assume that the spatial average across the coherence map provides a scalar estimate of the ''global'' phase coherence for the entire cortex; these global coherence values are particularly useful for tracking changes in Turing-Hopf interaction dynamics arising from variations in D 2 inhibitory diffusion (Figs. 7 and 12).

III. RESULTS
In Fig. 3, we identify two contrasting states for the homogeneous model cortex: an awake state on the upper branch of the equilibrium manifold and an anesthetically suppressed coma state on the lower branch at i ¼ 1:018.For both cortical states, we now investigate the effect of variation in the strength of inhibitory diffusion D 2 , and also the impact of direct axonal connections linking two separated points on the cortical grid.

A. Effect of inhibitory diffusion on resting wakefulness
We define resting wakefulness as a background idling state of the brain that spontaneously emerges in the absence of structured external inputs.We have suggested previously that the slow beating dynamics (& 0:1 Hz) observed in BOLD functional MRI recordings during default network activation [5] arises from a competitive interaction between Turing (spatial) and Hopf (temporal) instabilities: The inhibitory diffusion strength D 2 (via interneuronal gap junctions) mediates pattern formation, and the inhibitory rate constant i (governing inhibitory PSP charge transfer via chemical synapses) mediates network oscillation [29].Their coaction can generate smallamplitude oscillations that are phase synchronized across separated cortical regions, with the degree of phase synchrony being strongly dependent on D 2 .We selected i and D 2 as parameters of interest because variations in their values provide the simplest, physiologically justified way of generating the two types of instability; while these parameters have a controlling influence on model dynamics, they have no effect on the shape or distribution of the manifold of equilibrium states shown in Fig. 2. Figure 5 illustrates the dynamical effect of stepped reductions in gap-junction diffusion for a cortex initialized at the upper-branch wake state indicated in Figs. 2 and 3.The six rows correspond to six numerical experiments for different values of diffusion, and, for each experiment, the three columns show (a) a 4-s time-series extract of excitatory firing rates Q e ðtÞ for two representative points along the y ¼ 60 vertical midline of the 120 Â 120 cortical grid; (b) a 20-s space-time image Q e ðt; xÞ that captures the complete time evolution of cortical activity along the y ¼ 60 midline strip; and (c) a bird's-eye snapshot Q e ðy; xÞ of the state of the cortex at the conclusion of the 20-s simulation.
Scanning Fig. 5 from top to bottom, we observe that the cortical dynamics for the i ¼ 1:0 wake state is remarkably sensitive to changes in D 2 diffusion strength.From the linear stability predictions of Fig. 4, it is evident that this sensitivity arises from diffusion-dependent alterations in the balance between competing spatial (Turing) and temporal (Hopf) instabilities: for strong diffusion (e.g., D 2 * 0:7 cm 2 ), the Turing instability dominates, so spatial patterning takes precedence over temporal oscillations, while for weak diffusion (e.g., D 2 & 0:3 cm 2 ), the Hopf instability drives large-scale oscillations and traveling waves that erase any static Turing patterns.For intermediate diffusion strengths (0:4 & D 2 =cm 2 & 0:6), neither instability dominates, and the result is a turbulent mix of oscillating and translating neural structures whose boundaries and shapes evolve continuously in space and time.
Although the space-time images of Fig. 5 capture aspects of the spatiotemporal dynamics, the onset and time evolution of the neural activity structures is better visualized in a movie sequence.See Video 1 for a set of movies for three representative settings of inhibitory gap-junction diffusion.
We identify the top row of Fig. 5 (D 2 ¼ 0:7 cm 2 ) as the resting state of the awake brain.The cortex settles into an irregular mazelike labyrinthine Turing structure of elevated and suppressed activity, on which rides an approximately 4-Hz temporal oscillation that is not only phase correlated across separated Turing limbs, but also slowly modulated on time scales of 5-10 s.The spatiotemporal spectrum and phase coherence map for the noncognitivewake strip chart is shown in the top row of Fig. 6.The spatial spectrum in Fig. 6(a) shows peaks at wave number k % 10=L x waves/cm, consistent with the approximately 2.5-cm spacing between Turing limbs evident in the top row of Fig. 5(c); the coherence map of Fig. 6(b) confirms the claim that separated patches of cortex are oscillating in a phase-coherent manner.
Successive reductions in inhibitory diffusion strength tend to destabilize the Turing structures: Small-scale oscillations about the ''up'' (activated) and ''down'' (inactivated) Turing states [top row of Fig. 5(a)] are replaced by lower-frequency large-scale oscillations between up and down states (third and subsequent rows), with oscillations becoming progressively more synchronized in time and space, eventually leading to an extreme oscillatory state that we label as seizure (bottom two rows of Figs. 5 and 6).The passage from coherent wake (D 2 ¼ 0:7 cm 2 ) to coherent seizure (D 2 & 0:3 cm 2 ) is marked by the emergence of intermediate states whose patterns of neural activity are incoherent and turbulent.Following the method described by Destexhe [47], we ran a series of paired noise-free simulations (see Fig. S1   These model results allow two predictions.First, if the cortex is in the idling noncognitive-wake state, complete closure of inhibitory gap junctions should tend to promote seizure.Second, the passage from wake to seizure should be observable as an initial decrease, followed by increase, in phase coherence between separated cortical electrodes.

B. Effect of inhibitory diffusion on anesthetic coma
We now shift attention from the awake state to the coma state and ask the following question: What is the effect of variations in gap-junction diffusion for an anesthetized cortex?For our comatose state, we select a point on the lower branch of the homogeneous equilibrium manifold, just beyond the region of roots mapped in Figs. and 3, at coordinate ð i ;ÁV rest e Þ¼ð1:018;1:5mVÞ.Figures 8 and 9-respectively analogous to Figs. 5 and 6 for the wake state-illustrate the dynamical impact of stepped changes in gap-junction diffusion for an initially low-firing comatose cortex that includes the same two-point connectivity ( ¼ 200) and is perturbed continuously by the same intensity of spatiotemporal white noise as that used for the wake experiments described above.
If the diffusion is weak (D 2 & 0:1 cm 2 : bottom rows of Figs. 8 and 9), an initially silent, suppressed cortex remains electrically silent, and thus generates a ''flatline'' EEG trace.To ''rouse'' the cortex, it is sufficient to increase the diffusive-coupling strength.For example, when D 2 is raised to 0:3 cm 2 (fifth rows of Figs. 8 and 9), the two-point axonal connections at x 1;2 ¼ 20, 60 become spontaneous sources of coherent, periodic outgoing wave fronts that invade the entire cortex, producing synchronous patterns that may correspond to the bursting phase of burst suppression seen in very deep anesthesia.
Further increases in D 2 (to 0.4, 0.5, 0:6 cm 2 : rows 4, 3, 2 of Figs. 8 and 9) cause the seizurelike bursting patterns to break up into highly irregular, relatively incoherent, turbulent structures.This chaotic behavior arises from competitive interference between Hopf (temporal ordering) and Turing (spatial ordering) instabilities and may provide a natural generator for the irregular slow-wave cortical oscillations of inductive anesthesia (e.g., see Fig. 5E of [48]).Boosting D 2 to 0:8 cm 2 (top rows of Figs. 8 and 9) tilts the balance in favor of spatially ordered structures carrying small-amplitude oscillations that are strongly coherent across the cortex, and very similar to the noncognitive-wake state in Fig. 5 (top row).We might refer to this state as ''delirium'' since the is still the influence of anesthetic.Nevertheless, the notion that a sleeping cortex can be roused to a wakelike state by administering an agent (e.g., modafinil) that opens gap junctions does have clinical support [49].
The spectra and coherence maps of Fig. 9 provide additional insights into comatose dynamics.Although delirium (top row: D 2 ¼ 0:8 cm 2 ) and bursting (fifth row: D 2 ¼ 0:3 cm 2 ) both have elevated average phase coherence, their wave-number-frequency spectra and coherence maps are strikingly different.The delirium state is dominated by an oscillatory but nonpropagating Turing pattern, so the wave-number-vs-frequency dispersion trends are vertical, implying a zero group velocity.In contrast, the spatially periodic quilting in the bursting coherence map arises from periodic traveling waves whose group velocity can be deduced from the inverse slope of the k-vs-f spectral trend, giving a value of around 7:0 cm=s, consistent with the slope of the x-vs-t zigzagged wave-front chevrons in Fig. 8(b) (note that the group velocity is much slower than the axonal conduction speed of 140 cm/s: see Table I).
For the anesthetized cortex, we surveyed the full range of diffusion-controlled dynamics from zero diffusion (all gap junctions closed: D 2 ¼ 0) to strong inhibitory coupling (gap junctions open: D 2 $ 1:0 cm 2 ).For each 20-s experiment, we computed the spatially averaged global coherence for the final 2 s of recording; the results appear in Fig. 7(b).Although the wake [Fig.7(a)] and coma [Fig.7(b)] global coherence graphs are broadly similar, the effect of anesthetic is to completely eliminate seizures at the lowest diffusion values (''suppression'': 0 & D 2 =cm 2 & 0:16) and to shift the activated ''noncognitive-wake'' coherence peak of Fig. 7(a) to the right.This latter change is consistent with a comatose cortex requiring stronger Turing instability to support an activated state.
Lying between the Hopf-dominated coherence peak [on the left of Fig. 7(b)] and the Turing-dominated peak (on the right) is a broad intermediate zone (0:4 & D 2 =cm 2 & 0:7) of reduced coherence whose oscillations are chaotic: large amplitude, low frequency, but very irregular-in marked contrast to the regular periodic traveling waves for anesthetic bursting.This contrast, visualized earlier in the space-time strip charts of Fig. 8(b) (rows 4 and 5), is reinforced in the full 20-s time series shown in Fig. 10 and quantified in the spatially averaged power spectra of Fig. 11.
We identify the mixed-mode chaotic dynamics as ''anesthetic slow-wave'' because its time series and spectrum are very similar to the distorted slow-wave sleep patterns seen in general anesthesia (e.g., compare the delta waves of natural and anesthetic slow-wave sleep illustrated in Fig. 1A, B of [50] and in Fig. 2A, B of [51]).The modelgenerated anesthetic slow-wave spectrum of Fig. 11 shows a very broad energy distribution over the low-frequency domain 0:05 & f & 2:0 Hz, whereas for the seizure state, spectral energy is bunched to show strong resonant peaks at around 1.75 Hz and its harmonics.Both spectra exhibit a steep power-law decay at higher frequencies.

C. Effect of two-point connectivity on wake and coma dynamics
We investigate the impact of varying both orientation (horizontal, vertical, diagonal) and separation Ár¼jr 1 À r2 j of the two directly connected grid points, but generally found no significant differences in the fully developed cortical dynamics across the range of gap-junction diffusion settings (variations in D 2 ) in either cortical state (wake vs coma).Although the connected grid points act as a pair of excitatory seed locations that hasten symmetry breaking of the initial homogeneous equilibrium state, the resulting outward-propagating target patterns are transients that dissipate rapidly.
But this is not the case for seizure and anesthetic bursting: For both wake and coma experiments (row 5 of Figs. 5 and 8: D 2 ¼ 0:3 cm 2 ), the heterogeneous source points act as foci to ''pin'' outgoing seizure wave fronts.In a clinical context, having identified a focus of epileptic seizure events in cortical tissue, a clinician might plan to surgically excise or lesion the focus in order to suppress the seizures.We can perform a virtual excision of the seizure foci by setting ¼ 0 to eliminate the heterogeneity, then rerunning our simulations.The outcome can be deduced by comparing the orange ( ¼ 200) and black ( ¼ 0) global coherence curves of Fig. 12: For both wake (upper panel) and coma (lower panel), removal of the heterogeneity eliminates the seizure-bursting state at D 2 ¼ 0:3 cm 2 .
The surprising feature of Fig. 12 is that, apart from the Hopf-dominated seizure end of the curve, the presence of two-point heterogeneity makes little or no difference to the global coherence statistics.It appears that, provided there is sufficient inhibitory diffusion (D 2 * 0:4 cm 2 ), the Turing instability is able to activate normal cortical dynamics, and, at the global scale, the small local excitation from two-point connections (evident as parallel contrails of activity at x 1;2 ¼ 20, 60 in Figs. 5 and 8) is neither needed nor noticed.Thus, adding a point source of activation to an 11.Spatially averaged j Qe ðfÞj 2 power spectra for anesthetic slow-wave (thick-black curve: D 2 ¼ 0:4 cm 2 ) and bursting (gray curve: D 2 ¼ 0:3 cm 2 ) time series illustrated in Figs.10(a) and 10(b), respectively.Power is expressed in dB relative to the value at 100 Hz.Bursting, but not slow-wave, shows pronounced resonances at f 1 ¼ 1:75 Hz and its harmonics nf 1 , n ¼ 2 . . .12. The red line shows the best-fit trend for the power law P $ f À4:4 .Frequency resolution for both spectra is 0.05 Hz.  already-activated cortex has little effect on cortical dynamics, though the presence of a point source can seed the onset and precipitation of that dynamics from an initially homogeneous equilibrium state.

IV. DISCUSSION
In this paper, we have investigated the symmetrybreaking impact of bringing direct long-range connectivity into our previously isotropic mean-field cortical model.We have deliberately selected the simplest case of a single pair of two-way axonal connections so that we can identify the essential changes to cortical dynamics wrought by imposition of spatial heterogeneity; generalization to anatomically realistic maps of regional interconnectivity (e.g., using the CoCoMac primate connectivity database [52]) will form the basis of future work.We selected two reference states for detailed dynamical analysis: a noncognitive-wake activated state lying on the high-firing upper branch of the Fig. 2 equilibrium manifold, and an anesthetic-coma state on the quiescent lower branch.Both states lie in the vicinity of a c2TH point that allows dynamic interaction between spatial (Turing) and temporal (Hopf) instabilities.We were able to alter the relative balance between these competing instabilities by tuning D 2 , the diffusive-coupling strength for the gap-junction connections between adjoining inhibitory neurons: Strong diffusion encourages formation of spatial patterns (Turing structures); weak diffusion favors whole-of-cortex Hopf oscillations.Intermediate values of diffusion generate turbulent mixed states exhibiting low-frequency spatiotemporal chaos.
It should be noted that for each of our numerical simulations, we have treated the inhibitory diffusion strength D 2 as a fixed control parameter.In reality, gap-junction conductivity is not fixed, but is dynamically modulated by a plethora of physiological and pathological stimuli.For example, states of conscious attention are associated with increased brain concentrations of neuromodulator amines which tend to open gap junctions; conversely, states of somnolence are associated with increased concentration of the inhibitory neurotransmitter GABA which is known to close gap junctions [53,54].We have not attempted to model these dynamic feedback effects here.

A. Impact of two-point connections on cortical dynamics
Provided the level of Turing activation is moderately high (e.g., D 2 * 0:4 cm 2 ), we find that the inclusion of two-point connectivity makes little difference to either noncognitive-wake or ''anesthetic-coma'' dynamics.This apparent insensitivity to the presence of two-point connections arises because the cortex is already active, so the incremental contribution to cortical excitation from these point sources is swamped by ambient network activity.However, this is not the case at lower levels of cortical activation (e.g., D 2 & 0:3 cm 2 ).Without the point-source connections, an anesthetized cortex would remain in a low-firing suppressed state [see Fig. 12(b)], but with direct connections enabled, the incremental excitation is sufficient to flip the cortical state from suppression to seizurelike bursting.And for the noncognitive-wake case, enabling direct connections can switch the dynamics from a disordered, chaotic state to an ordered, coherent seizure state [Fig.12(a)].
Our cortical model pictures the cortex as a continuum of cortical tissue containing a single two-point axonal connection.In contrast, other workers have built network models of the cortex consisting of neural hubs that are linked via an anatomically informed point-to-point connectivity matrix [33,34,55,56] derived from diffusiontensor imaging of the human brain or fiber tracing of the primate cerebral cortex.Despite this obvious difference in modeling philosophies-a single fiber embedded in a cortical continuum in the present work versus networks of neural hubs in those other works-we see some unexpected and interesting similarities in model behaviors when we compare the dynamics of our lower (quiescent) branch with anatomically motivated network models for the brain rest state.Specifically, our finding that the lowactivity quiescent state of the brain is most sensitive to large-scale connectivity aligns nicely with the notion of criticality proposed by Ghosh et al. [33], later refined by Deco and colleagues [55,57]: namely, that, during rest, the brain operates at the edge of instability in order to maximize noise-evoked explorations of its state space.Deco and Jirsa [55] reported that their network model gave the best agreement with observed fMRI functional connectivity when the model was close to the transition region separating low-and high-firing states, and that the resulting dynamics was influenced by the presence of latent ''ghost'' attractors.In our work, spatiotemporal patterns are most likely to evolve from the lower equilibrium branch when the model is close to criticality, and these can arise even when there is no stable upper branch (e.g., at the coma state of Fig. 3), implying ghost attraction to the latent upper state.If two-point connections are enabled, a larger region of the resting branch is able to support pattern formation [compare black and orange curves of Fig. 12(b)], consistent with the idea that large-scale connections can provide anatomical guidance for rest-state dynamics.In future work, the single long-range connection will need to be replaced by an anatomically informed hierarchical connectivity map.Kaiser et al. [58] have examined the spreading of activation in neural networks and shown that hierarchical cluster networks-such as those found in the mammalian cortex-are more resistant to runaway network activation (''epileptic seizure'') than are random networks or small-world networks that lack hierarchical modules.This suggests that imposition of a modular topology onto our cortical model might tend to limit the spreading of spatiotemporal patterns to local modules unless there is strong coupling between neural modules.

B. Effect of inhibitory diffusion on noncognitive-wake dynamics
Our noncognitive-wake state at D 2 % 0:7 cm 2 is characterized by small-scale coherent oscillations that modulate a labyrinthine Turing pattern of elevated (up-state) and suppressed (down-state) activity (top row of Fig. 5).Boosting the diffusion strength to D 2 ¼ 0:8 cm 2 strengthens the spatial patterning while simultaneously extinguishing the Hopf oscillations and their phase coherence [Fig.7(a)], resulting in a low-coherence ''frozen Turing'' structure.Under normal circumstances, this frozen state is probably pathological, but may possibly be relevant to early brain development when the immature brain is richly endowed with gap junctions [38] and the conditions for formation of permanent Turing patterns are most favorable.Cartwright [59] has proposed that the morphogenesis of the mazelike features of the convoluted ridges (gyri) and valleys (sulci) of the cortex is a Turing solution, but his speculation assumes an axonal guidance competition between (unknown) inhibitor and activator chemical species rather than a gap-junction process in which inhibitory diffusion D 2 dominates excitatory diffusion D 1 .
Lowering the diffusion strength to D 2 ¼ 0:5 cm 2 strengthens the Hopf instability while destabilizing Turing structuring, resulting in the spontaneous emergence of a turbulent dynamics that is chaotic in space and time (third row of Fig. 5).Further reduction in D 2 eventually leads to a highly coherent seizurelike state.Thus, the model predicts that the path from coherent noncognitivewake to coherent seizure should cross an intermediate chaotic regime of low coherence.This prediction is in agreement with the clinical report of Mormann et al. [60] who observed decreased synchronization in EEG recordings during the period leading up to seizure.

C. How do anesthetic agents affect cortical coherence?
As well as boosting the strength of the IPSP [43,44], most general anesthetics also inhibit gap-junction communication [61].Consequently, the transition from awake to anesthetic coma on Fig. 3 should correspond to a movement from the noncognitive-wake high-coherence peak of Fig. 12(a) to a region of lower diffusion (and diminished coherence) labeled ''anesthetic slow-wave'' on the Fig. 12(b) graph.Slow-wave spatiotemporal dynamics, pictured in Fig. 8 (rows 2-5) and Fig. 10(a), is characterized by low-frequency, chaotic transitions between low-firing down-states and activated upstates.The notion of chaotic-state transitions is consistent with clinical observation of changes in phase coupling between pairs of EEG channels during propofol-induced anesthesia: For most channel pairs, Koskinen et al. [36] observed a systematic loss in phase synchrony in the subdelta band (0.05-1.0 Hz) during induction, followed by an increase during recovery of consciousness.
We need to acknowledge an important limitation to our coherence predictions.Because our model is devoid of subcortical structures such as the thalamus, it cannot inform on the anesthetic effect on thalamocortically generated rhythms such as alpha-band (8-12 Hz) and spindle (12-14 Hz) oscillations.Recent papers (e.g., [62][63][64]) show that transition into unconsciousness induced by propofol is accompanied by a decrease in alpha-band coherence in posterior regions, concomitant with an increase in coherence in frontal areas.This boost in alpha-band coherence is almost certainly the result of anesthetic-induced thalamic hyperpolarization; inclusion of thalamocortical feedback will be the subject of future work.

D. Chaotic slow oscillations in nonREM sleep?
The slow oscillation is the defining feature of scalpmeasured EEG under general anaesthesia and also during the slow-wave stage of natural sleep.Preliminary model investigations (not shown here) suggest that a spatiotemporally chaotic slow-wave state can emerge in the nonREM sleep state, provided that gap-junction diffusivity is reduced sufficiently.This is plausible, since the GABA neurotransmitter is abundant in nonREM sleep [65] and is thought to close gap junctions [53].Furthermore, the so-called ''wake-up'' pills (e.g., modafinil) are known to open gap junctions [49].Our suggestion that nonREM sleep might be chaotic is concordant with clinical observations that the cortex is less functionally connected during slow-wave sleep [66], and is consistent with a recent study reporting that anesthesia-induced slow oscillations occur asynchronously across the cortex, fragmenting global network activity while maintaining local (< 4 mm) connectivity patterns [67].In our model, turbulent wave fronts generated by chaotic interactions between Turing and Hopf instabilities support localized patterns of correlated activity that is disrupted at larger length scales [Fig.10(a)].

E. Implications of chaotic oscillations for memory processing during slow-wave sleep
Tononi and Cirelli [68] argue that slow-wave sleep is essential to restore synaptic weights-that were elevated during prior wakeful learning-to an energetically sustainable baseline level; this is the synaptic homeostasis hypothesis.Recent studies by Liu et al. [69] provide direct experimental evidence for both wake-related increases and sleep-related decreases, in cortical synaptic strength.We speculate that the turbulent dynamics of slow-wave spatiotemporal chaos might provide a natural substrate for synaptic down-scaling during sleep.Suppose that, at a particular instant, the sleeping cortex is in a well-mixed turbulent state such as that illustrated in row 4 of Fig. 8(c).About one-third of the cortex is in the excited (red) upstate, and these activated limbs are surrounded by larger regions of the low-firing (blue) down-state neural populations.These inactive neighbors comprise the bulk of the instantaneous presynaptic input to the up-state populations, and similarly, the majority of the up-state connections will be presynaptic to the down-state populations.Therefore, most synaptic activity is between populations in opposed states.Assuming a population-scale version of Hebbian ''unlearning'' applies here (''neurons out of sync fail to link''), then we might expect the unlearning to take the form of a gradual reduction in the synaptic weights at the interface between the opposed populations.At each successive instant, the turbulent mixing creates new up-down partitions of the cortex; thus, the unlearning rule is able to propagate throughout in a random and unbiased fashion, carried on the chaotic wave fronts.In addition, the lack of fixed timing relationships in a turbulent network could prevent the induction of spike-timing-dependent plasticity.

F. Mechanisms for the slow oscillation
What is the underlying mechanism for the slow cortical delta-wave oscillations that characterize deep sleep and general anesthesia?Physiological measurements show slow cyclic changes in ionic concentrations of extracellular Ca 2þ and K þ that are correlated with the transitions between the up and down phases of the slow oscillation [70][71][72], but it is presently unknown whether these ionic changes cause the oscillation-via a slow homeostasis that increases excitation in the down-state and suppresses excitation in the up-state-or are caused by the oscillation.Typically, cortical modelers have attempted to capture the slow ionic dynamics by incorporating an additional equation (with a suitably large time constant) to modulate, for example, the sigmoidal firing-rate function [73] or the excitatory synaptic strength [74].
Our approach is quite different.In our model, the slow oscillation is generated spontaneously by the same Hopf temporal instability that-in its uncontrolled pathological phase-has been linked to whole-of-cortex seizure [24][25][26][75][76][77].But, as shown here, provided there is sufficiently strong inhibitory gap-junction interneuronal coupling, a Turing spatial instability can emerge to interact with-and moderate-the Hopf oscillations, breaking the global coherence patterns of seizure by creating turbulent structures and chaotic slow oscillations.The notion that seizure and the slow oscillation are intimately linked is consistent with the experimental observations by Steriade et al. that spontaneous seizures at 2-4 Hz can develop without discontinuity from the slow (0.5-0.9 Hz) sleeplike oscillations of general anesthesia [78].
We have shown that-depending on neuromodulatory conditions-the mean-field cortex provides a substrate that supports either spatial patterns or global oscillations; we postulate that normal brain function relies on a dynamic balance between these two competing instabilities.As seen in Figs.7 and 12, inhibitory diffusion strength D 2 provides a sensitive control parameter that alters cortical dynamics from Hopf-dominated seizure for weak diffusion to Turing-dominated ultraslow beating patterns for strong diffusion settings.We find that incorporating a nonlocal direct connection between two points in the cortical grid makes no appreciable difference to dynamics or coherence patterns when diffusion is strong and the cortex is well activated, but does serve as an additional source of excitation in the low-diffusion state, tending to increase propensity to seizure.We do acknowledge, however, that this idealized nonlocal connection cannot represent real anatomical structure, a shortcoming that will be addressed in future work.

FIG. 2 .
FIG.2.Manifold of steady-state firing rates Q e across the ð; ÁV rest Þ anesthesia domain.The control parameter sets the anesthetic effect; ÁV rest is an additive offset representing background cortical excitation (ÁV rest > 0) or suppression (ÁV rest < 0).The yellow curve marks the edge of the reentrant ''fold''; the dashed black curve shows the projection of this edge onto the lower and upper surfaces, bounding the zone of multiple steady states.The red-green-blue curve shows the distribution of steady states for varying anesthetic inhibition at constant cortical excitation (see Fig.3).

FIG. 5 .
FIG. 5. Effect of inhibitory diffusion on wake state.Grid simulations, initialized at the high-firing wake state of Fig. 3, show dynamical effect of stepped reductions in diffusion strength from D 2 ¼ 0:7 cm 2 (top row) to D 2 ¼ 0 (bottom).(a) Time-series extracts showing excitatory firing rates Q e ðtÞ for the final 4 s of a 20-s run at x ¼ 30 ( black lines), x ¼ 85 ( gray lines) on the midline (y ¼ 60) of the 120 Â 120 cortical grid.(b) Q e ðt; xÞ space-time strip charts showing x-axis activity along the y ¼ 60 midline strip for a full 20-s simulation.(c) Q e ðy; xÞ bird's-eye view of cortical activity at t ¼ 20 s; the color bar indicates the firing rate in spikes/s.Simulation settings: toroidal boundary conditions; integration time step Át ¼ 0:4 ms.Stimulus: spatiotemporal white noise with scale factor ¼ 4; direct two-way connection joins x 1 ¼ 20 to x 2 ¼ 60 on the y ¼ 60 midline.

FIG. 6 .
FIG. 6. Fourier spectra and phase coherence maps for the wake strip charts of Fig. 5. (a) Qe ðf; kÞ shows the 2D discrete Fourier transform of Q e ðt; xÞ images of Fig. 5(b); f is the temporal frequency (Hz), and k is the spatial frequency in cycles=L x , where L x ¼ 25 cm is the x-axis side length of the cortical grid.Inset line graphs in black compare the average temporal spectrum Qe ðfÞ (horizontal format: summed over wave number k) against the average 1D spatial spectrum Qe ðkÞ (vertical format: summed over frequency f).(b) Phase coherence maps Rðx; x 0 Þ showing the degree of firing-rate coherence for all time-series pairs Q e ðt; xÞ, Q e ðt; x 0 Þ along the y ¼ 60 midline of a 120 Â 120 cortical grid; coherence maps are computed from Hilbert transforms of the final 2 s (from t ¼ 18 to 20 s) of Fig. 5(b).The color bar shows the mapping between color and phase coherence.
of Ref.[42] for sample divergence trends) and demonstrated early exponential divergence of adjacent state-space trajectories; thus, these low-coherence intermediate states (for 0:4 & D 2 =cm 2 & 0:6) have a positive dominant Lyapunov exponent and chaotic dynamics.The simulation results presented in Figs. 5 and 6 are six representative samples drawn from over 2000 numerical experiments that surveyed the full range of cortical dynamics from pure Hopf oscillation (no diffusion: D 2 ¼ 0) to frozen Turing patterns (strong diffusion: D 2 ¼ 1:0 cm 2 ).For each experiment, we extracted a scalar measure of the global (i.e., spatially averaged) phase coherence across the cortical grid by computing the mean value of its Rðx 0 ; xÞ [Fig.6(b)] phase coherence map; the survey results are presented in Fig. 7(a).Although some regions of the graph show considerable scatter, the underlying trend is clear: high global coherence in the wake state (D 2 =cm 2 % 0:7); diminished coherence in the intermediate mixed states (0:4 & D 2 =cm 2 & 0:6); and back to even higher coherence values in the seizure state (D 2 =cm 2 & 0:3).[The ''frozen Turing'' states (D 2 =cm 2 * 0:8) have the lowest global coherence values; however, with permanently elevated and suppressed limbs of unmodulated cortical activity, these states are likely to be pathological, but may have significance in cortical morphogenesis; this idea is briefly developed in Sec.IV.] Within the seizure regime [top-left corner of Fig. 7(a): D 2 =cm 2 & 0:2],we see a splitting in coherence values, suggesting that, when diffusion is weak, the seizing cortex has access to multiple and distinct stable dynamic states, leading to the possibility of hysteresis effects.An examination of the detailed bifurcation structure and neurophysiological implications will be reported in future work.

VIDEO 1 .FIG. 7 .
FIG. 7. Phase coherence trends for two-point-connected cortex.Spatially averaged phase coherence measurements for (a) high-firing wake ( i ¼ 1:00) and (b) low-firing coma ( i ¼ 1:018) states of the noise-driven cortex for default two-point connectivity strength ( ¼ 200) are shown.Gray circles represent the outcomes of 20-s numerical experiments at a given of inhibitory in the 0.0 to cm 2 in steps of 0.01 cm; simulations were repeated 20 times at each D 2 value for a total of 2020 data points per panel.For each experiment, we computed the Rðx 0 ; xÞ coherence matrix for final 2 s of the space-time record [e.g., see Figs. 5(b) and 6(b)], then extracted its upper-triangular matrix mean as an estimate of global phase coherence.Bold-black curves show the data median to guide the eye.Solid-(striped-) red bars indicate regions of a strongly (weakly) positive Lyapunov exponent as extracted from paired noise-free simulation runs.

FIG. 8 .
FIG. 8. Effect of inhibitory diffusion on the anesthetic-coma state.Grid simulations are given, initialized at the low-firing coma state of Fig. 3, showing the effect of stepped reductions in inhibitory diffusion strength from D 2 ¼ 0:8 (top) to D 2 ¼ 0:1 cm 2 (bottom).(a) Time-series extracts showing excitatory firing rates Q e ðtÞ for the final 4 s of a 20-s run at x ¼ 30 (black lines) and x ¼ 105 (gray lines) on the y ¼ 60 midline of the 120 Â 120 cortical grid.(b) Q e ðt; xÞ space-time charts showing x-axis activity along the y ¼ 60 midline for a full 20-s simulation.(c) Q e ðy; xÞ bird's-eye image of cortical activity at t ¼ 20 s; the color bar indicates the firing rate in spikes/s.Points x 1 ¼ 20 and x 2 ¼ 60 on the y ¼ 60 midline are directly connected; see Fig. 5 for the remaining simulation details.

FIG. 9 .
FIG. 9. Fourier spectra and phase coherence maps for the coma strip charts of Fig. 8. (a) Qe ðf; kÞ shows the 2D discrete Fourier transform of Q e ðt; xÞ images of Fig. 8(b); f is the temporal frequency (Hz), and k is the spatial frequency (cycles=L x ).Inset line graphs in black compare the average temporal spectrum Qe ðfÞ (horizontal format) against the average 1D spatial spectrum Qe ðkÞ (vertical format).(b) Phase coherence maps Rðx; x 0 Þ showing the degree of firing-rate coherence for timeseries pairs Q e ðt; xÞ; Q e ðt; x 0 Þ (with fx; x 0 g 2 f1; 2; . . .120g) along the y ¼ 60 midline; coherence maps are computed from the final 2 s of Fig. 8(b) strip charts.The color bar shows the mapping between color and phase coherence.
FIG.11.Spatially averaged j Qe ðfÞj 2 power spectra for anesthetic slow-wave (thick-black curve: D 2 ¼ 0:4 cm 2 ) and bursting (gray curve: D 2 ¼ 0:3 cm 2 ) time series illustrated in Figs.10(a) and 10(b), respectively.Power is expressed in dB relative to the value at 100 Hz.Bursting, but not slow-wave, shows pronounced resonances at f 1 ¼ 1:75 Hz and its harmonics nf 1 , n ¼ 2 . . .12. The red line shows the best-fit trend for the power law P $ f À4:4 .Frequency resolution for both spectra is 0.05 Hz.

FIG. 12 .
FIG.12.Phase coherence trends for (a) wake i ¼ 1:00) (b) coma ( i ¼ 1:018) for a homogeneous cortex devoid of two-point connections (i.e., ¼ 0).Bold-black curves show the median trend through data points (gray circles) obtained from cortex experiments; these ¼ 0 boldblack trends are superimposed on the ¼ thin-orange median fits from 7(a) and 7(b) to illustrate the impact of suppressing the two-point axonal connection.

TABLE I .
Symbol definitions and standard values for the cortical model.