Impact Factor 2.073

The world's most-cited Neurosciences journals

Original Research ARTICLE

Front. Comput. Neurosci., 15 June 2016 | https://doi.org/10.3389/fncom.2016.00053

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

dθidt=1-cosθi+(1+cosθi)(Ii+gr)    (1)
τdrdt=1Nj=1NPn(θj)-r    (2)

for i = 1, 2, …N where θi ∈ [0, 2π) is the phase of the ith neuron, Pn(θ)=an(1-cosθ)n,n+ and an is a normalization factor such that 02πPn(θ)dθ=2π. The function Pn 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 Ii 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

dθidt=1-cosθi+(1+cosθi)(Ii+gEEvi-gEIyi)    (3)
dϕidt=1-cosϕi+(1+cosϕi)(Ji+gIEui-gIIzi)    (4)
τdvidt=ri-vi    (5)
τduidt=qi-ui    (6)
τidyidt=si-yi    (7)
τidzidt=wi-zi    (8)

for i = 1, 2…N, where

qi=1Nj=-MIEMIEPn(θi+j)    (9)
ri=1Nj=-MEEMEEPn(θi+j)    (10)
si=1Nj=-MEIMEIPn(ϕi+j)    (11)
wi=1Nj=-MIIMIIPn(ϕi+j)    (12)

where Pn is as in Section 2.1. The positive integers MIE, MEE, MEI, and MII 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 gEE, gEI, gIE and gII 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 vi (when multiplied by gEE) gives the excitatory input to the ith excitatory neuron, and whose dynamics are driven by ri, which depends on the activity of the excitatory neurons with indices between iMEE and i + MEE. Similarly, ui (when multiplied by gIE) gives the excitatory input to the ith inhibitory neuron, and is driven by qi, which depends on the activity of the excitatory neurons with indices between iMIE and i + MIE. gEIyi is the inhibitory input to the ith excitatory neuron, driven by si, which depends on the activity of the inhibitory neurons with indices between iMEI and i + MEI. Lastly, gIIzi is the inhibitory input to the ith inhibitory neuron, driven by wi, which depends on the activity of the inhibitory neurons with indices between iMII and i + MII.

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. gII = 0. Thus (8) and (12) become irrelevant and from (7) we have that yi = si in (3).

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

h(I)=Δ/π(I-I0)2+Δ2    (13)

and the currents Ji from the Lorentzian

g(J)=Δ/π(J-J0)2+Δ2.    (14)

I0 and J0 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
www.frontiersin.org

Figure 1. A bump solution of Equations (3)–(6). Top: sin θi. Bottom: sin ϕi. Parameter values: N = 1024, Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, MIE = 40, MEE = 40, MEI = 60 and τ = 10.

FIGURE 2
www.frontiersin.org

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

qi=1Nj=1NAijIEPn(θj);ri=1Nj=1NAijEEPn(θj);                                            si=1Nj=1NAijEIPn(ϕj);    (15)

where

AijIE={1withprobability{1[1(2MIE+1)/N]p1,|ij|MIE(2MIE+1)p1/N,|ij|>MIE0otherwise    (16)
AijEE={1withprobability{1[1(2MEE+1)/N]p2,|ij|MEE(2MEE+1)p2/N,|ij|>MEE0 otherwise    (17)

and

AijEI={1withprobability{1[1(2MEI+1)/N]p3,|ij|MEI(2MEI+1)p3/N,|ij|>MEI0 otherwise    (18)

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

FIGURE 3
www.frontiersin.org

Figure 3. Typical realizations of AIE for p1 = 0 (left) 0.5 (middle) and 1 (right). N = 1024, MIE = 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 p1, p2, and p3 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, MEI, MEE, MIE → ∞ such that MEI/N → αEI, MEE/N → αEE and MIE/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

GIE(x,p1)={1(12αIE)p1,|x|<αIE2αIEp1,otherwise    (19)
GEE(x,p2)={1(12αEE)p2,|x|<αEE2αEEp2,otherwise    (20)
GEI(x,p3)={1(12αEI)p3,|x|<αEI2αEIp3,otherwise    (21)

where GIE(x, p1) 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 pj, j = 1, 2, 3, on one of the functions Equations (19)–(21) is shown in Figure 4. Taking GIE for example, we see that -1/21/2GIE(x,p1) dx=2αIE independent of p1, i.e., the expected total number of connections is preserved, and similarly for the other two functions.

FIGURE 4
www.frontiersin.org

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 FE(θ, x, I, t) and FI(ϕ, 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. FE satisfies the continuity equation (Luke et al., 2013)

FEt+θ{FE[1-cosθ+(1+cosθ)(I+gEEv-gEIs)]}=0    (22)

and similarly FI satisfies

FIt+ϕ{FI[1-cosϕ+(1+cosϕ)(J+gIEu)]}=0    (23)

where

τvt=r-v    (24)
τut=q-u    (25)

and

q(x,t)=01GIE(|x-y|,p1)-02π                FE(θ,y,I,t)an(1-cosθ)ndθ dI dy    (26)
r(x,t)=01GEE(|x-y|,p2)-02π                FE(θ,y,I,t)an(1-cosθ)ndθ dI dy    (27)
s(x,t)=01GEI(|x-y|,p3)-02π                FI(ϕ,y,J,t)an(1-cosϕ)ndϕ dJ dy    (28)

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 FE and FI, 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 eiθj (Strogatz, 2000). For the network studied here we can define the analogous spatially-dependent order parameters for the excitatory and inhibitory networks as

zE(x,t)=-02πFE(θ,x,I,t)eiθ dθ dI    (29)

and

zI(x,t)=-02πFI(ϕ,x,J,t)eiϕ dϕ dJ    (30)

respectively. For fixed x and t, zE(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 zI(x, t) and ϕ (Laing, 2014a, 2015). We can also determine from zE and zI 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 zE and zI

zEt=(iI0-Δ)(1+zE)2-i(1-zE)22             +i(1+zE)2(gEEv-gEIs)2    (31)
zIt=(iJ0-Δ)(1+zI)2-i(1-zI)22+i(1+zI)2gIEu2    (32)

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

q(x,t)=01GIE(|x-y|,p1)H(zE(y,t);n)dy    (33)
r(x,t)=01GEE(|x-y|,p2)H(zE(y,t);n)dy    (34)
s(x,t)=01GEI(|x-y|,p3)H(zI(y,t);n)dy    (35)

and

H(z;n)=an[C0+q=1nCq(zq+z¯q)]    (36)

where

Cq=k=0nm=0kn!(-1)kδk-2m,q2k(n-k)!m!(k-m)!    (37)

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

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 p1, p2 and p3 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 AIE, AEE and AEI, 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 Ii and Ji (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

fE(θ1,θ2,,θN;I1,I2,IN;t)fE({θ};{I};t)    (38)

and that of the inhibitory one by

fI(ϕ1,ϕ2,,ϕN;J1,J2,JN;t)fI({ϕ};{J};t)    (39)

which satisfy the continuity equations

fEt+j=1Nθj[fE(dθjdt)]=0    (40)

and

fIt+j=1Nϕj[fI(dϕjdt)]=0    (41)

where j/dt and j/dt are given by Equations (3) and (4).

Performing the manipulations in the Appendix we obtain

dzjEdt=(iI0-Δ)(1+zjE)2-i(1-zjE)22          +i(1+zjE)2(gEEvj-gEIsj)2    (42)
dzjIdt=(iJ0-Δ)(1+zjI)2-i(1-zjI)22+i(1+zjI)2gIEuj2    (43)

for j = 1, 2, …N where

qi=1Nj=1NAijIEH(zjE(t);n)    (44)
ri=1Nj=1NAijEEH(zjE(t);n)    (45)
si=1Nj=1NAijEIH(zjI(t);n)    (46)

and

τdvidt=ri-vi    (47)
τduidt=qi-ui    (48)

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 AIE, AEE and AEI. 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 AIE, AEE and AEI depend on the parameters p1, p2 and p3 respectively we now investigate how solutions of Equations (42)–(48) depend on these parameters. One difficulty in trying to vary, say, p1, is that the entries of AIE do not depend continuously on p1. Indeed, as presented, one should recalculate AIE each time p1 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 AIE (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 AIE(p1) as follows:

AijIE(p1)={Θ[rijp1(1(2MIE+1)/N)],|ij|MIEΘ[rij(1p1(2MIE+1)/N)],|ij|>MIE    (49)

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

The effects of quasistatically increasing p1 and p3 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 pi 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 p1, p2 and p3 independently, keeping the other two parameters fixed at zero. The results of varying p1 are shown in Figure 5, where we plot the firing rate of the two populations, derived as Re(wi)/π where wi=(1-z¯i)/(1+z¯i) for i = I, E, as in Montbrió et al. (2015), where the zi are fixed points of Equations (31) and (32). We see an increase and then decrease in bump width as p1 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 p1 is increased, such that when p1 = 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 p1 = 1 all inhibitory neurons receive the same input from the excitatory population.

FIGURE 5
www.frontiersin.org

Figure 5. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p1, with p2 = p3 = 0. There is a Hopf bifurcation on both white vertical lines and the bump is unstable between these. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, αIE = 40/1024, αEE = 40/1024, αEI = 60/1024 and τ = 10.

Increasing p2 while keeping p1 = p3 = 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 p2 ≈ 0.48, as shown in Figure 6. The behavior of the bumps for 0 ≤ p2 ≤ 0.48 is shown in Figure 7.

FIGURE 6
www.frontiersin.org

Figure 6. Maximum (over x) of the firing rate for the excitatory population as a function of p2 with p1 = p3 = 0. Solid: stable; dashed: unstable. The Hopf bifurcation is marked with a circle. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, αIE = 40/1024, αEE = 40/1024, αEI = 60/1024 and τ = 10.

FIGURE 7
www.frontiersin.org

Figure 7. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p2, with p3 = p1 = 0. There is a Hopf bifurcation at the white vertical line and the bump is destroyed in saddle-node bifurcation at p2 ≈ 0.48. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, αIE = 40/1024, αEE = 40/1024, αEI = 60/1024 and τ = 10.

Varying p3 we obtain Figure 8, where there are no bifurcations as p3 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 p3 is increased.

FIGURE 8
www.frontiersin.org

Figure 8. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p3, with p2 = p1 = 0. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 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 p1, p2 and p3, 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 p1 = p2 = p3 = 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 p1, while keeping p2 = p3 = 0. We initially set p1 = 0 and integrated Equations (42)–(48) to a steady state, using initial conditions that give a bump solution. We then increased p1 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 p1 = 1. The firing rate for the jth excitatory neuron is Re(wj)/π where wj=(1-z¯jE)/(1+z¯jE), 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 p1 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 AIE. 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
www.frontiersin.org

Figure 9. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p1, with p2 = p3 = 0. Compare with Figure 5. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, N = 1024, MIE = 40, MEE = 40, MEI = 60 and τ = 10.

Repeating this process as p3 is varied with p1 = p2 = 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 p3, as in Figure 9. Varing p2 with p1 = p3 = 0 we obtain similar results to those in Figure 7 (not shown). Typical behavior of (42)–(48) with p2 = 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
www.frontiersin.org

Figure 10. Firing rate for (A): excitatory population and (B): inhibitory population, as a function of p3, with p2 = p1 = 0. Compare with Figure 8. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, N = 1024, MIE = 40, MEE = 40, MEI = 60 and τ = 10.

FIGURE 11
www.frontiersin.org

Figure 11. Instantaneous firing rate for (A): excitatory population and (B): inhibitory population, with p2 = 0.3 and p3 = p1 = 0. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, N = 1024, MIE = 40, MEE = 40, MEI = 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 AEE and AEI) to calculate Equation (15). The frequency was measured directly from simulations. Varying p1 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 p3 we obtained Figure 13 (compare with Figure 10). Figure 14 shows oscillatory behavior at p2 = 0.3, p1 = p3 = 0 (the same parameter values as used in Figure 11).

FIGURE 12
www.frontiersin.org

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 p1 with p2 = p3 = 0. Compare with Figure 9. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, N = 1024, MIE = 40, MEE = 40, MEI = 60 and τ = 10.

FIGURE 13
www.frontiersin.org

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 p3 with p2 = p1 = 0. Compare with Figure 10. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, N = 1024, MIE = 40, MEE = 40, MEI = 60 and τ = 10.

FIGURE 14
www.frontiersin.org

Figure 14. Behavior of the full network Equations (3)–(6) with p2 = 0.3 and p3 = p1 = 0. (A): 1 − cos θj, (B): vj. Compare with Figure 11. Other parameters: Δ = 0.02, I0 = −0.16, J0 = −0.4, n = 2, gEE = 25, gIE = 25, gEI = 7.5, N = 1024, MIE = 40, MEE = 40, MEI = 60 and τ = 10.

Note that the results in Figures 914 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 p2) 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

qi=1Nj=1NAijIE···Pn(θj)fE({θ};{I};t) dθ1 dθ2dθN dI1 dI2dIN    =1Nj=1NAijIE-02πPn(θj)fjE(θj,Ij,t)dθj dIj    (A1)
ri=1Nj=1NAijEE···Pn(θj)fE({θ};{I};t) dθ1 dθ2dθN dI1 dI2dIN   =1Nj=1NAijEE-02πPn(θj)fjE(θj,Ij,t)dθj dIj    (A2)
si=1Nj=1NAijEI···Pn(ϕj)fI({ϕ};{J};t) dϕ1 dϕ2dϕN dJ1 dJ2dJN   =1Nj=1NAijEI-02πPn(ϕj)fjI(ϕj,Jj,t)dϕj dJj    (A3)

where fjE(θj,Ij,t) is the marginal distribution for θj, given by

fjE(θj,Ij,t)=···fE({θ};{I};t)kjdθk dIk    (A4)

and similarly

fjI(ϕj,Jj,t)=···fI({ϕ};{J};t)kjdϕk dJk    (A5)

Multiplying the continuity equation (40) by ∏kj k dIk and integrating we find that each fjE satisfies

fjEt+θj[fjE(dθjdt)]=0    (A6)

Similarly each fjI satisfies

fjIt+ϕj[fjI(dϕjdt)]=0    (A7)

Using the Ott/Antonsen ansatz we write

fjE(θj,Ij,t)=h(Ij)2π{1+n=1[αjE(Ij,t)]neinθj+c.c.}    (A8)

and

fjI(ϕj,Jj,t)=g(Jj)2π{1+n=1[αjI(Jj,t)]neinϕj+c.c.}    (A9)

for some functions αjE(Ij,t) and αjI(Jj,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

αjEt=-i[Ij+gEEvj-gEIsj-12+(Ij+gEEvj-gEIsj+1)αjE+Ij+gEEvj-gEIsj-12(αjE)2]    (A10)

and

αjIt=-i[Jj+gIEuj-12+(Jj+gIEuj+1)αjI+Jj+gIEuj-12(αjI)2]    (A11)

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

qi=1Nj=1NAijIE-h(Ij)H(αjE(Ij,t);n) dIj    (A12)
ri=1Nj=1NAijEE-h(Ij)H(αjE(Ij,t);n) dIj    (A13)
si=1Nj=1NAijEI-g(Jj)H(αjI(Jj,t);n) dJj    (A14)

where H is given by Equation (36). Using standard properties of the Lorentzian one can perform the integrals in Equations (A12)–(A14) and defining zjE(t)α¯jE(I0+iΔ,t) and zjI(t)α¯jI(J0+iΔ,t) we have

qi=1Nj=1NAijIEH(zjE(t);n)    (A15)
ri=1Nj=1NAijEEH(zjE(t);n)    (A16)
si=1Nj=1NAijEIH(zjI(t);n)    (A17)

Evaluating Equation (A10) at Ij = I0 + iΔ and Equation (A11) at Jj = J0 + iΔ we obtain

dzjEdt=(iI0-Δ)(1+zjE)2-i(1-zjE)22+i(1+zjE)2(gEEvj-gEIsj)2    (A18)
dzjIdt=(iJ0-Δ)(1+zjI)2-i(1-zjI)22+i(1+zjI)2gIEuj2    (A19)

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, Japan

Reviewed by:

Robert Kozma, University of Memphis, USA
Vladimir 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