Complex Dynamics of Propagating Waves in a Two-Dimensional Neural Field

Propagating waves with complex dynamics have been widely observed in neural population activity. To understand their formation mechanisms, we investigate a type of two-dimensional neural field model by systematically varying its recurrent excitatory and inhibitory inputs. We show that the neural field model exhibits a rich repertoire of dynamical activity states when the relevant strength of excitation and inhibition is increased, ranging from localized rotating and traveling waves to global waves. Particularly, near the transition between stable states of rotating and traveling waves, the model exhibits a bistable state; that is, both the rotating and the traveling waves can exist, and the inclusion of noise can induce spontaneous transitions between them. Furthermore, we demonstrate that when there are multiple propagating waves, they exhibit rich collective propagation dynamics with variable propagating speeds and trajectories. We use techniques from time series analysis such detrended fluctuation analysis to characterize the effect of the strength of excitation and inhibition on these collective dynamics, which range from purely random motion to motion with long-range spatiotemporal correlations. These results provide insights into the possible contribution of excitation and inhibition toward a range of previously observed spatiotemporal wave phenomena.

1. Introduction

Propagating waves have been observed at different neural levels within multiple recording techniques, including multi-electrode arrays (Freeman and Barrie, 2000; Rubino et al., 2006; Muller et al., 2014; Townsend et al., 2015; Zanos et al., 2015), voltage sensitive dye (VSD) imaging (Wu et al., 2008; Huang et al., 2010; Muller et al., 2014), electroencephalography (EEG), electrocorticography (ECoG), magnetoencephalography (MEG) and functional magnetic resonance imaging (fMRI) (Goldman et al., 1949; Ribary et al., 1991; Lee et al., 2005; Patten et al., 2012; Alexander et al., 2013). A growing body of evidence indicates that these waves are not just an epiphenomenon but have functional significance: in spontaneous activity, propagating waves have been shown to follow repeated temporal motifs instead of occurring randomly (Mohajerani et al., 2013; Townsend et al., 2015). They are postulated to facilitate information transfer across brain regions (Rubino et al., 2006; Sato et al., 2012) and carry out distributed dynamical computation (Gong and Van Leeuwen, 2009). In sensory and motor cortices, stimuli can elicit repeatable propagating waves and the properties of these waves can be linked to stimulus features (Prechtl et al., 1997; Xu et al., 2007; Wu et al., 2008; Sato et al., 2012; Muller et al., 2014). For instance, the phase and amplitude of traveling waves in the motor cortex and visual cortex correlate with reach target location (Rubino et al., 2006) and with saccade size (Zanos et al., 2015), respectively, and the propagation direction of moving waves in the visual cortex is sensitive to visual movement orientation (Townsend et al., 2017). Waves in this visual cortex are also implicated in perceptual phenomena like binocular rivalry (Lee et al., 2005), reinforcement of recent visual experience (Han et al., 2008) and visual hallucinations (Ermentrout and Cowan, 1979). More broadly, they have been observed during sleep with a possible role in memory consolidation (Botella-Soler et al., 2012; Muller et al., 2016). This body of evidence indicates that understanding the mechanisms behind the formation and modulation of propagating neural waves is essential for uncovering the principled dynamics of neural population activity and for understanding the working mechanisms of neural circuits (Muller et al., 2018; Townsend and Gong, 2018).

Propagating neural waves in vivo exhibit complex dynamics with variable propagation speeds and trajectories (Freeman and Barrie, 2000; Ferezou et al., 2007; Luczak et al., 2007; Han et al., 2008). It has been shown that such propagating neural waves with complex dynamics can emerge from spatially-extended, spiking neural networks (Keane and Gong, 2015), and the collective complex dynamics of these waves provides a mechanistic explanation for a range of irregular neural dynamics, including the variability of spike timing, slow firing rate fluctuations (Churchland et al., 2010), and correlated membrane potential fluctuations (DeWeese and Zador, 2006; Tan et al., 2016). The relative strength of excitation and inhibition is important for the emergence of such complex dynamics (Keane and Gong, 2015). In this study, we use a recently developed 2-D neural field model (Qi and Gong, 2015) to further investigate the dynamical impact of the neurophysiological mechanism of the relative E/I strength on wave dynamics. By varying the relative strength of excitation and inhibition, we find that a variety of propagating waves can emerge from the neural field, including localized rotating and traveling waves, splitting waves, and global waves. Based on their dynamical properties, we map out their presence on the E/I parameter space. Near the transition from the rotating state to the traveling state, there are co-existing propagating waves (i.e., plane and rotating waves); the noise-induced switching dynamics between these waves are systematically characterized. We further demonstrate that when there are multiple propagating neural waves, their interactions give rise to a range of pattern dynamics, including waves of purely random, Brownian motion and waves with great fluctuating dynamics that have long-range spatiotemporal correlations. These collective dynamics are systematically characterized using a range of methods including detrended fluctuation analysis and spectral analysis.

2. Materials and Methods

2.1. Neural Field Model

We consider a previously introduced two-dimensional neural field model with refractoriness (Qi and Gong, 2015) which is derived from a three-state spiking neural circuit model (Gong and Robinson, 2012). The refractoriness in the neural field model is in the form of non-linear negative feedback which is analogous to models that feature linear recovery variables (Pinto and Ermentrout, 2001), refractoriness (Meijer and Coombes, 2014) and synaptic depression (Kilpatrick and Bressloff, 2010; Bressloff and Kilpatrick, 2011). This particular model supports a wide variety of localized spatiotemporal patterns such as bumps and propagating waves. In this study, we investigate the effects of changing excitation and inhibition on the dynamics of individual pattern formation and propagation as well as their collective interactions.

The model takes the form of a set of scalar integro-differential equations:

$\begin{array}{l}\tau \frac{\partial f}{\partial t}=-f+gH(u-\kappa )\\ \tau \frac{\partial g}{\partial t}=-gH(u-\kappa )+ph\\ \tau \frac{\partial h}{\partial t}=-ph+f,& (1)\end{array}$

where, f(r, t), g(r, t), and h(r, t) are temporal coarse-grained variables that describe the probability that a neuron located at r ∈ R² at time t is found in either a firing, refractory, or resting state, respectively. The parameter τ is the membrane time constant which we set to τ = 10ms, a typical value for cortical neurons (Koch, 2004). The parameter p expresses the refractoriness of the field: the rate at which a neuron transitions from a refractory state to a resting state. H(x) is the Heaviside step function: H(x) = 1 if x ≥ 0 and 0 otherwise. Here, it expresses the response of a neuron at r when its synaptic inputs u(r, t) exceed the threshold firing parameter κ:

$\begin{array}{ll}u(\text{r},t)=(w*f)(\text{r},t)={\displaystyle {\int }_{{\text{R}}^2}}w\left(|\text{r}-{\text{r}}^{\prime }\right)f({\text{r}}^{\prime },t)d{\text{r}}^{\prime },& (2)\end{array}$

where w(r) is a coupling function with short-range excitation and longer-range inhibition; as in Folias and Bressloff (2005), it is constructed from Bessel functions due to their analytic tractability:

$\begin{array}{l}w(r)=W_E{w}_K(r/{\sigma }_E)-W_I{w}_K(r/{\sigma }_I),\\ w_K(r)=\frac{2}{3\pi }\left[K_0(r)-K_0(2r)\right],& (3)\end{array}$

where K_i is the modified Bessel function of the second kind; W_E/W_I are the excitatory/inhibitory strength parameters, respectively; and σ_E/σ_I are their corresponding spatial scales. Figure 1 shows the weight function vs. r for a set of fixed parameter values which produce a shape that approximately matches that of the original spiking model (Gong and Robinson, 2012). As in that model, spatial scales σ_E, σ_I are chosen such that the size of the localized waves generated in the neural field matches that of waves observed in the cortex (Han et al., 2008).

FIGURE 1

Figure 1. (A) Synaptic feedback to the rest of the neural field due to activity at the origin, as expressed by the weight function (Equation 3) using typical parameter values W_E = 144.4, W_I = 73.7, σ_E = 0.187 mm, σ_I = 0.324 mm. The labeled circle with radius 2.79 represents the transition between the excitatory region (w(r) > 0) and the inhibitory region (w(r) < 0). Although inhibitory feedback is weaker in magnitude than excitatory feedback, it affects a much larger region of the field. We account for this integrating over the excitatory and inhibitory region to obtain our excitation and inhibition parameters (Equation 5). (B) Weight function w(r) (solid line) as a function of distance r for the same parameter values as in (A). The dotted line represents r × w(r), showing the synaptic feedback adjusted for the portion of the neural field receiving that feedback. The magnitudes of g_(±) are given by the areas of the red and blue shaded regions (Equation 5).

For convenience, t is rescaled to the membrane time τt → t. Further, by using the normalization condition f + g + h = 1, the neural field equations can be simplified to a pair of equations:

$\begin{array}{l}\frac{\partial f}{\partial t}=-f+(1-f-h)H(u-\kappa )\\ \frac{\partial h}{\partial t}=-ph+f,& (4)\end{array}$

2.1.1. Excitatory/Inhibitiory Parameters

In order to characterize the effect of varying excitation and inhibition on the activity within the neural field, we must first quantify the excitatory and inhibitory feedback as expressed by the weight function (Equation 3). Because this feedback varies with distance, we require a variable that reflects both the strength of the feedback and its spatial extent. A natural place to start is to consider the 2D integral of the weight function. For the limits, we inspect the weight function, denoted by the solid line in Figure 1B. We note that feedback is positive (excitatory) up to a certain distance we label r₀ and negative (inhibitory) beyond that. We thus introduce the following variables to quantify the excitatory and inhibitory feedback by integrating over two separate regions in space r:

$\begin{array}{l}{g}_{(+)}={\displaystyle \underset{\left|\text{r}\right|<r_0}{\iint }}w(r)d\text{r}\\ g_{(-)}={\displaystyle \underset{\left|\text{r}\right|>r_0}{\iint }}w(r)d\text{r}& (5)\end{array}$

The dotted line in Figure 1B represents the value of r × w(r). Hence, the area of the shaded regions represents the total excitatory and inhibitory inputs to the rest of the field. These quantities vary non-linearly and in a coupled way with respect to all four weight function parameters (the coupling strengths W_E, W_I and the spatial scales σ_E, σ_I). We proceed by using g_(±) as our excitatory/inhibitory parameters.

In order to numerically characterize the dynamical effects of excitatory and inhibitory feedback on the neural field, we need to properly sample the g_(±) parameter space. Since the parameters g_(±) vary non-linearly with respect to the weight function parameters, it is not possible to evenly sample the parameter space with a scheme of linear variation. Furthermore, a given pair of values g₍₊₎ and g₍₋₎ does not correspond to a unique weight function; different combinations of W_E, W_I, σ_E, and σ_I can produce similar but not identical weight functions with the same values of g₍₊₎ and g₍₋₎. In order to eliminate the sampling biases that can arise from a particular scheme of variation, we first use a Monte Carlo approach to sample the g_(±) parameter space by generating random values of W_E, W_I, σ_E and σ_I. We use this to explore general trends and validate the use of g_(±) as a measure of excitation/inhibition in the model. For more demanding calculations and simulations that would be unfeasible to do for all of the g_(±) pairs generated using the Monte Carlo approach, we instead use a smaller set of parameters, found by partitioning the g_(±) parameter space into an evenly-spaced grid. For each pair of g_(±) values, we then find corresponding values of W_E, W_I, σ_E, σ_I by applying a non-linear optimization algorithm to Equation (5). To do so, we introduce an objective function $f_{obj}=g_{(+)}^2+g_{(-)}^2$ which can then be minimized using one of many widely-available non-linear numerical solvers. We use the algorithm outlined in Byrd et al. (2000), although the results do not depend on the choice of algorithm.

For the Monte Carlo approach, we generate random weight function parameters by uniformly sampling each parameter (W_E, W_I, σ_E and σ_I) in the interval (0, 150). For each combination, we evaluate Equation (3) and discard combinations which do not result in a profile consisting of short-range excitation and long-range inhibition. To explore the dynamics of the field under varying excitation and inhibition, we first consider the effects on the simplest form of spatiotemporal activity: localized, radially symmetric bumps.

3. Results

3.1. Bump Solutions

3.1.1. Dependence of Radius on Excitation and Inhibition

Localized radially symmetric bumps are solutions to the neural field equations which have a constant, uniform value (f, h) within their circular boundaries. Their simple form results in a closed expression that is analytically tractable. For a bump of radius ρ centered at the origin, these solutions have the following form (Qi and Gong, 2015):

$\begin{array}{l}f(r)=\{\begin{array}{cc}\frac{p}{1+2p}& \text{for }r<\rho \\ 0\hfill & \text{for }r>\rho \end{array},\\ h(r)=\{\begin{array}{cc}\frac{1}{1+2p}& \text{for }r<\rho \\ 0\hfill & \text{for }r>\rho \end{array}.& (6)\end{array}$

At the boundary of the bump, the sum of synaptic input u(r, t) (Equation 2) must approach κ. This leads to the following existence condition for all bumps:

$\begin{array}{ll}\frac{p}{1+2p}\left[{\displaystyle {\int }_0^{2\pi }}{\displaystyle {\int }_0^{\rho }}w\left(|\text{r}-{\text{r}}^{\prime }|\right)r^{\prime }d{r}^{\prime }d{\theta }^{\prime }-2\kappa \right]-\kappa =0& (7)\end{array}$

We solve this equation numerically to find the bump radius ρ for a given weight function. Figure 2A shows the synaptic input at the bump boundary u as a function of the bump radius ρ, with the red line showing the value of the firing threshold κ. Bump solutions occur where the curve intersects the line. Due to the concave profile of u(ρ), there can be at most two solutions. The gradient of u(ρ) at the solutions implies that the smaller bump is unstable while the larger bump is stable: if the synaptic input at the boundary is negative, this results in a decrease in activity and thus a reduction in the bump radius; if the synaptic input is positive, this results in an increase in the bump radius. The positive gradient of u(ρ) at the smaller bump indicates that a small perturbation will result in either extinction or expansion to the larger bump solution, while the negative gradient of u(ρ) at the larger bump solution suggests that small perturbations will tend to return the bump to its original solution. This is in agreement with direct numerical simulations and is presented in Figures 2B,C.

FIGURE 2

Figure 2. (A) Synaptic input at the boundary of a bump as a function of bump radius (as given by Equation 7). Solutions are indicated by the arrows. The blue shaded region represents bumps whose synaptic input at the boundary is inhibitory, leading to contraction of the bump, while the red region represents bumps whose synaptic input at the boundary is positive, leading to expansion. (B,C) Representation of the boundaries of two bump solutions from (A) in the neural field (black circles). The arrows show the value of synaptic input at the boundary of bumps with different radii. Positive values (red arrows) result in expansion, while negative values (blue arrows) result in contraction. This implies that the smaller bump in (C) is unstable while the larger bump in (B,C) is stable.

Since the smaller bump solutions are unstable, we proceed by focusing on the larger branch of solutions. As shown in Figure 3A, the radius of the larger bumps depends on the parameters g_(±). Decreasing inhibition g₍₋₎ or increasing excitation g₍₊₎ results in an increase in the radius of the bumps. As experimentally observed waves in the cortex have a radius of around 2–4 mm (Han et al., 2008), we only explore the region of the g_± parameter space that corresponds to patterns of this size.

FIGURE 3

Figure 3. (A) Radius of large bump solutions as a function of g₍₊₎ and g₍₋₎. We fix p = 0.42 as this was approximately the critical value identified by Qi and Gong (2015) at which bump solutions were found to become unstable for the weight function parameters used in Figure 1. Increasing excitation or decreasing inhibition results in larger bumps. (B) Eigenvalues of large bump solutions with respect to shift perturbations are positive for all values of g₍₊₎ and g₍₋₎, becoming larger as excitation is increased or inhibition is lowered.

3.1.2. Dependence of Stability on Excitation and Inhibition

To quantify the influence that the level of excitation/inhibition has on these solutions, we perform a linearly stability analysis on the larger branch of solutions. We consider linear shift perturbations, as these correspond to the transition from stationary bumps to traveling waves. We follow the method of Qi and Gong (2015) to numerically evaluate the largest eigenvalue of the perturbation, as shown in Figure 3B. We find that eigenvalues are positive for all values of g_(±), increasing significantly as excitation is increased or inhibition is weakened. This indicates that for p = 0.42, radial bumps will always be unstable even if the level of excitation and inhibition is changed. Although Qi and Gong (2015) showed that lower values of p do result in stable bump solutions with negative eigenvalues for shift perturbations, we do not consider them as we wish to explore the collective dynamics of interacting traveling waves, which are only seen at higher values of p.

3.2. Propagating Waves

We now study the behavior of more complex patterns such as traveling waves under varying excitation and inhibition. Though we are interested in the behavior of waves in a two-dimensional neural field, we begin by considering the one-dimensional case in order to develop insights that are applicable to the 2D case and to validate subsequent observations made using direct numerical simulations. By considering the mathematical form of traveling waves we can quantify and semi-analytically characterize the variation we observe when varying the level of excitation and inhibition.

3.2.1. Traveling Waves in One-Dimension

When considering the 1D case, we can describe a traveling wave with just two parameters: it's length L and speed c. From Qi and Gong (2015), the equations of a traveling wave whose front lies at the origin and is traveling in the positive x direction are:

$\begin{array}{l}F(x)\text{ =}\{\begin{array}{ll}0,\hfill & \text{if }x>0\hfill \\ \frac{p}{1+2p}\left\{1-\exp (\alpha x/c)[\cos (\beta x/c)-A\sin (\beta x/c)]\right\},\hfill & \text{if }-L\le x\le 0\hfill \\ F(-L)\exp [(x+L)/c],\hfill & \text{if }x<-L\hfill \end{array}\\ H(x)\text{ =}\{\begin{array}{ll}0,\hfill & \text{if }x>0\hfill \\ \frac{1}{1+2p}\left\{1-\exp (\alpha x/c)[\cos (\beta x/c)-B\sin (\beta x/c)]\right\},\hfill & \text{if }-L\le x\le 0\hfill \\ H(-L)+\frac{F(-L)}{1-p}\exp [p(x+L)/c]-\frac{F(-L)}{1-p}\exp [(x+L)/c],\hfill & \text{if }x<-L,\hfill \end{array}& (8)\end{array}$

where $\alpha =(2+p)/2,\beta =\sqrt{p(4-p)}/2,A=(2-\alpha +1/p)/\beta$ and B = {α − 2− p[(α − 2)² + β²]}/β. Since the synaptic input u(x) at the two boundary points x ∈ {0, −L} must equal the threshold κ, we can find solutions by numerically solving the following equations over c and L:

$\begin{array}{l}u(0)={\displaystyle {\int }_{-\infty }^{0}w(x)F(x)dx}=\kappa \\ u(-L)={\displaystyle {\int }_{-\infty }^{0}w(x+L)F(x)dx}=\kappa & (9)\end{array}$

The two curves described by Equation (9) are shown in Figure 4A. They intersect at a single point, which we verify to be a solution using direct numerical simulation. By simultaneously solving both equations numerically for c and L, we can find traveling wave solutions for any value of g_±. To do so, we introduce a single objective function formed by summing the squares of Equation (9):

$f_{obj}={(u(0)-\kappa )}^2+{(u(-L)-\kappa )}^2)$

Using a multi-parameter optimization algorithm (Byrd et al., 2000), we find values of the wave speed c and length L that minimize f_obj. Figures 4B,C show the length and speed of traveling waves in 1 dimension as a function g₍₊₎ and g₍₋₎. We find that there exists a strong relationship between the strength of excitation/inhibition and the behavior of traveling waves: the speed and length of traveling waves increases as excitation increases or inhibition decreases. Figure 4D shows how the speed of the waves is correlated with their length. Most of the solutions lie on a curve, with a strong positive correlation between the two parameters. All the (c, L) solutions also lie on one side of this curve, which suggest that there is a frontier bounding possible (c, L) values and that there may be more than one solution for a given value of g_(±). This is also observed in direct numeric simulations, where there can be several traveling wave solutions for a given value of g_(±), although one solution will have greater velocity and length than all others.

FIGURE 4

Figure 4. (A) Method of finding the speed and length of waves in 1 dimension by simultaneously solving Equation (9). (B) Length of traveling waves in 1 dimension vs. the parameters g_(±) using p = 0.3. (C) Speed of traveling waves in 1 dimension vs. g_(±) using p = 0.3. (D) Wave speed vs. wave velocity in both 1D (blue) with p = 0.3 and 2D with p = 0.42 (red). The method used to obtain the results for 2D waves is outlined in section 3.3.

3.2.1.1. Two-dimensional case

We now consider the case of propagating waves in a two-dimensional neural field. For this, we must use direct numerical simulation as the inclusion of refractoriness in the neural field model means there are no simple methods for finding a closed-form traveling wave solution which would be necessary for an analytical approach. Direct simulation, however, still requires an appropriate initial condition for the neural field in order to simulate a wave, as they do not always emerge spontaneously. As with the one-dimensional case, the form of these waves (their shape and propagation speed) varies with the parameter values. Hence, we cannot use the same initial state for every simulation, particularly as for some parameter values, the neural field can support other types of activity such as rotating waves. In these cases, it is important that the initial state be very close to the traveling wave solution.

In order to generate suitable initial states, we develop a method to numerically find an approximate traveling wave solution for any given set of parameter values. We note that we can apply Equation (9) along the direction of propagation to find the value of (f, h) inside the wave boundary, but we must know both the speed of the wave and its shape (boundary) beforehand. We proceed by trying to find an empirical equation for the general boundary of a traveling wave. As a start, we use the shapes of waves that emerge spontaneously from an initial state consisting of a deformed bump solution. Varying p leads to various types of waves; these are presented in Figure 5A.

FIGURE 5

Figure 5. (A) The changing appearance of patterns as the neural field parameters are varied, either by changing p or the excitation/inhibition parameters g_±. As patterns become more unstable, they transition from rotating waves to traveling waves and unstable waves which undergo periodic division. (B) Empirically fitting a function curve to the wave front (green line) and tail (red line) for a traveling wave. Crosses show boundary points obtained from numeric simulation.

We first note that traveling waves have a characteristic teardrop or comet shape with a distinct front and tail which we separate by a line through their widest point. We thus propose the following general parametric form for the boundary of a traveling wave centered at the origin in polar coordinates, with ρ_trav(θ) being the boundary of the traveling wave and $\theta \in \left[-\frac{\pi }{2},\frac{3\pi }{2}\right)$:

$\begin{array}{ll}{\rho }_{trav}(\theta )=\{\begin{array}{cc}{\rho }_{front}(\theta )& \text{if }\theta \le \pi /2\\ {\rho }_{tail}(\theta )\hfill & \text{if }\theta >\pi /2\end{array}& (10)\end{array}$

This gives three global parameters: the front length L_f, tail length L_t and the width L_w. We find that the front is well described by the equation of an ellipse:

${\rho }_{front}(\theta )=L_f{L}_w/\sqrt{L_w^2{cos}^2\theta +L_f^2{sin}^2\theta }$

The tail, however, has a more complex shape. From testing various families of equations, we find that the following equation is a good fit to the boundary obtained from numeric simulations and only requires one additional parameter:

$\begin{array}{ll}{\rho }_{tail}(\theta )=L_w-L_a\left[1-1/(1-L_b{cos}^2\theta )\right].& (11)\end{array}$

L_a and L_b are parameters that describe the curvature of the tail and are related to the tail length L_t by: $L_t=L_w+L_a{L}_b/(1-L_b^2)$ and 0 < L_b < 1. Using Equations (8) and (10), we have a total of five parameters that describe a traveling wave uniquely: c, L_f, L_w, L_a, L_b. One such wave is shown in Figure 5B. In order to find traveling wave solutions for a given value of g_±, we must find parameters such that the synaptic input along the boundary of the wave equals the threshold value κ:

$\begin{array}{ll}u(\text{r}({\rho }_{trav}(\theta ),\theta ))={\displaystyle {\int }_{{\text{R}}^2}}w\left(|\text{r}-{\text{r}}^{\prime }|\right)F({\text{r}}^{\prime })d{\text{r}}^{\prime }=\kappa ,& (12)\end{array}$

where F is a traveling wave moving in the positive x direction:

$\begin{array}{l}F(\text{r})=F({\sigma }_{+}(y)){|}_{(L={\sigma }_{+}(y)-{\sigma }_{-}(y))},\end{array}$

and σ_±(y) is the x-coordinate of the front (+) or tail (−) of the wave for a given y-coordinate.

Unlike the one-dimensional case where we only had two discrete boundary conditions, here we have a continuous boundary. In order to solve Equation (12) numerically, we can discretize the boundary $\theta \in \left[-\frac{\pi }{2},\frac{3\pi }{2}\right)$ over n points and minimize the following objective function over a finite number of points θ_i on the boundary:

$\begin{array}{ll}{f}_{obj}(c,L_f,L_w,L_a,L_b)={\displaystyle \sum _{i=1}^n}{\left[u(\text{r}({\rho }_{trav},{\theta }_i))-\kappa \right]}^2.& (13)\end{array}$

This function is highly nonlinear with many local minima, so we apply a pattern search algorithm as outlined in Audet and Dennis (2002) in order to find values of c, L_f, L_w, L_a, L_b that minimize Equation (14). What we observe, however, is that this does not give us a unique solution: it is possible to find several similar but distinct traveling waves that minimize Equation (14). We verify the validity of these solutions using direct numerical simulation. We discretize spatially and temporally and apply a numerical scheme (Fourth-order Runge-Kutta as in Butcher, 2007) to the integro-differential equations (Equation 4). We use a spatial resolution of 0.1 mm with a grid size of 601 × 601 and a temporal resolution of 0.1 ms. We set κ to 1 and p to 0.42. In order to explore the parameter space, we linearly sample g_± pairs as outlined earlier in section 2.1. We simulate the neural field and then plot the final length and speed of the traveling waves on Figure 4D, alongside that of the one-dimensional waves. As with 1D waves, we observe that 2D traveling waves have a positive correlation between wave speed and velocity, with a similar “frontier” bounding possible values. When the wave velocity is low, we find that traveling waves tend to transition to rotating waves; when wave velocity is higher, they tend to split into pairs of waves moving in opposite directions. The length of the wave appears to be a crucial to this behavior: shorter waves are more likely to transition to rotating waves, while longer waves are stable up to a critical length at which local excitation in the tail of the wave becomes stronger than the inhibition from the front, resulting in the formation of a secondary front moving in the opposite direction. This critical length is approximately twice the radius of a stable bump solution.

3.3. Segmentation of Parameter Space Based on Behavior of Propagating Waves

Using the method outlined in section 3.2.1.1, we are now able to generate a traveling wave for any set of parameters. For many values of g_±, we find that the traveling waves are unstable and either periodically divide or transition to a rotating wave. We use this property of how individual waves behave in order to separate the parameter space into four distinct regions based on the dominant form of activity, as shown in Figure 6A. Region I - Rotating waves: Asymmetrical patterns which move in a circular orbit about a fixed point as in Figure 6D. Region II - Traveling waves: Symmetric patterns which move in linear trajectory with constant speed and length as in Figure 6E. Region III - Unstable waves: Traveling waves whose length and speed keeps increasing until they split into two waves propagating in opposite directions as seen in Figure 5A. Due to this, over time the field becomes filled with patterns that are spatially confined to a small region of the field as shown in Figure 6C. Region IV - Global waves: Extended, slowly-moving patterns with multiple “fronts” or bands of activation that fill the neural field, shown in Figure 6B.

FIGURE 6

Figure 6. (A) Characterization of the parameter space based on the dominant type of pattern seen in numeric simulations. (B–E) Snapshot images of the neural field at various parameter values. Center of mass trajectories (x_t, y_t) for isolated waves are overlaid in (D,E). (F,G) Velocity time series v_x, v_y, v_tot of the trajectories of individual waves as plotted in (D–E), found by differentiating x_t, y_t numerically. (H,I) Total acceleration time series a_tot found by differentiating v_tot numerically. The total velocity and acceleration of individual waves are both constant and positive for rotating waves while stable traveling waves have v_tot > 0 and a_tot = 0.

3.4. Characterization of Stability Using Direct Numerical Simulations

The non-linearity present in the neural field equations (Equation 4) makes it difficult to determine the stability of the traveling wave solutions analytically, so we use a computational approach to produce a phase portrait and identify its equilibria, similar to methods that have been applied to experimental data from unstable dynamical systems (Wiebe and Virgin, 2016). For our state variables, we track the trajectories of the waves using their center of mass coordinates, which we then use to calculate their velocity and acceleration over time. The instantaneous values of these variables can be used to distinguish between the different types of patterns. The variables also show a smooth transition between them in phase space. We outline these variables and how they describe different pattern types below.

The center-of-mass coordinates $\overrightarrow{\rho }(t_k)$ of the wave S at each time step t_k, are found using:

$\begin{array}{ll}{\displaystyle \sum _{\text{r}\in S}}(f(\text{r},t_k)-\overrightarrow{\rho }(t_k))=0& (14)\end{array}$

The wave velocities are numerically calculated using the finite difference method: v(t_k) = 0.5 × [r(t_k+1) − r(t_k−1)]/Δt. Smoothing is then applied using a 10-point moving average. Acceleration a(t_k) is calculated similarly from the velocity. The norm of the velocity and acceleration vectors is taken to calculate their total velocity v(t_k) = |r(t_k)| and total acceleration a(t_k) = |a(t_k)|.

Different pattern types have unique trajectories which give rise to varying velocity and acceleration profiles. For a single wave, the velocity in a fixed direction (e.g., along the x-axis) along with total velocity and acceleration are the most useful quantities for classifying them numerically. In Figure 6D, we show the center-of-mass trajectory of a rotating wave, with the figures below showing its velocity and acceleration. Isolated rotating waves have constant positive total acceleration and total velocity, i.e., a_tot > 0 and v_tot > 0. In Figure 6E, we repeat this for traveling waves, which have constant velocity in the direction of propagation. Hence, they have a_tot = 0 and v_tot > 0.

We now focus on the region near the boundary separating rotating waves and traveling waves to explore their stability and the nature of the transition between these two regimes. Figures 6D,E show a snapshot of the neural field at two values of g_(±) around the transition region. Figures 6F–I show the corresponding wave velocity and acceleration. We see that after a transient period approximately 50 ms, both waves stabilize with v_tot > 0 and a_tot > 0 for the rotating wave and v_tot > 0 and a_tot = 0 for the traveling wave. In Figure 7A, we show the trajectories of the waves shown in Figures 6D,E on a phase diagram with v_tot and a_tot on the x and y axes, respectively. To better explore the dynamical nature of the wave behavior, we initialize the neural field with a traveling wave and apply a single perturbation in the form of uniform Gaussian noise. From their initial state, we see a smooth transition toward their final state with several oscillations before they settle into a steady state.

FIGURE 7

Figure 7. (A) Velocity-acceleration phase portrait showing the trajectories of traveling waves after a single perturbation. The different paths correspond to different levels of inhibition near the critical point of transition between region I and II as shown in Figure 6A. Stationary bumps are represented by a single point at the origin; rotating waves have non-zero acceleration, hence they lie in the positive region of the y-axis; stable traveling waves have zero acceleration and lie along the x-axis. (B–D) Phase portraits of traveling waves with continuous noisy inputs. The green shaded regions indicate the location of stable solutions in the absence of noise. (B) At parameter values lying in region I (as outlined in Figure 6), the trajectory of waves in phase space oscillates around the region corresponding to stable rotating waves (shown in green). (C) At parameter values in between region I and II, the trajectory of waves in phase space oscillates between the regions corresponding to stable traveling waves and rotating waves. (D) At parameter values in region II, the trajectory only oscillates around the region corresponding to stable traveling waves.

At the critical point of transition between the two regimes, both types of waves coexist as bistable equilibria; rotating waves and traveling waves can be observed and a sufficiently strong perturbation can induce a transition from one to the other. To demonstrate this, we apply continuous perturbations to the traveling waves at each integration step, represented as an added noise term μ to Equation (4):

$\begin{array}{l}\frac{\partial f}{\partial t}=-f+(1-f-h)H(u-\kappa )+\mu \\ \frac{\partial h}{\partial t}=-ph+f-\mu \end{array}$

Where μ = f × ξ, and ξ is Gaussian white noise with zero mean and variance ${\sigma }_{\xi }^2$. Modeling noise in this way allows us to ignore the effect of noise on other aspects of the neural field, such as modifying the radius of bumps or spontaneously generating new patterns in other regions of the field.

Figures 7B–D show the resulting trajectories in phase space for the same parameter values as in Figure 7A as well as an additional parameter value that lies between these. For the noise term, we use σ_ξ = 0.2. The region of phase space corresponding to the stable value in the absence of noise highlighted in green. For parameter values drawn from the “Rotating waves” region of Figure 6A, the trajectory oscillates around the equilibrium corresponding to a stable rotating wave (Figure 7B). For a wave in the “Traveling Wave” region of parameter space (Figure 7D), it oscillates about the equilibria for traveling waves. If we adjust the parameters to lie in between these two regimes (Figure 7C), then the trajectory in phase space passes through both equilibria, demonstrating the bistability of the system.

The probability of a transition occurring is dependent on σ_ξ, as shown in Figure 8, which shows the average lifetime of each state when the neural field is simulated for 60 ms for the same parameter value as Figure 7C. At low levels of σ_ξ, traveling waves do not transition to rotating waves, but as σ_ξ is increased above 0.26, transitions happen with increasing frequency, although rotating waves appear more often than traveling waves.

FIGURE 8

Figure 8. Average time spent in either traveling or rotating wave state in terms of the strength of noise inputs. The initial state was a traveling wave. When noise is weak (σ_ξ < 0.026), no transitions occur. As noise is increased (σ_ξ > 0.026), traveling waves transition to rotating waves, which have a longer lifetime. When noise is very strong (σ_ξ > 0.06), the waves do not settle into single configuration and their trajectories explore the full state space.

3.5. Characterization of Collective Dynamics of Interacting Patterns

In addition to the behavior of isolated patterns, we study the impact of varying excitatory and inhibitory strength on the collective behavior of multiple interacting patterns and quantify these differences using the statistical properties of their motion. In order to reduce the variance that may arise due to the spatial discretization, field size, and periodic boundary conditions, we fix the size of the neural field and the size of patterns while varying the excitation and inhibition. As the value g_± influences the size of the patterns, as shown earlier in Figure 3A, we simulate only a subset of all the possible g_± values. This subset was found by calculating the upper bump radius for the g_± values obtained using the Monte Carlo method using Equation (7) and then filtering only those whose bump radius corresponds to around (3.00±0.01) mm, consistent with the size of experimentally observed waves in the cortex (Han et al., 2008). These values lie near a curved region in the parameter space, shown in Figure 9. As the size of rotating and traveling waves is very close to that of radial bumps, this ensures that the pattern size remains constant over the range of parameters.

FIGURE 9

Figure 9. Different types of collective behavior emerge when the excitatory/inhibitory strength is varied. (A) Dots indicate parameter values used for simulations of collective activity. Values were chosen so that the size of patterns remained relatively constant (as measured by the radius of bump solutions given by Equation 7). (B–G) Snapshots of neural field activity at various parameter values show the change from global waves and Turing patterns (B) to interacting bumps of irregular sizes (C), bumps of similar size aligned to a grid (D), interacting traveling waves (E) and isolated rotating waves (F,G). (H–L) Velocity time series v_x of a single wave along a fixed direction shows the changing characteristics of wave movement as parameter values are changed, from irregular movement driven by frequent collisions (H–J) to regular oscillations (F,G). (M–Q) Periodogram (Fourier transform) of the time series v_x shows a transition from a relatively flat power spectrum (M) to one with a clear peak corresponding to the frequency of the rotating waves' oscillation (K,L).

For the initial state of the simulation, we fix the number of patterns present. If there are too few patterns, they interact rarely and the resulting behavior is similar to that of isolated patterns. There is also a maximum number of patterns that can exist in the neural field at the same time, which depends on the size of the patterns and the strength of the long-range inhibition. For patterns of length 6 mm with a square grid length of 60 mm, this limit is in the range of 9–14 patterns. At this limit, the patterns will be spatially confined to a grid due to repulsive interactions arising from long-range inhibition. In order to observe and quantify the changing nature of the interactions between the patterns, we simulate a constant, fixed amount; we chose 7 so that interactions between patterns would occur regularly but they would initially not be spatially confined to a grid.

We initialize the patterns as traveling waves with random orientations and positions. After simulating for sufficiently long time to eliminate transient activity (4s), we inspect the activity visually and calculate their trajectories through their center-of-mass coordinates. Figures 9C–G shows snapshots of the neural field and the motion of the patterns as g_(±) is varied. Shown in Figures 9F,G is the behavior at high levels of inhibition, with the spatiotemporal patterns following periodic, isolated rotating trajectories. As inhibition is weakened, the trajectories become more varied and dynamic due to increased interactions between the patterns and the increasing stability of traveling waves (Figure 9E). As inhibition is weakened further, past the transition point where traveling waves become unstable (Figures 9C,D), the waves multiply and divide, filling the neural field and becoming confined to isolated positions within a hexagonal grid. As inhibition is weakened further, near the region where global waves dominate (Figure 9C), the patterns begin to show much more diversity in their size and appearance and are no longer confined to a grid. In the region where global waves dominate (Figure 9B), we no longer see isolated bump-like patterns but long, slowly-moving wave fronts.

To quantify these differences, we use the trajectories of the waves over time for time series analysis. We calculate the center-of-mass coordinates and velocities as previously described using Equation (15). These are shown in Figures 9H–L. We apply a Hanning window to the velocity data then take a Discrete Fourier Transform (DFT), after a sufficient transient period (3s). We calculate the N normalized coefficients of the power spectrum S_i and use it to produce the periodograms in Figures 9M–Q. The regular oscillations seen in the states in Figures 9F,G are reflected in the peak in their respective periodograms around 9 Hz. This disappears as inhibition is reduced and the motion becomes more chaotic, as shown in Figure 9E. The other components of the spectrum, however, remain fairly consistent until the onset of chaotic behavior as seen in Figure 9C. For global wave activity as in Figure 9B, we did not calculate the wave velocity as the large spatial extent of the wave fronts made the center of mass coordinates a less useful metric for quantifying the wave behavior. To quantify the changes in spectral profile and the velocity time series, we use three different statistics: The degrees of freedom (DoF) of the power spectrum, the approximate entropy (ApEn) of the wave velocity time series, and the Detrended Fluctuation Analysis (DFA) exponent of the wave velocity. The DoF statistic measures the uniformity of the power spectrum (Vaillancourt and Newell, 2003) and is defined as:

$\begin{array}{ll}\text{DoF}=\frac{N{\left({\displaystyle {\sum }_i^N}{S}_i\right)}^2}{{\displaystyle {\sum }_i^N}{S}_i^2}& (15)\end{array}$

For a uniform power spectrum (white noise), the DoF is equal to 1, whereas for a delta-function shaped power spectrum, it approaches 0. The DoF statistic increases as inhibition is reduced (Figure 10A), indicating that the power spectrum of the motion of the waves becomes more uniform as the periodic components associated with the motion of isolated rotating waves disappear and the motion becomes more random.

FIGURE 10

Figure 10. (A) Spectral degrees of freedom parameter vs. the inhibitory strength for the motion of waves in Figure 9. This parameter quantifies the uniformity of the spectral distribution and shows that the distribution is uniform at low levels of inhibition and non-uniform when inhibition is high, reflecting the transition from random motion to periodic oscillation. (B) Approximate Entropy vs. the inhibitory strength. This quantifies the degree of irregularity in the motion of the waves which peaks for unstable waves, around the point of transition between similarly-shaped waves aligned to a grid (Figure 9D) and waves of varying size (Figure 9C). (C) Detrended fluctuation analysis of the time series of wave velocity for two parameter values, the exponent α is represented by the gradient of the trend line. A value of 1.5 corresponds to a Brownian motion process (random walk), while a value of 1 corresponds to 1/f noise; a value of 1/2 is associated with white noise. (D) DFA exponents α of the wave velocity for various parameter values. For most parameter values, the DFA lies in the range 1.5–2, consistent with a non-stationary process that can be modeled as fractional Brownian motion. At decreasing levels of inhibition, the motion becomes increasingly random, reflect by the DFA exponent decreasing and reaching a minimum of 0.5.

The ApEn quantifies the amount of regularity and unpredictability of the fluctuations in a time series. The ApEn returns a value between 0 and 2 which reflects the predictability of future values of a time series based on previous values. To calculate ApEn, we followed the method outlined in Pincus (1991). As shown in Figure 10B, the Approximate Entropy peaks around the transition between similarly-shaped waves aligned to a grid and waves of varying size, indicating that the motion at these parameter values is the most chaotic.

We next perform a Detrended Fluctuation Analysis (DFA) on the time series v_x(t_k) to analyze the scale-free decay of temporal (auto)correlations, also known as long-range temporal correlations (LRTC). The DFA measures the power-law scaling of the root-mean-square fluctuation of the integrated and linearly detrended signals as a function of time window size. It is calculated as follows (Chen et al., 2002):

1. The velocity time series v_x(t_k) is integrated and the mean subtracted using $y(k)=\sum _{i=1}^k\left[v_x(t_i)-\langle v_x(t_i)\rangle \right]$.

2. The integrated signal is divided into boxes of equal length n.

3. In each box we fit a linear trend to the time series which we denoted y_n(k).

4. The integrated signal is detrended by subtracting the local trend y_n.

5. For a given box size n, the root-mean-square (rms) fluctuation for this signal is calculated:

$\begin{array}{ll}{F}_{DFA}(n)=\sqrt{\frac{1}{N_{max}}{\displaystyle \sum _{k=1}^{N_{max}}}{\left[y(k)-y_n(k)\right]}^2}.& (16)\end{array}$

6. The above calculation is repeated for a broad range of n to give a relationship between F_DFA(n) and n.

7. A power-law exponent α is fitted to the fluctuation function $F_{DFA}(n)n^{\alpha }$.

The DFA exponent α is the slope of the fluctuation function and can be related to the power-law scaling exponent of the autocorrelation function. For a particular time series, a DFA exponent of 0.5 indicates a signal that is without autocorrelations (white noise), whereas a DFA exponent between 0.5 and 1.0 indicates the presence of scale-free temporal correlations (autocorrelations) (Linkenkaer-Hansen et al., 2001; Chen et al., 2002; Gao et al., 2006). A value of 1 is associated with 1/f noise (pink noise); a value above 1 is associated with an unbounded, non-stationary signal, with a value near 1.5 associated with Brownian noise. We follow the method outlined in Linkenkaer-Hansen et al. (2001) to calculate the DFA exponents. The results for two distinct states - rotating waves and traveling waves - are presented in Figure 10C to highlight how this metric quantifies the changing behavior of the neural field.

We find that most states where characterized by an exponent slightly above 1.5, however this decreased suddenly for unstable states (Figure 10D), indicating that the motion becomes increasingly random and uncorrelated. The decrease in α as inhibition decreases can be explained by considering the impact of decreasing inhibition on the stability of patterns. This decrease in stability allows for a much wider range of accessible states and interactions. At any one point in time t_k, a diverse range of pattern shapes and sizes can be seen, as seen in Figure 6E.

3.6. Characterization of Long-Range Spatiotemporal Correlations Using Spectral Analysis

We then study the spatiotemporal spectrum of the neural field activity. To do so, we simulate the neural field for 4 s (3,000 time steps) in order to eliminate transient states and apply a Discrete Fourier Transform (DFT) to the final portion of the neural field output using the FFTW library (Frigo and Johnson, 2005). As before, our choice of initial condition is 7 randomly distributed traveling waves. This configuration allows us to observe the collective dynamics of patterns in the regime where traveling waves interact regularly but still have some freedom to move (i.e., not confined to a grid due to their repulsive interactions). As we have periodic boundary conditions in the spatial dimensions, we only apply a Hann (Hanning) window function along the temporal dimension, as outlined in Oppenheim et al. (1999). We take the absolute value of the complex output and, assuming rotational invariance, we average the spatial dimensions radially similar to Dong and Atick (1995). As a result, the final spectrum has one spatial and one temporal dimension.

We present the spectrum for two parameter values in Figures 11A,C. Alongside each spectrum is a snapshot image of the corresponding neural field (Figures 11B,D). The parameter values chosen lie on either side of the transition from stable to unstable traveling waves. In both spectra, the power peaks at the origin and is concentrated along lines of constant spatial/temporal frequency ratios, i.e., $\frac{\omega }{f}$. This can be explained by considering the Fourier transform of a single object moving with velocity v. As demonstrated by Watson and Ahumada (1985), the resulting spectrum is a distortion of the power spectrum of a stationary object which lies along a line described by:

$\begin{array}{ll}vf+\omega =0.& (17)\end{array}$

For multiple objects moving in two dimensions with a range of velocities, the resulting spectrum reflects the distribution of velocities among the objects (Dong and Atick, 1995). Following the approach of Rivait and Langer (2007), we explore this by calculating the total power along lines of constant velocity i.e., $v=\frac{\omega }{f}$. That is, for a range of velocity values v and temporal frequencies ω, we linearly interpolate the power spectrum Ŝ above and below the line $f=\frac{\omega }{v}$ and sum the resulting values:

$\begin{array}{ll}{Ŝ}_{tot}={\displaystyle \sum _{\omega }}Ŝ(v,\omega )& (18)\end{array}$

We show the result in Figure 12A and observe that for both parameter values, the power peaks at a particular value of v and decays away from this value. The peak value and the corresponding v value varies between different parameter values, reflecting the changing distribution of velocities due to interactions between patterns and their maximum velocity. We characterize this by simulating the neural field for the same range of parameter values as in the previous section and calculating their spectra. In Figure 12B, we show the peak power values as a function of the parameter g₍₋₎, with the color scale showing the corresponding velocity. We observe that the power and velocity peaks near the previously identified transition point between stable and unstable traveling waves. This reflects that near this point, waves in the neural field are most mobile; their velocity is maximized and they are able to propagate with fewer interactions with other waves than when the neural field is filled with patterns.

FIGURE 11

Figure 11. Spatiotemporal spectra for two parameter values alongside a snapshot of the neural field activity. (A,B) Unstable waves confined to grid. Most of the spectral power is concentrated at lower frequencies and falls off at higher frequencies. (C,D) Stable traveling waves. Spectral power is also concentrated at lower frequencies, however, the power drops off rapidly once it crosses a line of constant f/ω, which is consistent with the maximal velocity of the waves.

FIGURE 12

Figure 12. (A) Total spectral power vs. velocity, expressed as the ratio of temporal and spatial frequency. (B) The peak value of velocity-power curves as in (A) vs. the parameter g₍₋₎. The color scale represents the velocity value of the peak, which increases near the transition between stable and unstable traveling waves.

4. Discussion

In this study we have investigated the neural mechanism of the formation of complex propagating wave dynamics in a 2D neural field model. We have found that, as excitation increases or inhibition decreases, the dominant form of activity transitions from stationary bumps to rotating waves, to traveling waves, to unstable waves, and finally to global waves. These transitions can be used to define distinct regions in the parameter space. Using linear stability analysis, we verified these results by showing that an increase in excitation and a decrease in inhibition both resulted in the system moving further away from the region where stationary bumps are stable. At the transition between these regions, we have observed bi-stability of different kinds of spatio-temporal patterns: under isolated conditions, each pattern is stable, but perturbations can induce transitions between the states. Furthermore, continuous switching between stable states is possible in the presence of constant external noise.

This provides insight into the underlying mechanisms of the behavior of spatiotemporal patterns seen in experiments, such as switching between different types of activity patterns in local cortical circuits. This type of switching has been observed between spiral, plane, ring and irregular waves in disinhibited mammalian cortex (Huang et al., 2004), as well as a range of source/sink, spiral, saddle, and other patterns in primate and mouse cortex (Townsend et al., 2015). The exact mechanisms behind these changes are still unclear, however our results show that two possible factors are the effects of external noise and short-term modulations in the level of excitation or inhibition which could shift the phase landscape and induce a transition between states. In our neural field model, the level of excitation and inhibition is given by the weight function (Equation 3), which expresses the synaptic inputs to a point in the field based on activity in the rest of the field. This is affected by the synaptic connectivity, as well as the action of neurotransmitters, external inputs, and neural “noise.” There are several mechanisms by which this could be modified, including synaptic plasticity (Sussillo and Abbott, 2009; Vogels et al., 2012), the influence of endogenous substances and drugs (Behrends and ten Bruggencate, 2017), and changes in external inputs due to shifts in behavioral states (Taub et al., 2013; Zhou et al., 2014). These various mechanisms operate across a range of timescales that vary over several orders of magnitude, which could help account for the diverse range of spatiotemporal phenomena, such as the stable, rhythmic waves seen in human primary motor cortex (Takahashi et al., 2011). Future studies could experimentally verify these mechanisms and compare to the output of our model, similar to the work of Golomb and Amitai (1997) and Richardson et al. (2005).

In addition to single patterns, we have also studied the collective dynamics of multiple interacting patterns and the influence of the level of excitation/inhibition. Analyzing the statistical properties of the collective motion reveals a range of collective propagation dynamics, including interacting propagating waves with Brownian motion, random motion within patterns spatially confined to a grid, irregular patterns with unstable dynamics and slowly moving global waves. The quality of the motion, as revealed by statistical methods such as DFA, is influenced by the level of excitation and inhibition and shows dramatic changes at particular transitions in state space, such as when traveling waves become unstable. Long-range spatiotemporal correlations are also revealed by the spatial and temporal power spectrum of the neural field activity. These spatiotemporal correlations arise from the unique distribution of traveling waves in space and their motion within the neural field. Consistent with our analysis using other methods, we observe a shift in the spectrum as the activity changes due to changing levels of excitation and inhibition. This is linked to a change in the distribution of velocities among patterns in the neural field. Our results thus advance existing studies of neural fields which have mainly focused on single waves or multiple waves with regular dynamics (Kilpatrick and Bressloff, 2010; Bressloff and Kilpatrick, 2011; Meijer and Coombes, 2014).

The presence of balanced excitation and inhibition in cortical activity has been well established experimentally in vivo (Haider, 2006; Okun and Lampl, 2008; Xue et al., 2014) and it is believed to be essential for normal brain function. Our results show that changing levels of excitation and inhibition could be a mechanism for regulating the dynamics of spatiotemporal activity in neural systems.

Data Availability

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

This study was co-designed by DN and PG. DN performed the analysis and wrote the drafts of the manuscript with PG.

Funding

This research was supported by Australian Research Council, Discovery Grant, Grant ID: DP160104368.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Alexander, D. M., Jurica, P., Trengove, C., Nikolaev, A. R., Gepshtein, S., Zvyagintsev, M., et al. (2013). Traveling waves and trial averaging: the nature of single-trial and averaged brain responses in large-scale cortical signals. NeuroImage 73, 95–112. doi: 10.1016/j.neuroimage.2013.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Audet, C., and Dennis, J. E. (2002). Analysis of generalized pattern searches. SIAM J. Optimizat. 13, 889–903. doi: 10.1137/S1052623400378742

CrossRef Full Text | Google Scholar

Behrends, J. C., and ten Bruggencate, G. (2017). Cholinergic modulation of synaptic inhibition in the guinea pig hippocampus in vitro: excitation of GABAergic interneurons and inhibition of GABA-release. J. Neurophysiol. 69, 626–629. doi: 10.1152/jn.1993.69.2.626

PubMed Abstract | CrossRef Full Text | Google Scholar

Botella-Soler, V., Valderrama, M., Crépon, B., Navarro, V., and Le Van Quyen, M. (2012). Large-scale cortical dynamics of sleep slow waves. PLoS ONE 7:e30757. doi: 10.1371/journal.pone.0030757

PubMed Abstract | CrossRef Full Text | Google Scholar

Bressloff, P., and Kilpatrick, Z. (2011). Two-dimensional bumps in piecewise smooth neural fields with synaptic depression. SIAM J. Appl. Math. 71, 379–408. doi: 10.1137/100799423

CrossRef Full Text | Google Scholar

Butcher, J. (2007). Runge-Kutta methods. Scholarpedia 2:3147. doi: 10.4249/scholarpedia.3147

CrossRef Full Text | Google Scholar

Byrd, R. H., Gilbert, J. C., and Nocedal, J. (2000). A trust region method based on interior point techniques for nonlinear programming. Math. Progr. Ser. B 89, 149–185. doi: 10.1007/PL00011391

CrossRef Full Text | Google Scholar

Chen, Z., Ivanov, P. C., Hu, K., and Stanley, H. E. (2002). Effect of nonstationarities on detrended fluctuation analysis. Phys. Rev. E 65:15. doi: 10.1103/PhysRevE.65.041107

PubMed Abstract | CrossRef Full Text | Google Scholar

Churchland, M. M., Yu, B. M., Cunningham, J. P., Sugrue, L. P., Cohen, M. R., Corrado, G. S., et al. (2010). Stimulus onset quenches neural variability: a widespread cortical phenomenon Mark. Nat. Neurosci. 13, 369–378. doi: 10.1038/nn.2501

CrossRef Full Text | Google Scholar

DeWeese, M. R., and Zador, A. M. (2006). Non-gaussian membrane potential dynamics imply sparse, synchronous activity in auditory cortex. J. Neurosci. 26, 12206–12218. doi: 10.1523/JNEUROSCI.2813-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, D. W., and Atick, J. J. (1995). Statistics of natural time-varying images. Network 6, 345–358. doi: 10.1088/0954-898X_6_3_003

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Cowan, J. D. (1979). A mathematical theory of visual hallucination patterns. Biol. Cybernet. 34, 137–150. doi: 10.1007/BF00336965

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferezou, I., Haiss, F., Gentet, L. J., Aronoff, R., Weber, B., and Petersen, C. C. (2007). Spatiotemporal dynamics of cortical sensorimotor integration in behaving mice. Neuron 56, 907–923. doi: 10.1016/j.neuron.2007.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Folias, S. E., and Bressloff, P. C. (2005). Breathers in two-dimensional neural media. Phys. Rev. Lett. 95, 1–4. doi: 10.1103/PhysRevLett.95.208107

PubMed Abstract | CrossRef Full Text | Google Scholar

Freeman, W. J., and Barrie, J. M. (2000). Analysis of spatial patterns of phase in neocortical gamma EEGs in rabbit. J. Neurophysiol. 84, 1266–1278. doi: 10.1152/jn.2000.84.3.1266

PubMed Abstract | CrossRef Full Text | Google Scholar

Frigo, M., and Johnson, S. G. (2005). The design and implementation of FFTW3. Proc. IEEE 93, 216–231. doi: 10.1109/JPROC.2004.840301

CrossRef Full Text | Google Scholar

Gao, J., Hu, J., Tung, W. W., Cao, Y., Sarshar, N., and Roychowdhury, V. P. (2006). Assessment of long-range correlation in time series: how to avoid pitfalls. Phys. Rev. E 73, 1–10. doi: 10.1103/PhysRevE.73.016117

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, S., Santelmann, W. F., Vivian, W. E., and Goldman, D. (1949). Traveling waves in the brain - comments and communications. Science 109, 1–19. doi: 10.1126/science.109.2838.524

CrossRef Full Text | Google Scholar

Golomb, D., and Amitai, Y. (1997). Propagating neuronal discharges in neocortical slices: computational and experimental study. J. Neurophysiol. 78, 1199–211. doi: 10.1152/jn.1997.78.3.1199

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, P., and Robinson, P. A. (2012). Dynamic pattern formation and collisions in networks of excitable elements. Phys. Rev. E 85, 1–5. doi: 10.1103/PhysRevE.85.055101

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, P., and van Leeuwen, C. (2009). Distributed dynamical computation in neural circuits with propagating coherent activity patterns. PLoS Comput. Biol. 5:e1000611. doi: 10.1371/journal.pcbi.1000611

PubMed Abstract | CrossRef Full Text | Google Scholar

Haider, B. (2006). Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition. J. Neurosci. 26, 4535–4545. doi: 10.1523/JNEUROSCI.5297-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, F., Caporale, N., and Dan, Y. (2008). Reverberation of recent visual experience in spontaneous cortical waves. Neuron 60, 321–327. doi: 10.1016/j.neuron.2008.08.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X., Troy, W. C., Yang, Q., Ma, H., Laing, C. R., Schiff, S. J., et al. (2004). Spiral waves in disinhibited mammalian neocortex. J. Neurosci. 24, 9897–9902. doi: 10.1523/JNEUROSCI.2705-04.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X., Xu, W., Liang, J., Takagaki, K., Gao, X., and Wu, J. y. (2010). Spiral wave dynamics in neocortex. Neuron 68, 978–990. doi: 10.1016/j.neuron.2010.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Keane, A., and Gong, P. (2015). Propagating waves can explain irregular neural dynamics. J. Neurosci. 35, 1591–1605. doi: 10.1523/JNEUROSCI.1669-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kilpatrick, Z. P., and Bressloff, P. C. (2010). Spatially structured oscillations in a two-dimensional excitatory neuronal network with synaptic depression. J. Comput. Neurosci. 28, 193–209. doi: 10.1007/s10827-009-0199-

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, C. (2004). Neurons: Biophysics of Computation: Information Processing in Single Neurons. New York, NY: Oxford University Press.

Google Scholar

Lee, D., Ezhkova, E., Li, B., Pattenden, S. G., Tansey, W. P., and Workman, J. L. (2005). The proteasome regulatory particle alters the SAGA coactivator to enhance its interactions with transcriptional activators. Cell 123, 423–436. doi: 10.1016/j.cell.2005.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Linkenkaer-Hansen, K., Nikouline, V. V., Palva, J. M., and Ilmoniemi, R. J. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370–1377. doi: 10.1523/JNEUROSCI.21-04-01370.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Luczak, A., Barthó, P., Marguet, S. L., Buzsaki, G., and Harris, K. D. (2007). Sequential structure of neocortical spontaneous activity in vivo. Proc. Natl. Acad. Sci. U.S.A. 104, 347–352. doi: 10.1073/pnas.0605643104

PubMed Abstract | CrossRef Full Text | Google Scholar

Meijer, H. G. E., and Coombes, S. (2014). Travelling waves in a neural field model with refractoriness. J. Math. Biol. 68, 1249–1268. doi: 10.1007/s00285-013-0670-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohajerani, M. H., Chan, A. W., Mohsenvand, M., Ledue, J., Liu, R., McVea, D. A., et al. (2013). Spontaneous cortical activity alternates between motifs defined by regional axonal projections. Nat. Neurosci. 16, 1426–1435. doi: 10.1038/nn.3499

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, L., Chavane, F., Reynolds, J., and Sejnowski, T. J. (2018). Cortical travelling waves: mechanisms and computational principles. Nat. Rev. Neurosci. 19, 255–268. doi: 10.1038/nrn.2018.20

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, L., Piantoni, G., Koller, D., Cash, S. S., Halgren, E., and Sejnowski, T. J. (2016). Rotating waves during human sleep spindles organize global patterns of activity that repeat precisely through the night. eLife 5, 1–16. doi: 10.7554/eLife.17267

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, L., Reynaud, A., Chavane, F., and Destexhe, A. (2014). The stimulus-evoked population response in visual cortex of awake monkey is a propagating wave. Nat. Commun. 5:3675. doi: 10.1038/ncomms4675

PubMed Abstract | CrossRef Full Text | Google Scholar

Okun, M., and Lampl, I. (2008). Instantaneous correlation of excitation and inhibition during ongoing and sensory-evoked activities. Nat. Neurosci. 11, 535–537. doi: 10.1038/nn.2105

PubMed Abstract | CrossRef Full Text | Google Scholar

Oppenheim, A. V., Schafer, R. W., and Buck, J. R. (1999). Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice Hall.

Google Scholar

Patten, T. M., Rennie, C. J., Robinson, P. A., and Gong, P. (2012). Human cortical traveling waves: dynamical properties and correlations with responses. PLoS ONE 7:e38392. doi: 10.1371/journal.pone.0038392

PubMed Abstract | CrossRef Full Text | Google Scholar

Pincus, S. M. (1991). Approximate entropy as a measure of system complexity. Mathematics 88, 2297–2301. doi: 10.1073/pnas.88.6.2297

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, D. J., and Ermentrout, G. B. (2001). Spatially structured activity in synaptically coupled neuronal networks: I. Traveling fronts and pulses. SIAM J. Appl. Math. 62, 206–225. doi: 10.1137/S0036139900346453

CrossRef Full Text | Google Scholar

Prechtl, J. C., Cohen, L. B., Pesaran, B., Mitra, P. P., and Kleinfeld, D. (1997). Visual stimuli induce waves of electrical activity in turtle cortex. Proc. Natl. Acad. Sci. U.S.A. 94, 7621–7626. doi: 10.1073/pnas.94.14.7621

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, Y., and Gong, P. (2015). Dynamic patterns in a two-dimensional neural field with refractoriness. Phys. Rev. E 92, 1–13. doi: 10.1103/PhysRevE.92.022702

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribary, U., Ioannides, A. A., Singh, K. D., Hasson, R., Bolton, J. P., Lado, F., et al. (1991). Magnetic field tomography of coherent thalamocortical 40-Hz oscillations in humans. Proc. Natl. Acad. Sci. U.S.A. 88, 11037–11041. doi: 10.1073/pnas.88.24.11037

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, K. A., Schiff, S. J., and Gluckman, B. J. (2005). Control of traveling waves in the mammalian cortex. Phys. Rev. Lett. 94, 1–4. doi: 10.1103/PhysRevLett.94.028103

PubMed Abstract | CrossRef Full Text | Google Scholar

Rivait, D., and Langer, M. S. (2007). “Spatiotemporal power spectra of motion parallax: the case of cluttered 3D scenes,” in Human Vision and Electronic Imaging XII (San Jose, CA), 649208–649208–9.

Google Scholar

Rubino, D., Robbins, K. A., and Hatsopoulos, N. G. (2006). Propagating waves mediate information transfer in the motor cortex. Nat. Neurosci. 9, 1549–1557. doi: 10.1038/nn1802

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, T. K., Nauhaus, I., and Carandini, M. (2012). Traveling Waves in Visual Cortex. Neuron 75, 218–229. doi: 10.1016/j.neuron.2012.06.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Sussillo, D., and Abbott, L. F. (2009). Generating coherent patterns of activity from chaotic neural networks. Neuron 63, 544–557. doi: 10.1016/j.neuron.2009.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, K., Saleh, M., Penn, R. D., and Hatsopoulos, N. G. (2011). Propagating waves in human motor cortex. Front. Hum. Neurosci. 5:40. doi: 10.3389/fnhum.2011.00040

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, A. Y., Chen, Y., Scholl, B., Seidemann, E., and Priebe, N. J. (2016). Sensory stimulation shifts visual cortex from synchronous to asynchronous states. Nature 509, 226–229. doi: 10.1038/nature13159

PubMed Abstract | CrossRef Full Text | Google Scholar

Taub, A. H., Katz, Y., and Lampl, I. (2013). Cortical balance of excitation and inhibition is regulated by the rate of synaptic activity. J. Neurosci. 33, 14359–14368. doi: 10.1523/JNEUROSCI.1748-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Townsend, R., Solomon, S. S., Martin, P. R., Solomon, S. G., and Gong, P. (2017). Visual motion discrimination by propagating patterns in primate cerebral cortex. J. Neurosci. 37, 1538–1517. doi: 10.1523/JNEUROSCI.1538-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Townsend, R. G., and Gong, P. (2018). Detection and analysis of spatiotemporal patterns in brain activity. PLoS Comput. Biol. 14:e1006643. doi: 10.1371/journal.pcbi.1006643

PubMed Abstract | CrossRef Full Text | Google Scholar

Townsend, R. G., Solomon, S. S., Chen, S. C., Pietersen, A. N. J., Martin, P. R., Solomon, S. G., et al. (2015). Emergence of complex wave patterns in primate cerebral cortex. J. Neurosci. 35, 4657–4662. doi: 10.1523/JNEUROSCI.4509-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaillancourt, D. E., and Newell, K. M. (2003). Aging and the time and frequency structure of force output variability. J. Appl. Physiol. 94, 903–912. doi: 10.1152/japplphysiol.00166.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogels, T. H., Sprekeler, H., Zenke, F., Clopath, C., and Gerstner, W. (2012). Inhibitory plasticity balances excitation and inhibition in sensory. Science 334, 1569–1573. doi: 10.1126/science.1211095

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, A. B., and Ahumada, A. J. (1985). Model of human visual-motion sensing. J. Opt. Soc. Am. A 2:322. doi: 10.1364/JOSAA.2.000322

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiebe, R., and Virgin, L. N. (2016). On the experimental identification of unstable static equilibria. Proc. R. Soc. A Math. Phys. Eng. Sci. 472:20160172. doi: 10.1098/rspa.2016.0172

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J. Y., Huang, X., and Zhang, C. (2008). Propagating waves of activity in the neocortex: what they are, what they do. Neuroscientist 14, 487–502. doi: 10.1177/1073858408317066

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, W., Huang, X., Takagaki, K., and young Wu, J. (2007). Compression and reflection of visually evoked cortical waves. Neuron 55, 119–129. doi: 10.1016/j.neuron.2007.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, M., Atallah, B. V., and Scanziani, M. (2014). Equalizing excitation-inhibition ratios across visual cortical neurons. Nature 511, 596–600. doi: 10.1038/nature13321

PubMed Abstract | CrossRef Full Text | Google Scholar

Zanos, T. P., Mineault, P. J., Nasiotis, K. T., Guitton, D., and Pack, C. C. (2015). A sensorimotor role for traveling waves in primate visual cortex. Neuron 85, 615–627. doi: 10.1016/j.neuron.2014.12.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, M., Liang, F., Xiong, X. R., Li, L., Li, H., Xiao, Z., et al. (2014). Scaling down of balanced excitation and inhibition by active behavioral states in auditory cortex. Nat. Neurosci. 17, 841–850. doi: 10.1038/nn.3701

PubMed Abstract | CrossRef Full Text | Google Scholar