Bumps in Small-World Networks

1. Introduction

Spatially-localized “bumps” of activity in neuronal networks have been studied for many years, as they are thought to play a role in short term memory (Camperi and Wang, 1998; Compte et al., 2000; Bressloff, 2012; Wimmer et al., 2014) and the head direction system (Redish et al., 1996; Zhang, 1996), among other phenomena. Some models of bump formation have used a firing rate description (Wilson and Cowan, 1973; Amari, 1977; Laing et al., 2002; Laing and Troy, 2003; Owen et al., 2007) while others have considered networks of spiking neurons (Gutkin et al., 2001; Laing and Chow, 2001; Wang, 2001). The simplest models typically have “Mexican-hat” connectivity in a single population of neurons, where nearby neurons are excitatorily coupled and more distant ones are inhibitorily coupled (Ermentrout, 1998; Bressloff, 2012). However, more realistic models consider both excitatory and inhibitory neurons with non-negative connectivity within and between populations (Pinto and Ermentrout, 2001; Blomquist et al., 2005). Almost all previous models have considered homogeneous and isotropic networks, which typically support a continuous family of reflection-symmetric bumps, parameterized by their position in the network. Some exceptions are (Brackley and Turner, 2009, 2014), in which a spatially-inhomogeneous coupling function is used, and (Thul et al., 2016), in which a spatially-varying random firing threshold is imposed.

In this paper we further investigate the effects of breaking the spatial homogeneity of neural networks which support bump solutions, by randomly adding long-range connections and simultaneously removing short-range connections in a particular formulation of small-world networks (Song and Wang, 2014). Small-world networks (Watts and Strogatz, 1998) have been much studied and there is evidence for the existence of small-worldness in several brain networks (Bullmore and Sporns, 2009). In particular, we are interested in determining how sensitive networks which support bumps are to this type of random rewiring of connections, and thus how precisely networks must be constructed in order to support bumps.

We will consider networks of heterogeneous excitatory and inhibitory theta neurons, the theta neuron being the canonical model for a Type I neuron for which the onset of firing is through a saddle-node on an invariant circle bifurcation (Ermentrout and Kopell, 1986; Ermentrout, 1996). In several limits such networks are amenable to the use of the Ott/Antonsen ansatz (Ott and Antonsen, 2008, 2009), and we will build on previous work using this ansatz in the study of networks of heterogeneous theta neurons (Luke et al., 2013; Laing, 2014a, 2015; So et al., 2014). We present the model in Section 2.2 and then consider two limiting cases: an infinite number of neurons (Section 2.3) and an infinite ensemble of finite networks with the same connectivity (Section 2.4). Results are given in Section 3 and we conclude in Section 4. The Appendix contains some mathematical manipulations relating to Section 2.4.

2. Materials and Methods

2.1. Introduction

First consider an all-to-all coupled network of N heterogeneous theta neurons whose dynamics are given by

$\begin{array}{ll}\frac{d{\theta }_i}{dt}=1-\text{cos}{\theta }_i+(1+\text{cos}{\theta }_i)(I_i+gr)& (1)\end{array}$

$\begin{array}{ll}\tau \frac{dr}{dt}=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{P}_n({\theta }_j)-r& (2)\end{array}$

for i = 1, 2, …N where θ_i ∈ [0, 2π) is the phase of the ith neuron, $P_n(\theta )=a_n{(1-\text{cos}\theta )}^n,n\in {\mathbb{N}}^{+}$ and a_n is a normalization factor such that ${\displaystyle {\int }_0^{2\pi }}{P}_n(\theta )d\theta =2\pi$. The function P_n is meant to mimic the action potential generated when a neuron fires, i.e., its phase increases through π; n controls the “sharpness” of this function. The I_i are input currents randomly chosen from some distribution, g is the strength of connectivity within the network (positive for excitatory coupling and negative for inhibitory), and τ is a time constant governing the synaptic dynamics. The variable r is driven up by spiking activity and exponentially decays to zero in the absence of activity, on a timescale τ.

The model (1)–(2) with τ = 0 (i.e., instantaneous synapses) was studied by Luke et al. (2013), who found multistability and oscillatory behavior. The case of τ > 0 was considered in Laing, unpublished and similar forms of synaptic dynamics have been considered elsewhere (Börgers and Kopell, 2005; Ermentrout, 2006; Coombes and Byrne, unpublished). The model presented below results from generalizing Equations (1) and (2) in several ways. Firstly, we consider two populations of neurons, one excitatory and one inhibitory. Thus, we will have two sets of variables, one for each population. Such a pair of interacting populations was previously considered by Luke et al. (2014); Börgers and Kopell (2005); Coombes and Byrne, unpublished; and Laing, unpublished. Secondly, we consider a spatially-extended network, in which both the excitatory and inhibitory neurons lie on a ring, and are (initially) coupled to a fixed number of neurons either side of them. Networks with similar structure have been studied by many authors (Redish et al., 1996; Compte et al., 2000; Gutkin et al., 2001; Laing and Chow, 2001; Laing, 2014a, 2015).

2.2. Model

We consider a network of 2N theta neurons, N excitatory and N inhibitory. Within each population the neurons are arranged in a ring, and there are synaptic connections between and within populations, whose strength depends on the distance between neurons, as in Laing and Chow (2002) and Gutkin et al. (2001) (In the networks we will consider, connection strengths are either 1 or 0, i.e., neurons are either connected or not connected). Inhibitory synapses act on a timescale τ_i, whereas the excitatory ones act on a timescale, τ. θ_i ∈ [0, 2π) is the phase of the ith excitatory neuron and ϕ_i ∈ [0, 2π) is the phase of the ith inhibitory one. The equations are

$\begin{array}{ll}\frac{d{\theta }_i}{dt}=1-\text{cos}{\theta }_i+(1+\text{cos}{\theta }_i)(I_i+g_{EE}{v}_i-g_{EI}{y}_i)& (3)\end{array}$

$\begin{array}{ll}\frac{d{\varphi }_i}{dt}=1-\text{cos}{\varphi }_i+(1+\text{cos}{\varphi }_i)(J_i+g_{IE}{u}_i-g_{II}{z}_i)& (4)\end{array}$

$\begin{array}{ll}\tau \frac{d{v}_i}{dt}=r_i-v_i& (5)\end{array}$

$\begin{array}{ll}\tau \frac{d{u}_i}{dt}=q_i-u_i& (6)\end{array}$

$\begin{array}{ll}{\tau }_i\frac{d{y}_i}{dt}=s_i-y_i& (7)\end{array}$

$\begin{array}{ll}{\tau }_i\frac{d{z}_i}{dt}=w_i-z_i& (8)\end{array}$

for i = 1, 2…N, where

$\begin{array}{ll}{q}_i=\frac{1}{N}{\displaystyle \sum _{j=-M_{IE}}^{M_{IE}}}{P}_n({\theta }_{i+j})& (9)\end{array}$

$\begin{array}{ll}{r}_i=\frac{1}{N}{\displaystyle \sum _{j=-M_{EE}}^{M_{EE}}}{P}_n({\theta }_{i+j})& (10)\end{array}$

$\begin{array}{ll}{s}_i=\frac{1}{N}{\displaystyle \sum _{j=-M_{EI}}^{M_{EI}}}{P}_n({\varphi }_{i+j})& (11)\end{array}$

$\begin{array}{ll}{w}_i=\frac{1}{N}{\displaystyle \sum _{j=-M_{II}}^{M_{II}}}{P}_n({\varphi }_{i+j})& (12)\end{array}$

where P_n is as in Section 2.1. The positive integers M_IE, M_EE, M_EI, and M_II give the width of connectivity from excitatory to inhibitory, excitatory to excitatory, inhibitory to excitatory, and inhibitory to inhibitory populations, respectively. The non-negative quantities g_EE, g_EI, g_IE and g_II give the overall connection strengths within and between the two populations (excitatory to excitatory, inhibitory to excitatory, excitatory to inhibitory, and inhibitory to inhibitory, respectively). The variable v_i (when multiplied by g_EE) gives the excitatory input to the ith excitatory neuron, and whose dynamics are driven by r_i, which depends on the activity of the excitatory neurons with indices between i − M_EE and i + M_EE. Similarly, u_i (when multiplied by g_IE) gives the excitatory input to the ith inhibitory neuron, and is driven by q_i, which depends on the activity of the excitatory neurons with indices between i − M_IE and i + M_IE. g_EIy_i is the inhibitory input to the ith excitatory neuron, driven by s_i, which depends on the activity of the inhibitory neurons with indices between i − M_EI and i + M_EI. Lastly, g_IIz_i is the inhibitory input to the ith inhibitory neuron, driven by w_i, which depends on the activity of the inhibitory neurons with indices between i − M_II and i + M_II.

For simplicity, and motivated by the results in Pinto and Ermentrout (2001), we assume that the inhibitory synapses act instantaneously, i.e. τ_i = 0, and that there are no connections within the inhibitory population, i.e. g_II = 0. Thus (8) and (12) become irrelevant and from (7) we have that y_i = s_i in (3).

The networks are made heterogeneous by randomly choosing the currents I_i from the Lorentzian

$\begin{array}{ll}h(I)=\frac{\Delta /\pi }{{(I-I_0)}^2+{\Delta }^2}\hfill & (13)\hfill \end{array}$

and the currents J_i from the Lorentzian

$\begin{array}{ll}g(J)=\frac{\Delta /\pi }{{(J-J_0)}^2+{\Delta }^2}.& (14)\end{array}$

I₀ and J₀ are the centers of these distributions, and for simplicity we assume that both have the same width, Δ. The heterogeneity of the neurons (i.e., the positive value of Δ) is not necessary in order for the network to support bumps, but it is necessary for the Ott/Antonsen ansatz, used extensively below, to be valid (Ott et al., 2011). Networks of identical phase oscillators are known to show non-generic behavior which can be studied using the Watanabe/Strogatz ansatz (Watanabe and Strogatz, 1993, 1994). We want to avoid non-generic behavior, and having a heterogeneous network is also more realistic. For typical parameter values we see the behavior shown in Figures 1, 2, i.e., a stable stationary bump in which the active neurons are spatially localized.

FIGURE 1

Figure 1. A bump solution of Equations (3)–(6). Top: sin θ_i. Bottom: sin ϕ_i. Parameter values: N = 1024, Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

FIGURE 2

Figure 2. Average frequency for excitatory population (blue) and inhibitory (red) for the solution shown in Figure 1.

While these bumps may look superficially like “chimera” states in a ring of oscillators (Abrams and Strogatz, 2004, 2006; Laing, 2009; Panaggio and Abrams, 2015) they are different in one important aspect. Chimera states in the references above occur in networks for which the dynamics depend on only phase differences. Thus these systems are invariant with respect to adding the same constant to all oscillator phases, and can be studied in a rotating coordinate frame in which the synchronous oscillators have zero frequency, i.e., only relative frequencies are meaningful. In contrast, networks of theta neurons like those studied here are not invariant with respect to adding the same constant to all oscillator phases. The actual value of phase matters, and the neurons with zero frequency in Figure 2 have zero frequency simply because their input is not large enough to cause them to fire.

We now want to introduce rewiring parameters in such a way that on average, the number of connections is preserved as the networks are rewired. This is different from other formulations of small-world networks in which additional edges are added (Newman and Watts, 1999; Medvedev, 2014; but see Puljic and Kozma, 2008 for an example in which the number of connections to a node is precisely conserved). The reason for doing this is to keep the balance of excitation and inhibition constant. If we were to add additional connections, for example, within the excitatory population, the results seen might just be a result of increasing the number of connections, rather than their spatial arrangement. We are interested in the effects of rewiring connections from short range to long range, and thus use the form suggested in Song and Wang (2014). We replace Equations (9) to (11) by

$\begin{array}{l}{q}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{IE}{P}_n({\theta }_j);r_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EE}{P}_n({\theta }_j);\hfill \\ s_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EI}{P}_n({\varphi }_j);\hfill & (15)\end{array}$

where

$\begin{array}{ll}{A}_{ij}^{IE}=\{\begin{array}{ll}1\hfill & \text{with}\hfill \\ \hfill & \text{probability}\{\begin{array}{ll}1-[1-(2M_{IE}+1)/N]p_1,\hfill & |i-j|\le M_{IE}\hfill \\ (2M_{IE}+1)p_1/N,\hfill & |i-j|>M_{IE}\hfill \end{array}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}& (16)\end{array}$

$\begin{array}{ll}{A}_{ij}^{EE}=\{\begin{array}{ll}1\hfill & \text{with}\hfill \\ \hfill & \text{probability}\{\begin{array}{ll}1-[1-(2M_{EE}+1)/N]p_2,\hfill & |i-j|\le M_{EE}\hfill \\ (2M_{EE}+1)p_2/N,\hfill & |i-j|>M_{EE}\hfill \end{array}\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array}& (17)\end{array}$

and

$\begin{array}{ll}{A}_{ij}^{EI}=\{\begin{array}{ll}1\hfill & \text{with}\hfill \\ \hfill & \text{probability}\{\begin{array}{ll}1-[1-(2M_{EI}+1)/N]p_3,\hfill & |i-j|\le M_{EI}\hfill \\ (2M_{EI}+1)p_3/N,\hfill & |i-j|>M_{EI}\hfill \end{array}\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array}& (18)\end{array}$

where |i − j| refers to the shortest distance between neurons i and j, measured on the ring. When p₁ = p₂ = p₃ = 0, Equation (15) reverts to Equations (9) to (11). Note that when p₁ = 1, the probability of $A_{ij}^{IE}$ being 1 is independent of i and j, and that the expected number of nonzero entries in a row of $A_{ij}^{IE}$ (i.e., the expected number of connections from the excitatory population to an inhibitory neuron) is independent of p₁. Similar statements apply for the other two matrices and their parameters p₂ and p₃. Typical variation of A^IE with p₁ is shown in Figure 3 and it is clear that increasing p₁ interpolates between purely local connections (p₁ = 0) and uniform random connectivity (p₁ = 1).

FIGURE 3

Figure 3. Typical realizations of A^IE for p₁ = 0 (left) 0.5 (middle) and 1 (right). N = 1024, M_IE = 40. Black corresponds to a matrix entry of 1, white to 0.

We could simply simulate Equations (3)–(6) with Equation (15) for particular values of p₁, p₂, and p₃ but we would like to gain a deeper understanding of the dynamics of such a network. The first approach is to take the continuum limit in which the number of neurons in each network goes to infinity, in a particular way.

2.3. Continuum Limit

We take the continuum limit: N, M_EI, M_EE, M_IE → ∞ such that M_EI/N → α_EI, M_EE/N → α_EE and M_IE/N → α_IE, where 0 < α_EI, α_EE, α_IE < 1/2, and set the circumference of the ring of neurons to be 1. In this limit the sums (Equation 15) are replaced by integrals (more specifically, convolutions) with the connectivity kernels

$\begin{array}{ll}{G}_{IE}(x,p_1)=\{\begin{array}{ll}1-(1-2{\alpha }_{IE})p_1,\hfill & |x|<{\alpha }_{IE}\hfill \\ 2{\alpha }_{IE}{p}_1,\hfill & \text{otherwise}\hfill \end{array}& (19)\end{array}$

$\begin{array}{ll}{G}_{EE}(x,p_2)=\{\begin{array}{ll}1-(1-2{\alpha }_{EE})p_2,\hfill & |x|<{\alpha }_{EE}\hfill \\ 2{\alpha }_{EE}{p}_2,\hfill & \text{otherwise}\hfill \end{array}& (20)\end{array}$

$\begin{array}{ll}{G}_{EI}(x,p_3)=\{\begin{array}{ll}1-(1-2{\alpha }_{EI})p_3,\hfill & |x|<{\alpha }_{EI}\hfill \\ 2{\alpha }_{EI}{p}_3,\hfill & \text{otherwise}\hfill \end{array}& (21)\end{array}$

where G_IE(x, p₁) is the probability that a point in the excitatory population is connected to a point in the inhibitory population a distance x away, and similarly for the other two kernels. The effect of varying p_j, j = 1, 2, 3, on one of the functions Equations (19)–(21) is shown in Figure 4. Taking G_IE for example, we see that ${\displaystyle {\int }_{-1/2}^{1/2}}{G}_{IE}(x,p_1)dx=2{\alpha }_{IE}$ independent of p₁, i.e., the expected total number of connections is preserved, and similarly for the other two functions.

FIGURE 4

Figure 4. One of the functions Equations (19)–(21) with (A): p = 0, (B): p = 0.5, and (C): p = 1. For this example, α = 0.1. Note the similarity with the middle row of the matrices shown in Figure 3.

Taking the continuum limit of Equations (3)–(6) we describe the dynamics of the θ_i and ϕ_i in terms of probability densities F_E(θ, x, I, t) and F_I(ϕ, x, J, t), respectively, where x and t are (continuous) space and time, and I and J are random variables with densities h(I) and g(J) respectively. F_E satisfies the continuity equation (Luke et al., 2013)

$\begin{array}{ll}\frac{\partial F_E}{\partial t}+\frac{\partial }{\partial \theta }\left\{F_E\left[1-\text{cos}\theta +(1+\text{cos}\theta )(I+g_{EE}v-g_{EI}s)\right]\right\}=0& (22)\end{array}$

and similarly F_I satisfies

$\begin{array}{ll}\frac{\partial F_I}{\partial t}+\frac{\partial }{\partial \varphi }\left\{F_I\left[1-\text{cos}\varphi +(1+\text{cos}\varphi )(J+g_{IE}u)\right]\right\}=0& (23)\end{array}$

where

$\begin{array}{ll}\tau \frac{\partial v}{\partial t}=r-v& (24)\end{array}$

$\begin{array}{ll}\tau \frac{\partial u}{\partial t}=q-u& (25)\end{array}$

and

$\begin{array}{l}q(x,t)={\displaystyle {\int }_0^1}{G}_{IE}(|x-y|,p_1){\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}\hfill \\ F_E(\theta ,y,I,t)a_n{(1-\text{cos}\theta )}^{n}d\theta dIdy\hfill & (26)\end{array}$

$\begin{array}{l}r(x,t)={\displaystyle {\int }_0^1}{G}_{EE}(|x-y|,p_2){\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}\hfill \\ F_E(\theta ,y,I,t)a_n{(1-\text{cos}\theta )}^{n}d\theta dIdy\hfill & (27)\end{array}$

$\begin{array}{l}s(x,t)={\displaystyle {\int }_0^1}{G}_{EI}(|x-y|,p_3){\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}\hfill \\ F_I(\varphi ,y,J,t)a_n{(1-\text{cos}\varphi )}^{n}d\varphi dJdy\hfill & (28)\end{array}$

The forms of Equations (22) and (23) mean that they are amenable to the use of the Ott/Antonsen ansatz (Ott and Antonsen, 2008, 2009). This ansatz states that if the neurons are not identical (i.e., Δ > 0 for the networks studied here), solutions of the continuity equation [Equations (22) and (23)] decay exponentially onto a lower-dimensional manifold on which the θ and ϕ dependence of F_E and F_I, respectively, have a particular form. This form is a Fourier series in θ (or ϕ) in which the nth coefficient is some function to the nth power. [See Equation (A8), for example]. Thus, we can restrict Equations (22) and (23) to this manifold, thereby simplifying the dynamics.

The standard Kuramoto order parameter for an all-to-all coupled network with phases {θ_j} is the expected value of $e^{i{\theta }_j}$ (Strogatz, 2000). For the network studied here we can define the analogous spatially-dependent order parameters for the excitatory and inhibitory networks as

$\begin{array}{ll}{z}_E(x,t)={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}{F}_E(\theta ,x,I,t)e^{i\theta }d\theta dI& (29)\end{array}$

and

$\begin{array}{ll}{z}_I(x,t)={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}{F}_I(\varphi ,x,J,t)e^{i\varphi }d\varphi dJ& (30)\end{array}$

respectively. For fixed x and t, z_E(x, t) is a complex number with a phase and a magnitude. The phase gives the most likely value of θ and the magnitude governs the “sharpness” of the probability distribution of θ (at that x and t), and similarly for z_I(x, t) and ϕ (Laing, 2014a, 2015). We can also determine from z_E and z_I the instantaneous firing rate of each population (see Section 3.1 and Montbrió et al., 2015).

Performing manipulations as in Laing (2014a, 2015), Luke et al. (2013), and So et al. (2014) we obtain the continuum limit of Equations (3)–(6): evolution equations for z_E and z_I

$\begin{array}{l}\frac{\partial z_E}{\partial t}=\frac{(i{I}_0-\Delta ){(1+z_E)}^2-i{(1-z_E)}^2}{2}\\ +\frac{i{(1+z_E)}^2(g_{EE}v-g_{EI}s)}{2}& (31)\end{array}$

$\begin{array}{ll}\frac{\partial z_I}{\partial t}=\frac{(i{J}_0-\Delta ){(1+z_I)}^2-i{(1-z_I)}^2}{2}+\frac{i{(1+z_I)}^2{g}_{IE}u}{2}& (32)\end{array}$

together with Equations (24) and (25), where

$\begin{array}{ll}q(x,t)={\displaystyle {\int }_0^1}{G}_{IE}(|x-y|,p_1)H(z_E(y,t);n)dy& (33)\end{array}$

$\begin{array}{ll}r(x,t)={\displaystyle {\int }_0^1}{G}_{EE}(|x-y|,p_2)H(z_E(y,t);n)dy& (34)\end{array}$

$\begin{array}{ll}s(x,t)={\displaystyle {\int }_0^1}{G}_{EI}(|x-y|,p_3)H(z_I(y,t);n)dy& (35)\end{array}$

and

$\begin{array}{ll}H(z;n)=a_n\left[C_0+{\displaystyle \sum _{q=1}^n}{C}_q(z^q+{\overline{z}}^q)\right]& (36)\end{array}$

where

$\begin{array}{ll}{C}_q={\displaystyle \sum _{k=0}^n}{\displaystyle \sum _{m=0}^k}\frac{n!{(-1)}^k{\delta }_{k-2m,q}}{2^k(n-k)!m!(k-m)!}& (37)\end{array}$

and where by |x − y| in Equations (26)–(28) we mean the shortest distance between x and y given that they are both points on a circle, i.e., |x − y| = min (|x − y|, 1 − |x − y|).

The advantage of this continuum formulation is that bumps like that in Figure 1 are fixed points of Equations (31) and (32) and Equations (24) and (25). Once these equations have been spatially discretized, we can find fixed points of them using Newton's method, and determine the stability of these fixed points by finding the eigenvalues of the linearization around them. We can also follow these fixed points as parameter are varied, detecting (local) bifurcations (Laing, 2014b). The results of varying p₁, p₂ and p₃ independently are shown in Section 3.1.

2.4. Infinite Ensembles

We now consider the case where N is fixed and finite, and so are the matrices A^IE, A^EE and A^EI, but we average over an infinite ensemble of networks with these connectivities, where each member of the ensemble has a different (but consistent) realization of the random currents I_i and J_i (Barlev et al., 2011; Laing et al., 2012). This procedure results in 4N ordinary differential equations (ODEs), 2N of them for complex quantities and the other 2N for real quantities. Thus, there is no reduction of dimension from the original system (Equations 3–6), but as in Section 2.3, bump states will be fixed points of these ODEs.

Letting the number of members in the ensemble go to infinity, we describe the state of the excitatory network by the probability density function

$\begin{array}{ll}{f}^E({\theta }_1,{\theta }_2,\dots ,{\theta }_N;I_1,I_2,\dots I_N;t)\equiv f^E(\{\theta \};\{I\};t)& (38)\end{array}$

and that of the inhibitory one by

$\begin{array}{ll}{f}^I({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_N;J_1,J_2,\dots J_N;t)\equiv f^I(\{\varphi \};\{J\};t)& (39)\end{array}$

which satisfy the continuity equations

$\begin{array}{ll}\frac{\partial f^E}{\partial t}+{\displaystyle \sum _{j=1}^N}\frac{\partial }{\partial {\theta }_j}\left[f^E\left(\frac{d{\theta }_j}{dt}\right)\right]=0& (40)\end{array}$

and

$\begin{array}{ll}\frac{\partial f^I}{\partial t}+{\displaystyle \sum _{j=1}^N}\frac{\partial }{\partial {\varphi }_j}\left[f^I\left(\frac{d{\varphi }_j}{dt}\right)\right]=0& (41)\end{array}$

where dθ_j/dt and dϕ_j/dt are given by Equations (3) and (4).

Performing the manipulations in the Appendix we obtain

$\begin{array}{l}\frac{d{z}_j^E}{dt}=\frac{(i{I}_0-\Delta ){\left(1+z_j^E\right)}^2-i{\left(1-z_j^E\right)}^2}{2}\\ +\frac{i{(1+z_j^E)}^2(g_{EE}{v}_j-g_{EI}{s}_j)}{2}& (42)\end{array}$

$\begin{array}{ll}\frac{d{z}_j^I}{dt}=\frac{(i{J}_0-\Delta ){\left(1+z_j^I\right)}^2-i{\left(1-z_j^I\right)}^2}{2}+\frac{i{\left(1+z_j^I\right)}^2{g}_{IE}{u}_j}{2}& (43)\end{array}$

for j = 1, 2, …N where

$\begin{array}{ll}{q}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{IE}H\left(z_j^E(t);n\right)& (44)\end{array}$

$\begin{array}{ll}{r}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EE}H\left(z_j^E(t);n\right)& (45)\end{array}$

$\begin{array}{ll}{s}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EI}H\left(z_j^I(t);n\right)& (46)\end{array}$

and

$\begin{array}{ll}\tau \frac{d{v}_i}{dt}=r_i-v_i& (47)\end{array}$

$\begin{array}{ll}\tau \frac{d{u}_i}{dt}=q_i-u_i& (48)\end{array}$

for i = 1, 2, …N. Equations (42)–(48) form a complete description of the expected behavior of a network with connectivities given by the matrices A^IE, A^EE and A^EI. Note the similarities with Equations (31)–(35) and Equations (24) and (25). As mentioned above, the advantage of this formulation is that states like that in Figure 1 will be fixed points of Equations (42)–(48), for the specified connectivities.

Recalling that the matrices A^IE, A^EE and A^EI depend on the parameters p₁, p₂ and p₃ respectively we now investigate how solutions of Equations (42)–(48) depend on these parameters. One difficulty in trying to vary, say, p₁, is that the entries of A^IE do not depend continuously on p₁. Indeed, as presented, one should recalculate A^IE each time p₁ is changed. In order to generate results comparable with those from Section 2.3 we introduce a consistent family of matrices, following Medvedev (2014). Consider A^IE (similar procedures apply for the other two matrices) and define an N × N matrix r, each entry of which is independently and randomly chosen from a uniform distribution on the interval (0, 1). The matrix r is now considered to be fixed, and we define $A^{IE}(p_1)$ as follows:

$\begin{array}{ll}{A}_{ij}^{IE}(p_1)=\{\begin{array}{cc}\Theta [r_{ij}-p_1(1-(2M_{IE}+1)/N)],& |i-j|\le M_{IE}\\ \Theta [r_{ij}-(1-p_1(2M_{IE}+1)/N)],& |i-j|>M_{IE}\end{array}& (49)\end{array}$

where Θ is the Heaviside step function and the indices are taken modulo N. Comparing this with Equation (16) we see that for a fixed p₁, generating a new r and using Equation (49) is equivalent to generating A^IE using Equation (16). The reason for using Equation (49) is that since the r_ij are chosen once and then fixed, an entry in A^IE will switch from 0 to 1 (or vice versa) at most once as p₁ is varied monotonically in the interval [0, 1].

The effects of quasistatically increasing p₁ and p₃ for Equations (42)–(48) are shown in Section 3.2.

3. Results

3.1. Results for Continuum Limit

For the system Equations (31) and (32) and Equations (24) and (25) we discretize the spatial domain into 1024 evenly spaced points and approximate the integrals in Equations (33)–(35) with Riemann sums. We numerically integrate the spatially-discretized evolution equations in time, using appropriate initial conditions, until a steady state is reached. This steady state is then continued using pseudo-arclength continuation, and the stability of the solutions found determined by examining the eigenvalues of the Jacobian evaluated at them (Laing, 2014b). The increment between successive values of the p_i found during continuation is not fixed and the numerical results found were interpolated to a uniform grid for plotting in Figures 5, 7, 8. We consider varying p₁, p₂ and p₃ independently, keeping the other two parameters fixed at zero. The results of varying p₁ are shown in Figure 5, where we plot the firing rate of the two populations, derived as Re(w_i)/π where $w_i=(1-{\overline{z}}_i)/(1+{\overline{z}}_i)$ for i = I, E, as in Montbrió et al. (2015), where the z_i are fixed points of Equations (31) and (32). We see an increase and then decrease in bump width as p₁ is increased. There is also a pair of supercritical Hopf bifurcations, between which the bump is unstable (It is only weakly unstable, with the rightmost eigenvalue of the Jacobian having a maximal real part of 0.015 in this interval). At the leftmost Hopf bifurcation the Jacobian has eigenvalues ±1.8191i and at the rightmost it has eigenvalues ±1.7972i, with all other eigenvalues having negative real parts. One notable aspect is the increase in firing rate of the inhibitory population “outside” the bump as p₁ is increased, such that when p₁ = 1 the firing rate in this population is spatially homogeneous. This is to be expected, as there are no inhibitory-to-inhibitory connections, and when p₁ = 1 all inhibitory neurons receive the same input from the excitatory population.

FIGURE 5

Figure 5. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p₁, with p₂ = p₃ = 0. There is a Hopf bifurcation on both white vertical lines and the bump is unstable between these. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, α_IE = 40/1024, α_EE = 40/1024, α_EI = 60/1024 and τ = 10.

Increasing p₂ while keeping p₁ = p₃ = 0 we find that the bump undergoes a Hopf bifurcation (Jacobian has eigenvalues ±0.3404i) and then is destroyed in a saddle-node bifurcation at p₂ ≈ 0.48, as shown in Figure 6. The behavior of the bumps for 0 ≤ p₂ ≤ 0.48 is shown in Figure 7.

FIGURE 6

Figure 6. Maximum (over x) of the firing rate for the excitatory population as a function of p₂ with p₁ = p₃ = 0. Solid: stable; dashed: unstable. The Hopf bifurcation is marked with a circle. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, α_IE = 40/1024, α_EE = 40/1024, α_EI = 60/1024 and τ = 10.

FIGURE 7

Figure 7. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p₂, with p₃ = p₁ = 0. There is a Hopf bifurcation at the white vertical line and the bump is destroyed in saddle-node bifurcation at p₂ ≈ 0.48. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, α_IE = 40/1024, α_EE = 40/1024, α_EI = 60/1024 and τ = 10.

Varying p₃ we obtain Figure 8, where there are no bifurcations as p₃ is increased all the way to 1, corresponding to the case where all excitatory neurons feel the same inhibition, just a weighted mean of the output from the inhibitory population. We again see an increase and then slight decrease in bump width as p₃ is increased.

FIGURE 8

Figure 8. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p₃, with p₂ = p₁ = 0. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, α_IE = 40/1024, α_EE = 40/1024, α_EI = 60/1024 and τ = 10.

While a Hopf bifurcation of a bump may seem undesirable from a neurocomputational point of view, it should be kept in mind that oscillations are an essential phenomenon in many different neural networks, and they are widely studied (Ashwin et al., 2016).

We have only varied one of p₁, p₂ and p₃, keeping the other two probabilities at zero. A clearer picture of the system's behavior could be obtained by simultaneously varying two, or all three, of these probabilities. We leave this as future work, but mention that for the special case p₁ = p₂ = p₃ = p, the bump persists and is stable up to p ≈ 0.49, where it undergoes a saddle-node bifurcation (not shown).

3.2. Results for Infinite Ensemble

This section refers to Equations (42)–(48). In Figure 9 we show the results of slowly increasing p₁, while keeping p₂ = p₃ = 0. We initially set p₁ = 0 and integrated Equations (42)–(48) to a steady state, using initial conditions that give a bump solution. We then increased p₁ by 0.01 and integrated Equations (42)–(48) again for 10,000 time units, using as an initial condition the final state of the previous integration. We continued this process up to p₁ = 1. The firing rate for the jth excitatory neuron is Re(w_j)/π where $w_j=(1-{\overline{z}}_j^E)/(1+{\overline{z}}_j^E)$, and similarly for an inhibitory neuron. Comparing Figure 9 with Figure 5 we see the same behavior, the main difference being that the bump now moves in an unpredictable way around the domain as p₁ is increased. This is due to the system no longer being translationally invariant, and the bump moving to a position in which it is stable (Thul et al., 2016). Unlike the situation shown in Figure 5 we did not observe any Hopf bifurcations, for this realization of the A^IE. Presumably this is also a result of breaking the translational invariance and the weakly unstable nature of the bump shown in Figure 5 between the Hopf bifurcations.

FIGURE 9

Figure 9. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p₁, with p₂ = p₃ = 0. Compare with Figure 5. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, N = 1024, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

Repeating this process as p₃ is varied with p₁ = p₂ = 0 we obtain the results in Figure 10. Comparing with Figure 8 we see very good agreement, although the bump does move considerably for small p₃, as in Figure 9. Varing p₂ with p₁ = p₃ = 0 we obtain similar results to those in Figure 7 (not shown). Typical behavior of (42)–(48) with p₂ = 0.3 (i.e., beyond the Hopf bifurcation shown in Figure 7) is shown in Figure 11, where the oscillations are clearly seen.

FIGURE 10

Figure 10. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p₃, with p₂ = p₁ = 0. Compare with Figure 8. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, N = 1024, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

FIGURE 11

Figure 11. Instantaneous firing rate for (A): excitatory population and (B): inhibitory population, with p₂ = 0.3 and p₃ = p₁ = 0. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, N = 1024, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

3.3. Results for Original Network

To verify the results obtained above we ran the full network Equations (3)–(6) but using the connectivity Equation (49) (and similar constructions for A^EE and A^EI) to calculate Equation (15). The frequency was measured directly from simulations. Varying p₁ we obtain the results in Figure 12; again, no oscillatory behavior associated with a Hopf bifurcation was observed and the results are similar to those in Figure 9. Varying p₃ we obtained Figure 13 (compare with Figure 10). Figure 14 shows oscillatory behavior at p₂ = 0.3, p₁ = p₃ = 0 (the same parameter values as used in Figure 11).

FIGURE 12

Figure 12. Firing rate for (A): excitatory population and (B): inhibitory population in the full network Equations (3)–(6), averaged over a time window of length 500, as a function of p₁ with p₂ = p₃ = 0. Compare with Figure 9. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, N = 1024, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

FIGURE 13

Figure 13. Firing rate for (A): excitatory population and (B): inhibitory population in the full network Equations (3)–(6), averaged over a time window of length 500, as a function of p₃ with p₂ = p₁ = 0. Compare with Figure 10. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, N = 1024, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

FIGURE 14

Figure 14. Behavior of the full network Equations (3)–(6) with p₂ = 0.3 and p₃ = p₁ = 0. (A): 1 − cos θ_j, (B): v_j. Compare with Figure 11. Other parameters: Δ = 0.02, I₀ = −0.16, J₀ = −0.4, n = 2, g_EE = 25, g_IE = 25, g_EI = 7.5, N = 1024, M_IE = 40, M_EE = 40, M_EI = 60 and τ = 10.

Note that the results in Figures 9–14 are each for a single realization of a typical (parameterized) small-world network. To gain insight into general small-world networks it would be of interest to study the statistics of the behavior of such networks.

4. Discussion

We have considered the effects of randomly adding long-range and simultaneously removing short-range connections in a network of model theta neurons which is capable of supporting spatially localized bump solutions. Such rewiring makes the networks small-world, at least for small values of the rewiring probabilities. By using theta neurons we are able to use the Ott/Antonsen ansatz to derive descriptions of the networks in two limits: an infinite number of neurons, and an infinite ensemble of finite networks, each with the same connectivity. The usefulness of this is that the bumps of interest are fixed points of the dynamical equations derived in these ways, and can thus be found, their stability determined, and followed as parameters are varied using standard dynamical systems techniques.

For the parameters chosen we found bumps to be surprisingly robust: in several cases a rewiring probability could be taken from 0 to 1 without destroying a bump. However, rewiring connections within the excitatory population (increasing p₂) was found to destabilize a bump through a Hopf bifurcation and later destroy the unstable bump in a saddle-node bifurcation. Simulations of the full network were used to verify our results.

The network studied has many parameters: the spatial spread of local couplings, the timescale of excitatory synapses, the connection strengths within and between populations, and the distributions of heterogeneous input currents. These were all set so that the network without rewiring supported a stable bump solution, but we have not investigated the effects of varying any of these parameters. However, even without considering rewiring, Equations (31)–(35) and Equations (24) and (25) provide a framework for investigating the effects of varying these parameters on the existence and stability of bump solutions, since these continuum equations are derived directly from networks of spiking neurons, unlike many neural field models.

Author Contributions

All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest Statement

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

Acknowledgments

I thank the referees for their useful suggestions.

References

Abrams, D., and Strogatz, S. (2004). Chimera states for coupled oscillators. Phys. Rev. Lett. 93, 174102. doi: 10.1103/PhysRevLett.93.174102

PubMed Abstract | CrossRef Full Text | Google Scholar

Abrams, D., and Strogatz, S. (2006). Chimera states in a ring of nonlocally coupled oscillators. Int. J. Bifurcat. Chaos 16, 21–37. doi: 10.1142/S0218127406014551

CrossRef Full Text | Google Scholar

Amari, S.-I. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybern. 27, 77–87. doi: 10.1007/BF00337259

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashwin, P., Coombes, S., and Nicks, R. (2016). Mathematical frameworks for oscillatory network dynamics in neuroscience. J. Math. Neurosci. 6, 1–92. doi: 10.1186/s13408-015-0033-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Barlev, G., Antonsen, T. M., and Ott, E. (2011). The dynamics of network coupled phase oscillators: an ensemble approach. Chaos 21:025103. doi: 10.1063/1.3596711

PubMed Abstract | CrossRef Full Text | Google Scholar

Blomquist, P., Wyller, J., and Einevoll, G. T. (2005). Localized activity patterns in two-population neuronal networks. Physica D 206, 180–212. doi: 10.1016/j.physd.2005.05.004

CrossRef Full Text | Google Scholar

Börgers, C., and Kopell, N. (2005). Effects of noisy drive on rhythms in networks of excitatory and inhibitory neurons. Neural Comput. 17, 557–608. doi: 10.1162/0899766053019908

PubMed Abstract | CrossRef Full Text | Google Scholar

Brackley, C., and Turner, M. S. (2009). Persistent fluctuations of activity in undriven continuum neural field models with power-law connections. Phys. Rev. E 79:011918. doi: 10.1103/PhysRevE.79.011918

PubMed Abstract | CrossRef Full Text | Google Scholar

Brackley, C. A., and Turner, M. S. (2014). “Heterogeneous connectivity in neural fields: a stochastic approach,” in Neural Fields, eds S. Coombes, P. beim Graben, R. Potthast, and J. J. Wright (Berlin; Heidelberg: Springer), 213–234.

Bressloff, P. C. (2012). Spatiotemporal dynamics of continuum neural fields. J. Phys. A 45:033001. doi: 10.1088/1751-8113/45/3/033001

CrossRef Full Text | Google Scholar

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

Camperi, M., and Wang, X.-J. (1998). A model of visuospatial working memory in prefrontal cortex: recurrent network and cellular bistability. J. Comput. Neurosci. 5, 383–405. doi: 10.1023/A:1008837311948

PubMed Abstract | CrossRef Full Text | Google Scholar

Compte, A., Brunel, N., Goldman-Rakic, P. S., and Wang, X.-J. (2000). Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cereb. Cortex 10, 910–923. doi: 10.1093/cercor/10.9.910

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (1996). Type i membranes, phase resetting curves, and synchrony. Neural Comput. 8, 979–1001. doi: 10.1162/neco.1996.8.5.979

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (1998). Neural networks as spatio-temporal pattern-forming systems. Rep. Prog. Phys. 61, 353–430. doi: 10.1088/0034-4885/61/4/002

CrossRef Full Text | Google Scholar

Ermentrout, B. (2006). Gap junctions destroy persistent states in excitatory networks. Phys. Rev. E 74:031918. doi: 10.1103/PhysRevE.74.031918

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1986). Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46, 233–253. doi: 10.1137/0146017

CrossRef Full Text | Google Scholar

Gutkin, B. S., Laing, C. R., Colby, C. L., Chow, C. C., and Ermentrout, G. B. (2001). Turning on and off with excitation: the role of spike-timing asynchrony and synchrony in sustained neural activity. J. Comput. Neurosci. 11, 121–134. doi: 10.1023/A:1012837415096

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C., and Chow, C. (2001). Stationary bumps in networks of spiking neurons. Neural Comput. 13, 1473–1494. doi: 10.1162/089976601750264974

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2009). The dynamics of chimera states in heterogeneous Kuramoto networks. Physica D 238, 1569–1588. doi: 10.1016/j.physd.2009.04.012

CrossRef Full Text | Google Scholar

Laing, C. R. (2014a). Derivation of a neural field model from a network of theta neurons. Phys. Rev. E 90:010901. doi: 10.1103/PhysRevE.90.010901

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2014b). Numerical bifurcation theory for high-dimensional neural models. J. Math. Neurosci. 4, 13. doi: 10.1186/2190-8567-4-13

CrossRef Full Text | Google Scholar

Laing, C. R. (2015). Exact neural fields incorporating gap junctions. SIAM J. Appl. Dyn. Syst. 14, 1899–1929. doi: 10.1137/15M1011287

CrossRef Full Text | Google Scholar

Laing, C. R., and Chow, C. C. (2002). A spiking neuron model for binocular rivalry. J. Comput. Neurosci. 12, 39–53. doi: 10.1023/A:1014942129705

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R., Rajendran, K., and Kevrekidis, I. G. (2012). Chimeras in random non-complete networks of phase oscillators. Chaos 22:013132. doi: 10.1063/1.3694118

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R., and Troy, W. (2003). PDE methods for nonlocal models. SIAM J. Appl. Dyn. Syst. 2, 487–516. doi: 10.1137/030600040

CrossRef Full Text | Google Scholar

Laing, C. R., Troy, W. C., Gutkin, B., and Ermentrout, G. B. (2002). Multiple bumps in a neuronal model of working memory. SIAM J. Appl. Math. 63, 62–97. doi: 10.1137/S0036139901389495

CrossRef Full Text | Google Scholar

Luke, T. B., Barreto, E., and So, P. (2013). Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons. Neural Comput. 25, 3207–3234. doi: 10.1162/NECO_a_00525

PubMed Abstract | CrossRef Full Text | Google Scholar

Luke, T. B., Barreto, E., and So, P. (2014). Macroscopic complexity from an autonomous network of networks of theta neurons. Front. Comput. Neurosci. 8:145. doi: 10.3389/fncom.2014.00145

PubMed Abstract | CrossRef Full Text | Google Scholar

Medvedev, G. S. (2014). Small-world networks of kuramoto oscillators. Physica D 266, 13–22. doi: 10.1016/j.physd.2013.09.008

CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5:021028. doi: 10.1103/PhysRevX.5.021028

CrossRef Full Text | Google Scholar

Newman, M. E., and Watts, D. J. (1999). Renormalization group analysis of the small-world network model. Phys. Lett. A 263, 341–346. doi: 10.1016/S0375-9601(99)00757-4

CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. M. (2008). Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18:037113. doi: 10.1063/1.2930766

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. M. (2009). Long time evolution of phase oscillator systems. Chaos 19:023117. doi: 10.1063/1.3136851

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., Hunt, B. R., and Antonsen, T. M. (2011). Comment on “Long time evolution of phase oscillator systems” [Chaos 19, 023117 (2009)]. Chaos 21:025112. doi: 10.1063/1.3574931

PubMed Abstract | CrossRef Full Text | Google Scholar

Owen, M., Laing, C., and Coombes, S. (2007). Bumps and rings in a two-dimensional neural field: splitting and rotational instabilities. New J. Phys. 9, 378. doi: 10.1088/1367-2630/9/10/378

CrossRef Full Text | Google Scholar

Panaggio, M. J., and Abrams, D. M. (2015). Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators. Nonlinearity 28:R67. doi: 10.1088/0951-7715/28/3/R67

CrossRef Full Text | Google Scholar

Pinto, D. J., and Ermentrout, G. B. (2001). Spatially structured activity in synaptically coupled neuronal networks: II. lateral inhibition and standing pulses. SIAM J. Appl. Math. 62, 226–243. doi: 10.1137/S0036139900346465

CrossRef Full Text | Google Scholar

Puljic, M., and Kozma, R. (2008). Narrow-band oscillations in probabilistic cellular automata. Phys. Rev. E 78:026214. doi: 10.1103/PhysRevE.78.026214

PubMed Abstract | CrossRef Full Text | Google Scholar

Redish, A. D., Elga, A. N., and Touretzky, D. S. (1996). A coupled attractor model of the rodent head direction system. Network 7, 671–685. doi: 10.1088/0954-898X_7_4_004

CrossRef Full Text | Google Scholar

So, P., Luke, T. B., and Barreto, E. (2014). Networks of theta neurons with time-varying excitability: macroscopic chaos, multistability, and final-state uncertainty. Physica D 267, 16–26. doi: 10.1016/j.physd.2013.04.009

CrossRef Full Text | Google Scholar

Song, H. F., and Wang, X.-J. (2014). Simple, distance-dependent formulation of the Watts-Strogatz model for directed and undirected small-world networks. Phys. Rev. E 90:062801. doi: 10.1103/PhysRevE.90.062801

PubMed Abstract | CrossRef Full Text | Google Scholar

Strogatz, S. (2000). From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D 143, 1–20. doi: 10.1016/S0167-2789(00)00094-4

CrossRef Full Text | Google Scholar

Thul, R., Coombes, S., and Laing, C. R. (2016). Neural field models with threshold noise. J. Math. Neurosci. 6, 1–26. doi: 10.1186/s13408-016-0035-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J. (2001). Synaptic reverberation underlying mnemonic persistent activity. Trends Neurosci. 24, 455–463. doi: 10.1016/S0166-2236(00)01868-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, S., and Strogatz, S. (1993). Integrability of a globally coupled oscillator array. Phys. Rev. Lett. 70, 2391–2394. doi: 10.1103/PhysRevLett.70.2391

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, S., and Strogatz, S. (1994). Constants of motion for superconducting Josephson arrays. Physica D 74, 197–253. doi: 10.1016/0167-2789(94)90196-1

CrossRef Full Text | Google Scholar

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of “small-world” networks. Nature 393, 440–442. doi: 10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80. doi: 10.1007/BF00288786

PubMed Abstract | CrossRef Full Text | Google Scholar

Wimmer, K., Nykamp, D. Q., Constantinidis, C., and Compte, A. (2014). Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory. Nat. Neurosci. 17, 431–439. doi: 10.1038/nn.3645

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K. (1996). Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory. J. Neurosci. 16, 2112–2126.

PubMed Abstract | Google Scholar

Appendix

Mathematical Details Relating to Section 2.4

In the limit of an infinite ensemble we have

$\begin{array}{l}{q}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{IE}{\displaystyle \int }···{\displaystyle \int }{P}_n({\theta }_j)f^E(\{\theta \};\{I\};t)d{\theta }_{1}d{\theta }_2\dots d{\theta }_{N}d{I}_{1}d{I}_2\dots d{I}_N\hfill \\ =\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{IE}{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}{P}_n({\theta }_j)f_j^E({\theta }_j,I_j,t)d{\theta }_{j}d{I}_j\hfill & (\text{A}1)\end{array}$

$\begin{array}{l}{r}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EE}{\displaystyle \int }···{\displaystyle \int }{P}_n({\theta }_j)f^E(\{\theta \};\{I\};t)d{\theta }_{1}d{\theta }_2\dots d{\theta }_{N}d{I}_{1}d{I}_2\dots d{I}_N\hfill \\ =\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EE}{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}{P}_n({\theta }_j)f_j^E({\theta }_j,I_j,t)d{\theta }_{j}d{I}_j\hfill & (\text{A}2)\end{array}$

$\begin{array}{l}{s}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EI}{\displaystyle \int }···{\displaystyle \int }{P}_n({\varphi }_j)f^I(\{\varphi \};\{J\};t)d{\varphi }_{1}d{\varphi }_2\dots d{\varphi }_{N}d{J}_{1}d{J}_2\dots d{J}_N\hfill \\ =\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EI}{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}{P}_n({\varphi }_j)f_j^I({\varphi }_j,J_j,t)d{\varphi }_{j}d{J}_j\hfill & (\text{A}3)\end{array}$

where $f_j^E({\theta }_j,I_j,t)$ is the marginal distribution for θ_j, given by

$\begin{array}{ll}{f}_j^E({\theta }_j,I_j,t)={\displaystyle \int }···{\displaystyle \int }{f}^E(\{\theta \};\{I\};t){\displaystyle \prod _{k\ne j}}d{\theta }_{k}d{I}_k\hfill & (\text{A}4)\hfill \end{array}$

and similarly

$\begin{array}{ll}{f}_j^I({\varphi }_j,J_j,t)={\displaystyle \int }···{\displaystyle \int }{f}^I(\{\varphi \};\{J\};t){\displaystyle \prod _{k\ne j}}d{\varphi }_{k}d{J}_k& (\text{A}5)\end{array}$

Multiplying the continuity equation (40) by ∏_k≠j dθ_k dI_k and integrating we find that each $f_j^E$ satisfies

$\begin{array}{ll}\frac{\partial f_j^E}{\partial t}+\frac{\partial }{\partial {\theta }_j}\left[f_j^E\left(\frac{d{\theta }_j}{dt}\right)\right]=0& (\text{A}6)\end{array}$

Similarly each $f_j^I$ satisfies

$\begin{array}{ll}\frac{\partial f_j^I}{\partial t}+\frac{\partial }{\partial {\varphi }_j}\left[f_j^I\left(\frac{d{\varphi }_j}{dt}\right)\right]=0& (\text{A}7)\end{array}$

Using the Ott/Antonsen ansatz we write

$\begin{array}{ll}{f}_j^E({\theta }_j,I_j,t)=\frac{h(I_j)}{2\pi }\left\{1+{\displaystyle \sum _{n=1}^{\infty }}{\left[{\alpha }_j^E(I_j,t)\right]}^n{e}^{in{\theta }_j}+c.c.\right\}& (\text{A}8)\end{array}$

and

$\begin{array}{ll}{f}_j^I({\varphi }_j,J_j,t)=\frac{g(J_j)}{2\pi }\left\{1+{\displaystyle \sum _{n=1}^{\infty }}{\left[{\alpha }_j^I(J_j,t)\right]}^n{e}^{in{\varphi }_j}+c.c.\right\}& (\text{A}9)\end{array}$

for some functions ${\alpha }_j^E(I_j,t)$ and ${\alpha }_j^I(J_j,t)$, where “c.c.” means the complex conjugate of the previous term. Substituting Equation (A8) into Equation (A6) and Equation (A9) into Equation (A7) we find that

$\begin{array}{ll}\frac{\partial {\alpha }_j^E}{\partial t}=-i\left[\frac{I_j+g_{EE}{v}_j-g_{EI}{s}_j-1}{2}+(I_j+g_{EE}{v}_j-g_{EI}{s}_j+1){\alpha }_j^E+\frac{I_j+g_{EE}{v}_j-g_{EI}{s}_j-1}{2}{\left({\alpha }_j^E\right)}^2\right]& (\text{A}10)\end{array}$

and

$\begin{array}{ll}\frac{\partial {\alpha }_j^I}{\partial t}=-i\left[\frac{J_j+g_{IE}{u}_j-1}{2}+(J_j+g_{IE}{u}_j+1){\alpha }_j^I+\frac{J_j+g_{IE}{u}_j-1}{2}{\left({\alpha }_j^I\right)}^2\right]& (\text{A}11)\end{array}$

Substituting Equations (A8) and (A9) into Equations (A1)–(A3) we obtain

$\begin{array}{ll}{q}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{IE}{\displaystyle {\int }_{-\infty }^{\infty }}h(I_j)H({\alpha }_j^E(I_j,t);n)d{I}_j& (\text{A}12)\end{array}$

$\begin{array}{ll}{r}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EE}{\displaystyle {\int }_{-\infty }^{\infty }}h(I_j)H({\alpha }_j^E(I_j,t);n)d{I}_j& (\text{A}13)\end{array}$

$\begin{array}{ll}{s}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EI}{\displaystyle {\int }_{-\infty }^{\infty }}g(J_j)H({\alpha }_j^I(J_j,t);n)d{J}_j& (\text{A}14)\end{array}$

where H is given by Equation (36). Using standard properties of the Lorentzian one can perform the integrals in Equations (A12)–(A14) and defining $z_j^E(t)\equiv {\overline{\alpha }}_j^E(I_0+i\Delta ,t)$ and $z_j^I(t)\equiv {\overline{\alpha }}_j^I(J_0+i\Delta ,t)$ we have

$\begin{array}{ll}{q}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{IE}H(z_j^E(t);n)& (\text{A}15)\end{array}$

$\begin{array}{ll}{r}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EE}H(z_j^E(t);n)& (\text{A}16)\end{array}$

$\begin{array}{ll}{s}_i=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{A}_{ij}^{EI}H(z_j^I(t);n)& (\text{A}17)\end{array}$

Evaluating Equation (A10) at I_j = I₀ + iΔ and Equation (A11) at J_j = J₀ + iΔ we obtain

$\begin{array}{ll}\frac{d{z}_j^E}{dt}=\frac{(i{I}_0-\Delta ){\left(1+z_j^E\right)}^2-i{\left(1-z_j^E\right)}^2}{2}+\frac{i{(1+z_j^E)}^2(g_{EE}{v}_j-g_{EI}{s}_j)}{2}& (\text{A}18)\end{array}$

$\begin{array}{ll}\frac{d{z}_j^I}{dt}=\frac{(i{J}_0-\Delta ){\left(1+z_j^I\right)}^2-i{\left(1-z_j^I\right)}^2}{2}+\frac{i{\left(1+z_j^I\right)}^2{g}_{IE}{u}_j}{2}& (\text{A}19)\end{array}$

for j = 1, 2, …N. Equations (A15)–(A19) are Equations (42)–(46) in Section 2.4.