Spinodal decomposition in a mean-field model of the cortex: Emergence of hexagonally symmetric activation patterns

Spinodal decomposition is a well-known pattern-forming mechanism in metallurgic alloys, semiconductor crystals, and colloidal gels. In metallurgy, if a heated sample of a homogeneous Zn-Al alloy is suddenly quenched below a critical temperature, then the sample can spontaneously precipitate into inhomogenous textures of Zn- and Al-rich regions with significantly altered material properties such as ductility and hardness. Here we report on our recent discovery that a two-dimensional model of the human cortex with inhibitory diffusion can, under particular homogeneous initial conditions, exhibit a form of nonconserved spinodal decomposition in which regions of the cortex self-organize into hexagonally distributed binary patches of activity and inactivity. Fine-scale patterns precipitate rapidly, and then the dynamics slows to render coarser-scale shapes which can ripen into a range of slowly evolving patterns including mazelike labyrinths, hexagonal islands and continents, nucleating “mitotic cells” which grow to a critical size then subdivide, and inverse nucleations in which quiescent islands are surrounded by a sea of activity. One interesting class of activity coalesces into a soliton-like narrow ribbon of depolarization that traverses the cortex at $\sim 4\mathrm{cm}/\mathrm{s}$. We speculate that this may correspond to the thus far unexplained interictal waves of cortical activation that precede grand-mal seizure in an epileptic event. We note that spinodal decomposition is quite distinct from the Turing mechanism for symmetry breaking in cortex investigated in earlier work by the authors [Steyn-Ross et al., Phys. Rev. E 76, 011916 (2007)].

Figure 1

Manifold of cortical equilibrium states. (A) Bird's-eye view showing the three-state regions in copper (light gray) and single-state surroundings (black) as a function of inhibitory gain ${\lambda }_i$ and excitatory offset $\mathrm{\Delta }{V}_e^{\text{rest}}$. The coexistence region is bounded by the “binodal coastlines” in (A) and by the solid-yellow (light gray) and dot-dashed black curves in (B). The red (medium-gray) circle marks critical point CP. The three ($\oplus$) markers labeled (a), (b), and (c) locate the $(\mathrm{\Delta }{V}_e^{\text{rest}},{\lambda }_i)$ coordinate resulting in the spinodal decompositions illustrated in the meander, honeycomb, soliton, and nucleation panels. (B) Three-dimensional view of distribution of excitatory firing rates $Q_e$ across the ${\lambda }_i\text{–}\mathrm{\Delta }{V}_e^{\text{rest}}$ domain. The (a), (b), and (c) simulations all lie within the binodal multiroot region characterized by a stable dark-red (midgray) upper “active phase” and deep-blue (dark-gray) lower “quiescent” branch, separated by an unstable midbranch. All simulations are initialized to commence at the unstable midbranch equilibrium point.

Figure 2

Linear stability predictions as a function of wave number for the four decomposition simulations indicated in Fig. 1 and captured in the (a) meander (Fig. 5), (b) nucleation (Fig. 6), (c) honeycomb (Fig. 4), and (d) soliton (Fig. 10) snap-shot galleries. From top to bottom, panel rows display, respectively, the dispersion trends for the top, middle, and bottom homogeneous equilibrium states. Dispersion graphs have been ordered from left to right to show progressive increases in the degree of instability of the middle branch. Solid-blue (black) curves show growth rate $\alpha$ (in ${\mathrm{s}}^{-1}$) predicted from linear eigenvalue analysis, with $\alpha \left(q\right)>0$ indicating instability at wave number $q$; dashed-red (dashed-gray) curves show corresponding $\omega \left(q\right)/2\pi$ (Hz) eigenfrequency trends. In all cases, top- and bottom-branch equilibria are stable across all wave numbers, although the top branch for (d) soliton is very close to a 3.5-Hz Hopf bifurcation ($\alpha =-0.001{\mathrm{s}}^{-1}$ at zero wave number). The emergent Turing-like labyrinth structure in (a) seems to be guided by the presence of a damped Turing bifurcation $(\oplus )$ on the bottom branch near $q/2\pi \approx 0.5{\mathrm{cm}}^{-1}$.

Figure 3

Hypothetical free energy curves and phase diagram for an idealized “temperature”-driven spinodal decomposition of cortical activity. Here $Q$ represents the firing rate of the excitatory neural population with $0\le Q\le Q_e^{\text{max}}$. $T_c$ is the critical “temperature” below which homogeneous cortical activity $Q_{\text{mid}}$ spontaneously separates into high- and low-activity states $Q_{\text{hi}}$ and $Q_{\text{low}}$, respectively. (a) Black upper free energy curve for elevated temperature $T_2>T_c$ shows a single minimum at activity $Q_{\text{mid}}$. Curve deforms downward into the blue (dark-gray) double-well potential at reduced temperature $T_1<T_c$ with minima corresponding to separated phases (binodal points) $Q_{\text{hi}}$ and $Q_{\text{low}}$. Points of inflexion $Q_{-}$, $Q_{+}$ mark the bounds of the unstable spinodal region. (b) Activity–temperature phase diagram showing a rapid temperature quench (vertical flow) from $T_2$ to $T_1$ inside the spinodal region, leading to separation into distinct phases (horizontal flows). (Panel B modified from Fig. 1 of Ref. [20].)

Figure 4

Spinodal decomposition into a “honeycomb” structure: red (light gray with white borders) hexagonal islands of activity surrounded by narrow blue (dark-gray) canals of inactivity. Precipitation of jagged bulk structures occurs rather rapidly ($\sim 1\mathrm{s}$), and then the fine coastal details disappear as the pattern ripens on a much slower timescale (approximately tens of seconds). Parameter settings: $\left[\mathrm{\Delta }{V}_e^{\text{rest}}/\text{mV},{\lambda }_i,D_2/{\text{cm}}^2,{\gamma }_i^0/{\mathrm{s}}^{-1}\right]=\left[-1.85,0.7843,0.3,80\right]$; homogeneous model cortex was started on the midbranch equilibrium point for location (b) in Fig. 1. See Video 1 in the Supplemental Material [35] for a 60-s movie to accompany honeycomb structure emerging from homogeneous midbranch equilibrium.

Figure 5

Spinodal decomposition into a “meander” structure: a mazelike labyrinth of continuously linked blue (black) canals of quiescence isolating red (white-bordered gray) islets of activity. As for the honeycomb of Fig. 4, the bulk pattern is already apparent after $\sim 1\mathrm{s}$ and then matures more slowly. Settings: $\left[\mathrm{\Delta }{V}_e^{\text{rest}}/\text{mV},{\lambda }_i,D_2/{\text{cm}}^2,{\gamma }_i^0/{\mathrm{s}}^{-1}\right]=\left[1.3,1.0,0.35,80\right]$; initial state: midbranch steady state for location-(a) in Fig. 1. See Video 2 in the Supplemental Material [35] for a 100-s movie to accompany emergence of a Turing-like labyrinthine cortical activation.

Figure 6

Spinodal decomposition into an almost regular “nucleation” pattern: a cellular array of red (white-bordered gray) activity nuclei in a sea of blue (black) quiescence. Larger structures pinch off and separate into daughter cells; smaller structures grow until they become unstable, whereupon they elongate into a dumbbell shape that thins and divides into pairs of daughter cells in a process reminiscent of biological mitosis. Cell population multiplies until the evident repulsive force between nearest neighbors suppresses the tendency of each cell to grow and divide. Settings: $\left[\mathrm{\Delta }{V}_e^{\text{rest}}/\text{mV},{\lambda }_i,D_2/{\text{cm}}^2,{\gamma }_i^0/\text{s}{}^{-1}\right]=\left[-2.5,0.8,0.45,45\right]$; initial state: midbranch steady state for location-(c) in Fig. 1. See Video 3 in the Supplemental Material [35] for a 60-s movie to accompany emergence of an array of activity cells undergoing mitosis.

Figure 7

Impact of noise intensity on emergence of nucleation spinodal patterns. Each column corresponds to one of four simulation experiments $k=[{10}^{-4},1,5,80]$, where $k$ scales the amplitude of the default setting for white-noise coefficient $s_0=0.125$ [see Eq. (4) and Table 1]; time increases vertically downward. Each simulation uses the identical sequence of random numbers. In general, stronger stochastic stimulation promotes faster nucleation, but patterns became unstable if the noise is too intense (fourth column: $k=80$).

Figure 8

Time development of the cellular nucleation pattern of Fig. 6 expressed as relative proportions of spectral content above and below a reference wave number, $q/2\pi =0.38{\mathrm{cm}}^{-1}$. At $t=0$, the unstable midbranch of the homogeneous model cortex is perturbed with small-amplitude spatial white noise. The upper red (gray) curve shows that the initially high spatial-frequency fraction falls precipitously as large amorphous structures emerge spontaneously to break symmetry, reaching a minimum at $t_1\approx 0.4$ s. The lower blue (black) curve tracks the low-frequency spatial energy which peaks at this point, signaling the onset of a structural phase transition. The high-wave-number fraction then grows rapidly as the differentiation into high- and low-firing cortical patches proceeds and then matures more slowly after $t_2\approx 1.25$ s.

Figure 9

Hexagonal structuring of (a) cellular and (d) labyrinthine decomposition patterns revealed by 2D spatial Fourier transforms in panels (b) and (c), respectively; panels (c) and (f) show the amplitude spectra obtained after averaging over 25 sets of 100-s grid simulations. The radius of the primary ring structures ($q/2\pi =$ 0.44 and 0.40 ${\mathrm{cm}}^{-1}$) correspond to the inverse of the characteristic repetition lengths of the cell and maze patterns of (a) and (d). In (a), some representative cells have been joined (white lines) to aid the eye; while most show the expected six-nearest-neighbor hexagonal links, structure defects are evident with some seven- and five-neighbor arrangements.

Figure 10

Destabilization of spinodal phase separation by close proximity to a Hopf bifurcation. Cortex was initialized at the same $(\mathrm{\Delta }{V}_e^{\text{rest}},{\lambda }_i)$ starting coordinate that produced the stationary honeycomb pattern of Fig. 4, but the inhibitory rate-constant was reduced by a factor of $\sim 4$ to induce a marginally damped Hopf instability in the upper branch (see Fig. 2). Although initial coastline shapes are similar (compare the 0.2-s panel here with the 0.3-s panel of Fig. 4), temporal evolution is now qualitatively different with vigorous competition between red (white-bordered gray) activated and blue (black) quiescent patches that eventually resolves into a soliton-like wave of activity that traverses the cortex at $\sim 4\mathrm{cm}/\mathrm{s}$. See Video 4 in the Supplemental Material [35] for a 38-s movie to accompany emergence of a soliton shockwave in the cortical model. Settings: $\left[\mathrm{\Delta }{V}_e^{\text{rest}}/\text{mV},{\lambda }_i,D_2/{\text{cm}}^2,{\gamma }_i^0/\text{s}{}^{-1}\right]=\left[-1.85,0.7843,0.4,22\right]$.

Figure 11

Time series for a seizurelike precursor event in mouse cortex. The electrical activity of a 400-$\mu \mathrm{m}$-thick coronal slice of mouse brain tissue is sampled with five numbered Ag/AgCl electrodes (left photo). A seizurelike event (SLE) impinges simultaneously on electrodes 1 and 2 (right-hand traces) and then traverses the cortical rind to reach electrodes 3, 4, and 5 at successively later times. The mesh spacing is $\sim 1.5\mathrm{mm}$, giving an SLE travel speed of 3 to 5 cm/s. [See the Appendix for experiment details.]

Figure 12

Gallery of representative spinodal decomposition patterns for 18 $(\mathrm{\Delta }{V}_e^{\text{rest}},{\lambda }_i)$ coordinates located slightly inland from the lower binodal coastline. Each snapshot shows early pattern development 10 s after launching from the unstable midbranch separatrix. In addition to precursor honeycomb (panels 1–3), red-on-blue (gray-on-black) nucleation (4–10), and meander (11–16) patterns, blue-on-red (black-on-gray) inverse nucleations are also evident (panels 17 and 18). Spatiotemporal settings: $\left[D_2/{\text{cm}}^2,{\gamma }_i^0/\text{s}{}^{-1}\right]=\left[0.35,80\right]$.

Figure 13

Spinodal patterns (first row) and spatial spectra (second row) developed after 10 s on a $6\times 6\text{−}\mathrm{cm}$ cortex. From left to right, respective inhibitory diffusion ($D_2$) and excitatory axonal decay rate (${\mathrm{\Lambda }}_{eb}$) are $D_2/{\text{cm}}^2=[0.20,0.05,0.02,0.015]$ and ${\mathrm{\Lambda }}_{eb}/{\text{cm}}^{-1}=[24,24,24,28]$. Other model settings are identical to those used for the nucleation patterns of Figs. 1 and 6 on a $25\times 25\text{−}\mathrm{cm}$ cortex. [(a)–(d)] As diffusive and axonal connectivities reduce, the characteristic length scales of the spinodal patterns shrink, and, consequently, the dominant spatial mode (second row) increases in frequency (units: ${\mathrm{cm}}^{-1}$). White ring indicates dominant spatial frequency obtained from radial averaging [see Fig. 14].

Figure 14

(a) Radially averaged spatial activity for spinodal pattern Fig. 13. (b) Wave-number trend with axonal decay-rate and inhibitory diffusion (right) across 20 spinodal simulations and two cortical domain sizes. Evidently the dominant wave number follows a linear trend $q_{\text{dom}}\propto {[{\mathrm{\Lambda }}_{eb}/\sqrt{D_2}]}^{1/2}$ across a sixfold change in characteristic length scale from $1/0.44{\text{cm}}^{-1}=2.3$ cm (lower-left datum) to $1/2.67{\text{cm}}^{-1}=0.37$ cm (upper-right datum). Error bars show $\pm 2/\mathrm{\Delta }x$ (blue circles: $25\times 25\text{−}\mathrm{cm}$ cortex) and $\pm 1/\mathrm{\Delta }x$ (red crosses: $6\times 6\text{−}\mathrm{cm}$ cortex), where $\mathrm{\Delta }x=1/L_x$ is the resolution limit of the spatial Fourier transform on an $L_x=L_y$ square grid. Crosses labeled (a)–(d) correspond to the four spectra of Fig. 13.