# Bumps in Small-World Networks

- Institute of Natural and Mathematical Sciences, Massey University, Auckland, New Zealand

We consider a network of coupled excitatory and inhibitory theta neurons which is capable of supporting stable spatially-localized “bump” solutions. We randomly add long-range and simultaneously remove short-range connections within the network to form a small-world network and investigate the effects of this rewiring on the existence and stability of the bump solution. We consider two limits in which continuum equations can be derived; bump solutions are *fixed points* of these equations. We can thus use standard numerical bifurcation analysis to determine the stability of these bumps and to follow them as parameters (such as rewiring probabilities) are varied. We find that under some rewiring schemes bumps are quite robust, whereas in other schemes they can become unstable via Hopf bifurcation or even be destroyed in saddle-node bifurcations.

## 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

for *i* = 1, 2, …*N* where θ_{i} ∈ [0, 2π) is the phase of the *i*th neuron, ${P}_{n}(\theta )={a}_{n}{(1-\text{cos}\theta )}^{n},n\in {\mathbb{N}}^{+}$ and *a*_{n} is a normalization factor such that ${\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 2*N* 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 *i*th excitatory neuron and ϕ_{i} ∈ [0, 2π) is the phase of the *i*th inhibitory one. The equations are

for *i* = 1, 2…*N*, where

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 *i*th 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 *i*th 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*_{EI}*y*_{i} is the inhibitory input to the *i*th 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*_{II}*z*_{i} is the inhibitory input to the *i*th 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

and the currents *J*_{i} from the Lorentzian

*I*_{0} and *J*_{0} 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. A bump solution of Equations (3)–(6)**. **Top:** sin θ_{i}. **Bottom:** sin ϕ_{i}. Parameter values: *N* = 1024, Δ = 0.02, *I*_{0} = −0.16, *J*_{0} = −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. 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

where

and

where |*i* − *j*| refers to the shortest distance between neurons *i* and *j*, measured on the ring. When *p*_{1} = *p*_{2} = *p*_{3} = 0, Equation (15) reverts to Equations (9) to (11). Note that when *p*_{1} = 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*_{1}. Similar statements apply for the other two matrices and their parameters *p*_{2} and *p*_{3}. Typical variation of *A*^{IE} with *p*_{1} is shown in Figure 3 and it is clear that increasing *p*_{1} interpolates between purely local connections (*p*_{1} = 0) and uniform random connectivity (*p*_{1} = 1).

**Figure 3. Typical realizations of A^{IE} for p_{1} = 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*_{1}, *p*_{2}, and *p*_{3} 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

where *G*_{IE}(*x, p*_{1}) 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 ${\int}_{-1/2}^{1/2}}{G}_{IE}(x,{p}_{1})\text{}dx=2{\alpha}_{IE$ independent of *p*_{1}, i.e., the expected total number of connections is preserved, and similarly for the other two functions.

**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)

and similarly *F*_{I} satisfies

where

and

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 *n*th coefficient is some function to the *n*th 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

and

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}

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

and

where

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*_{1}, *p*_{2} and *p*_{3} 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 4*N* ordinary differential equations (ODEs), 2*N* of them for complex quantities and the other 2*N* 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

and that of the inhibitory one by

which satisfy the continuity equations

and

where *dθ*_{j}/*dt* and *dϕ*_{j}/*dt* are given by Equations (3) and (4).

Performing the manipulations in the Appendix we obtain

for *j* = 1, 2, …*N* where

and

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*_{1}, *p*_{2} and *p*_{3} respectively we now investigate how solutions of Equations (42)–(48) depend on these parameters. One difficulty in trying to vary, say, *p*_{1}, is that the entries of *A*^{IE} do not depend continuously on *p*_{1}. Indeed, as presented, one should recalculate *A*^{IE} each time *p*_{1} 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:

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*_{1}, 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*_{1} is varied monotonically in the interval [0, 1].

The effects of quasistatically increasing *p*_{1} and *p*_{3} 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*_{1}, *p*_{2} and *p*_{3} independently, keeping the other two parameters fixed at zero. The results of varying *p*_{1} 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*_{1} 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.8191*i* and at the rightmost it has eigenvalues ±1.7972*i*, 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*_{1} is increased, such that when *p*_{1} = 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} = 1 all inhibitory neurons receive the same input from the excitatory population.

**Figure 5. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p_{1}, with p_{2} = p_{3} = 0**. There is a Hopf bifurcation on both white vertical lines and the bump is unstable between these. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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*_{2} while keeping *p*_{1} = *p*_{3} = 0 we find that the bump undergoes a Hopf bifurcation (Jacobian has eigenvalues ±0.3404*i*) and then is destroyed in a saddle-node bifurcation at *p*_{2} ≈ 0.48, as shown in Figure 6. The behavior of the bumps for 0 ≤ *p*_{2} ≤ 0.48 is shown in Figure 7.

**Figure 6. Maximum (over x) of the firing rate for the excitatory population as a function of p_{2} with p_{1} = p_{3} = 0**. Solid: stable; dashed: unstable. The Hopf bifurcation is marked with a circle. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p_{2}, with p_{3} = p_{1} = 0**. There is a Hopf bifurcation at the white vertical line and the bump is destroyed in saddle-node bifurcation at

*p*

_{2}≈ 0.48. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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*_{3} we obtain Figure 8, where there are no bifurcations as *p*_{3} 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*_{3} is increased.

**Figure 8. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p_{3}, with p_{2} = p_{1} = 0**. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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*_{1}, *p*_{2} and *p*_{3}, 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*_{1} = *p*_{2} = *p*_{3} = *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*_{1}, while keeping *p*_{2} = *p*_{3} = 0. We initially set *p*_{1} = 0 and integrated Equations (42)–(48) to a steady state, using initial conditions that give a bump solution. We then increased *p*_{1} 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} = 1. The firing rate for the *j*th 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*_{1} 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. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p_{1}, with p_{2} = p_{3} = 0**. Compare with Figure 5. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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*_{3} is varied with *p*_{1} = *p*_{2} = 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*_{3}, as in Figure 9. Varing *p*_{2} with *p*_{1} = *p*_{3} = 0 we obtain similar results to those in Figure 7 (not shown). Typical behavior of (42)–(48) with *p*_{2} = 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. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p_{3}, with p_{2} = p_{1} = 0**. Compare with Figure 8. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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. Instantaneous firing rate for (A): excitatory population and (B): inhibitory population, with p_{2} = 0.3 and p_{3} = p_{1} = 0**. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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*_{1} 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*_{3} we obtained Figure 13 (compare with Figure 10). Figure 14 shows oscillatory behavior at *p*_{2} = 0.3, *p*_{1} = *p*_{3} = 0 (the same parameter values as used in Figure 11).

**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_{1} with p_{2} = p_{3} = 0**. Compare with Figure 9. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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. 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_{3} with p_{2} = p_{1} = 0**. Compare with Figure 10. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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. Behavior of the full network Equations (3)–(6) with p_{2} = 0.3 and p_{3} = p_{1} = 0**.

**(A)**: 1 − cos θ

_{j},

**(B)**:

*v*

_{j}. Compare with Figure 11. Other parameters: Δ = 0.02,

*I*

_{0}= −0.16,

*J*

_{0}= −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*_{2}) 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## Appendix

### Mathematical Details Relating to Section 2.4

In the limit of an infinite ensemble we have

where ${f}_{j}^{E}({\theta}_{j},{I}_{j},t)$ is the marginal distribution for θ_{j}, given by

and similarly

Multiplying the continuity equation (40) by ∏_{k≠j} *dθ*_{k} *dI*_{k} and integrating we find that each ${f}_{j}^{E}$ satisfies

Similarly each ${f}_{j}^{I}$ satisfies

Using the Ott/Antonsen ansatz we write

and

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

and

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

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

Evaluating Equation (A10) at *I*_{j} = *I*_{0} + *i*Δ and Equation (A11) at *J*_{j} = *J*_{0} + *i*Δ we obtain

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

Keywords: Ott/Antonsen, theta neuron, bump, small-world, working memory, bifurcation

Citation: Laing CR (2016) Bumps in Small-World Networks. *Front. Comput. Neurosci.* 10:53. doi: 10.3389/fncom.2016.00053

Received: 25 February 2016; Accepted: 23 May 2016;

Published: 15 June 2016.

Edited by:

Tomoki Fukai, RIKEN Brain Science Institute, JapanReviewed by:

Robert Kozma, University of Memphis, USAVladimir Klinshov, Institute of Applied Physics, Russia

Copyright © 2016 Laing. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Carlo R. Laing, c.r.laing@massey.ac.nz