Patterns of synchronization in 2D networks of inhibitory neurons

1. Introduction

In different brain regions, neural firing patterns can be characterized as coordinated assembly activity, in which individual collections of cells fire synchronously and not at the same time as other collections of synchronously firing cells (Gray et al., 1989; Laurent and Davidowitz, 1994; Harris et al., 2003; Dragoi and Buzsáki, 2006; Muldoon et al., 2013; Barbera et al., 2016). Neural assemblies may consist of spatially localized neurons (Muldoon et al., 2013; Barbera et al., 2016) or of neurons that are widespread across one or more brain regions (Gray et al., 1989; Harris et al., 2003). For example, neural assemblies have been observed between neurons in different cortical columns (Gray et al., 1989), within regions of the hippocampus (Harris et al., 2003; Dragoi and Buzsáki, 2006), the dentate gyrus (Muldoon et al., 2013), and between cells in the striatum (Carrillo-Reid et al., 2008; Adler et al., 2013; Barbera et al., 2016) and the olfactory bulb (Laurent and Davidowitz, 1994). Understanding the dynamics and formation of synchronized assemblies within larger neural networks has gained increasing importance in neuroscience (Engel et al., 2001) and has been studied both experimentally (Laurent and Davidowitz, 1994; Harris et al., 2003; Galán et al., 2006; Muldoon et al., 2013; Barbera et al., 2016) and using computational modeling (Golomb et al., 1992; Golomb and Rinzel, 1994; Li et al., 2003; Achuthan and Canavier, 2009; Canavier et al., 2009; Ponzi and Wickens, 2010; Kilpatrick and Ermentrout, 2011; Angulo-Garcia et al., 2016).

One approach to understanding the formation of neural assemblies is through the study of cluster solutions in model networks of intrinsically oscillating neurons (Golomb et al., 1992; Golomb and Rinzel, 1994; Li et al., 2003; Galán et al., 2006; Kilpatrick and Ermentrout, 2011; Miller et al., 2015; Campbell and Wang, 2018; Ryu et al., 2021). Cluster solutions are solutions where the network of oscillators breaks into subgroups. Within each subgroup, the phases of the oscillators are the same, while oscillators in different subgroups are phase-locked with some non-zero phase difference. A useful mathematical framework for studying cluster solutions is the phase model reduction (Hoppensteadt and Izhikevich, 1997; Schwemmer and Lewis, 2012). This framework has been used to study synchronization and clustering in a variety of coupled oscillator networks (Ashwin and Swift, 1992; Okuda, 1993; Kopell and Ermentrout, 2002; Saraga et al., 2006; Mancilla et al., 2007). A complementary approach is to consider a continuum model representing the limit of an infinite number of oscillators (Ermentrout, 1992; Strogatz, 2000). In such models, which are represented as partial differential equations, cluster solutions correspond to wave-like solutions, sometimes called twisted states (Ermentrout, 1992; Wiley et al., 2006; Kazanci and Ermentrout, 2007; Girnyk et al., 2012; Heitmann et al., 2012; Heitmann and Ermentrout, 2015).

Cluster solutions have been extensively studied in models that assume homogeneous, all-to-all connectivity among oscillators (Sakaguchi and Kuramoto, 1986; Ashwin and Swift, 1992; Hansel et al., 1993; Okuda, 1993; Strogatz, 2000; Zanette, 2000; Kazanci and Ermentrout, 2008; Campbell and Wang, 2018). The number of studies considering cluster solutions in networks with structured connectivity is more limited and they primarily focus on networks with neurons arranged in a 1D ring, of arbitrary size N. In Kazanci and Ermentrout (2007, 2008), the phase model reduction is used to analyze how the addition of local (nearest neighbor) gap-junctional coupling in a network with all-to-all synaptic coupling affects cluster solutions. They found that gap-junctional coupling can induce a shift from synchronous or two cluster solutions to an N-cluster solution. The presence of time-delayed synaptic connections among all-to-all coupled oscillators was analyzed using a more general phase model in Campbell and Wang (2018). They found that time delays can lead to the coexistence of multiple stable clustering solutions. In a continuum network model, the existence and stability of cluster solutions when each oscillator has identical coupling to a subset of its nearest neighbors was analyzed in Wiley et al. (2006); Girnyk et al. (2012). In a study considering structured excitatory-inhibitory coupling in a continuum model (Heitmann and Ermentrout, 2015), the extent of inhibitory coupling was shown to influence the existence of synchronous or traveling wave solutions in both 1D and 2D networks. Recently, a phase model analysis of cluster solutions was extended to 2D networks in Culp (2021), in which the existence of cluster solutions and stability conditions of a particular type of clusters (that have identical phase difference between adjacent oscillators in the horizontal and vertical directions) were derived on a square torus.

In our previous work (Miller et al., 2015; Ryu et al., 2021), we used the phase model approach to determine existence and stability conditions for cluster solutions in 1D, discrete inhibitory neural networks with various connectivity schemes. In Miller et al. (2015) we investigated the clustering dynamics of a network of oscillating inhibitory neurons in which each neuron is coupled to its two nearest neighbors on each side. Our phase model analysis showed that changing the connection weights can change the stability of solutions, confirmed by numerical simulations. In Ryu et al. (2021), we considered simple non-monotonic, distance-dependent connectivity schemes in 1D inhibitory networks. We similarly used the phase model approach to identify conditions for the existence and stability of cluster firing solutions in which clusters consist of spatially adjacent neurons in inhibitory neural networks.

In this paper, we extend our previous work to a 2D, discrete inhibitory network model on a general torus, where neurons are coupled to two or more nearest neighbors in three directions (horizontal, vertical, and diagonal) with different coupling strengths in all directions. Using phase model analysis, we derive stability conditions for different patterns of clustered firing behavior in the networks, and further explore the effect of heterogeneous coupling strengths on the stability of cluster solutions. Stability conditions are derived for a general phase model, but applied to the specific phase model derived from a biophysical network model consisting of inhibitory interneurons (Wang and Buzsáki, 1996). We show that the phase model predictions give an accurate picture of the effect of various connectivity schemes on the solutions that occur in the full biophysical model.

Our paper is structured as follows: Section 2 provides the description of the methodology we employ, including the reduction from a biophysical model to a phase model (Section 2.1) and the specific biophysical neuron model used for numerical simulations (Section 2.2). Section 3 describes our analysis of the existence and stability of cluster solutions in a general phase model. Section 4 describes the application of these results to the specific phase model corresponding to the biophysical network and the verification of these results via numerical simulations. We conclude with a discussion of our results in Section 5.

2. Methods

2.1. Phase reduction method

In this work, we consider a general neural network model on a 2D lattice that consists of m×n identical, weakly coupled, inhibitory oscillating cells with periodic boundary conditions:

$\begin{array}{l}\frac{d{X}_{ij}}{dt}=F(X_{ij})+ϵ{\displaystyle \sum _{p,q}}{w}_{p,q}G(X_{ij},X_{i+p,j+q}),\\ 0\le i\le n-1,0\le j\le m-1.& (1)\end{array}$

Here, each X_ij defines a k-dimensional variable. F:ℝ^k → ℝ^k is the internal vector field of the isolated oscillator and G:ℝ^k × ℝ^k → ℝ^k is the synaptic coupling function. W = (w_{p, q}) is the connection matrix with w_{p, q} > 0 when a synapse exists from cell q to cell p. We apply the phase reduction method to the above model to study (i) the existence and stability of certain cluster solutions on a 2D general torus network, and (ii) how the coupling structure on the 2D networks can affect the stability of these cluster solutions.

In the remaining part of this subsection, we will review the phase reduction method. Interested readers can refer to Hoppensteadt and Izhikevich (1997); Ermentrout and Terman (2010); Schwemmer and Lewis (2012) for more details.

Each uncoupled single neuron is assumed to admit an exponentially asymptotically stable T-periodic orbit, $\left\{\widehat{X}(t):0\le t\le T=2\pi /\Omega \right\}$, and $\widehat{X}(t)$ satisfies

$\begin{array}{ll}\frac{dX}{dt}=F(X(t)),X\in {ℝ}^k.& (2)\end{array}$

By the theory of weakly coupled oscillators (see Hoppensteadt and Izhikevich, 1997; Ermentrout and Terman, 2010; Schwemmer and Lewis, 2012), the complete state of each neuron in the network can be approximately captured by its phase on the standardized limit cycle θ(t), for 0 ≤ t < 2π. This decreases the number of equations that describe an uncoupled single neuron from k equations to one, which in turn can significantly reduce the dimension of the neural network model. The phase of the ij^th neuronal oscillator is slowly varying and its dynamics is governed by

$\begin{array}{ll}\frac{d{\theta }_{ij}}{dt}=\Omega +ϵ{\displaystyle \sum _{p,q}}{w}_{p,q}H({\theta }_{i+p,j+q}-{\theta }_{ij}),& (3)\end{array}$

where

$\begin{array}{ll}H(\psi )=\frac{1}{T}{\displaystyle {\int }_0^T}Z(t)G[\widehat{X}(t),\widehat{X}(t+\psi /\Omega )]dt.& (4)\end{array}$

The function H, often called the interaction function, measures the modulation of the instantaneous phase of the i^th oscillator due to its coupling with other oscillators on the network. The function Z is the unique periodic solution of the linearized adjoint system

$\begin{array}{l}\frac{dZ}{dt}=-{\left[DF(\widehat{X}(t))\right]}^{T}Z,\end{array}$

satisfying the normalization condition

$\begin{array}{l}\frac{1}{T}{\displaystyle {\int }_0^T}Z(t)·F(\widehat{X}(t))dt=1.\end{array}$

2.2. Neural network model

We apply our theoretical results to a biophysical inhibitory neural network model, where each individual neuron is modeled by the conductance-based Wang and Buzsáki inhibitory interneuron model (Wang and Buzsáki, 1996). Networks consist of N neurons coupled in a 2D lattice structure through inhibitory synaptic currents with N = m×n. Membrane voltage of the ith cell, V_i (in mV) is modeled by:

$\begin{array}{cc}\begin{array}{l}C\frac{d{V}_i}{dt}=I_{app}\\ -g_{Na}{m}_{\infty }^3(V_i)h_i(V_i-V_{Na})-g_K{n}_i^4(V_i-V_K)\\ -g_L(V_i-V_L)\\ -g_{syn}(V_i-V_{syn}){\displaystyle \sum _{j=1}^N}{W}_{ij}{s}_j,\\ =I_{app}-I_{ion}(V_i,h_i,n_i)-I_{syn}(\overrightarrow{V},\overrightarrow{s})=f^V(\overrightarrow{V},h_i,n_i,\overrightarrow{s}),\\ \frac{d{h}_i}{dt}=\gamma ({\alpha }_h(V_i)(1-h_i)-{\beta }_h(V_i)h_i)=f^h(V_i,h_i),\\ \frac{d{n}_i}{dt}=\gamma ({\alpha }_n(V_i)(1-n_i)-{\beta }_n(V_i)n_i)=f^n(V_i,n_i),\\ \frac{d{s}_i}{dt}=-\frac{s_i}{{\tau }_{inh}}+{\alpha }_{inh}(V_i)(1-s_i),\end{array}& (5)\end{array}$

with

$\begin{array}{ll}\begin{array}{l}{m}_{\infty }^3(V)=\frac{-0.1(V+35)/(exp(-0.1(V+35))-1)}{4exp(-(V+60)/18)},\\ {\alpha }_h(V)=0.07exp(-(V+58)/20),\\ {\beta }_h(V)=1/(exp(-0.1(V+28))+1),\\ {\alpha }_n(V)=-0.01(V+34)/(exp(-0.1(V+34))-1),\\ {\beta }_n(V)=0.125exp(-(V+44)/80),\\ {\alpha }_{inh}(V)={\alpha }_0/(1+exp(-V/5)),\end{array}& (6)\end{array}$

where h_i is the gating variable governing inactivation of the inward sodium current, n_i is the gating variable governing activation of the outward potassium current and s_i is the gating variable for the synaptic current generated by presynaptic cell i. Coupling between neurons is dictated by the connectivity matrix W_ij and maximum synaptic strength is set to g_syn. The definition of model parameters and their default values are provided in Table 1.

TABLE 1

Table 1. Description of parameters and the values used for the neuron model.

3. Analysis results

We consider the network on a perfect m×n lattice. We can define coupling in the network by letting S={(1,0), (-1,0), (0, 1), (0, -1), (-1, 1), (1,1), (1,-1), (-1, -1), (2, 0), (-2, 0), (0, 2), (0, -2)}. A network of phase oscillators with S-coupling means that each neuron is connected to its 8 nearest neighbors (vertical, horizontal, and diagonal directions) and 4 more distant neighbors (second nearest neighbors in the vertical and horizontal directions). An example showing the 12 neighbors of one neuron is shown in Figure 1.

FIGURE 1

Figure 1. The connections to the 12 neighbors of the neuron i, j are shown by solid edges. Connections to nearest neighbors in the vertical and horizontal directions are shown in red. Connections to diagonal neighbors are blue. Connections to second nearest neighbors vertically and horizontally are green.

The network model takes the form of (3). In the above system of equations and what follows, we assume that all the horizontal (or vertical) indices are defined under the mod operation mod n (or m).

Phase-locked solutions of a phase model such as system (3), will be solutions where all oscillators evolve with the same frequency. Thus, we look for solutions of (3) in the form

$\begin{array}{ll}{\bar{\theta }}_{ij}=(\Omega +ϵ\omega )t+{\varphi }_{ij}.& (7)\end{array}$

Substituting (7) into (3) we find

$\begin{array}{ll}\omega ={\displaystyle \sum _{(p,q)\in S}}{w}_{p,q}H{(\varphi }_{(i+p),(j+q)}-{\varphi }_{ij}).& (8)\end{array}$

Since ω must be constant and (8) should be independent of all i, j, we consider a special type of solution with constant phase difference along the horizontal and vertical directions. Let ψ_h and ψ_v be the phase differences between adjacent cells in the horizontal and vertical directions, respectively. For any p_h that divides n, it follows from the periodic boundary condition that a p_h-cluster solution along the horizontal direction implies that p_hψ_h = 2ℓ_hπ for some ℓ_h < p_h and gcd(ℓ_h, p_h) = 1. Similarly, in the vertical direction, for any p_v that divides m, a p_v-cluster along the vertical direction implies that p_vψ_v = 2ℓ_vπ for some ℓ_v < p_v and gcd(ℓ_v, p_v) = 1. For the network as a whole this gives a p-cluster solution with p = lcm(p_h, p_v) (i.e., p is the least common multiple of p_h and p_v). Hence, we have the following existence result.

Proposition 1. Suppose that p_h is a factor of n, and p_v is a factor of m. Let p = lcm(p_h, p_v). Then a torus network of m×n neurons with S-coupling admits a p-cluster solution with constant phase differences ψ_h and ψ_v between adjacent cells in the horizontal and vertical directions, respectively. These phase differences satisfy ${\psi }_h=\frac{2{\ell }_h\pi }{p_h}$ and ${\psi }_v=\frac{2{\ell }_v\pi }{p_v}$ such that ℓ_h < p_h and gcd(ℓ_h, p_h) = 1, and ℓ_v < p_v and gcd(ℓ_v, p_v) = 1.

We define the solutions described in Proposition 1 as (p_h, p_v) p-cluster solutions. For such a p-cluster solution with phase differences ψ_h and ψ_v, we have that

$\begin{array}{l}\omega =w_{1,0}H({\psi }_h)+w_{0,1}H({\psi }_v)\\ +w_{-1,0}H(-{\psi }_h)+w_{0,-1}H(-{\psi }_v)+w_{-1,1}H(-{\psi }_h+{\psi }_v)\\ +w_{1,-1}H({\psi }_h-{\psi }_v)\\ +w_{-1,-1}H(-{\psi }_h-{\psi }_v)+w_{1,1}H({\psi }_h+{\psi }_v)\\ +w_{2,0}H(2{\psi }_h)+w_{0,2}H({\psi }_v)\\ +w_{-2,0}H(-2{\psi }_h)+w_{0,-2}H(-2{\psi }_v)\\ ={\displaystyle \sum _{(p,q)\in S}}{w}_{p,q}H(p{\psi }_h+q{\psi }_v)\end{array}$

is constant.

We now determine the stability of these (p_h, p_v) p-cluster solutions. Let ${\theta }_{ij}={\bar{\theta }}_{ij}+u_{ij}(t)$. Linearizing (3) about ${\bar{\theta }}_{ij}$ leads to

$\begin{array}{l}\frac{d{u}_{ij}}{dt}\\ =ϵ{\displaystyle \sum _{p,q\in S}{w}_{p,q}}[-H^{\prime }({\theta }_{(i+p),(j+q)}-{\theta }_{ij})u_{ij}\\ +H^{\prime }({\theta }_{(i+p),(j+q)}-{\theta }_{ij})u_{(i+p),(j+q)}].\end{array}$

Rewriting the above equation in the matrix form, we have:

$\begin{array}{ll}\frac{d\overrightarrow{u}}{dt}=ϵ(-c{I}_N+Ŵ)\overrightarrow{u},& (9)\end{array}$

where

$\begin{array}{l}\overrightarrow{u}=(u_{0,0},u_{1,0},\dots ,u_{n-1,0},u_{0,1},u_{1,1},\dots ,u_{n-1,1},\dots ,\\ u_{0,m-1},u_{1,m-1},\dots ,u_{n-1,m-1}{)}^T,\\ {\widehat{w}}_{p,q}=w_{p,q}{H}^{\prime }(p{\psi }_h+q{\psi }_v),n\\ c={\displaystyle \sum _{(p,q)\in S}{\widehat{w}}_{p,q}},\end{array}$

and Ŵ = bcirc(Ŵ⁰, Ŵ¹, …, Ŵ^m−1) where

$\begin{array}{l}{\widehat{W}}^0=\text{circ}(0,{\widehat{w}}_{1,0},{\widehat{w}}_{2,0},0,\dots ,0,{\widehat{w}}_{-2,0},{\widehat{w}}_{-1,0}),\\ {\widehat{W}}^1=\text{circ}({\widehat{w}}_{0,1},{\widehat{w}}_{1,1},0,\dots ,0,{\widehat{w}}_{-1,1}),{\widehat{W}}^2={\widehat{w}}_{0,2}I,\\ {\widehat{W}}^{m-2}={\widehat{w}}_{0,-2}I,\\ {\widehat{W}}^{m-1}=\text{circ}({\widehat{w}}_{0,-1},{\widehat{w}}_{1,-1},0,\dots ,0,{\widehat{w}}_{-1,-1}),\end{array}$

and all other Ŵ^k are n×n zero matrices.

Let ${\rho }_n=e^{2\pi \sqrt{-1}/n}$. The eigenvalues of the non-zero Ŵ^k are

$\begin{array}{l}{\lambda }_j^{{\widehat{W}}^0}={\widehat{w}}_{1,0}{\rho }_n^j+{\widehat{w}}_{-1,0}{\rho }_n^{-j},\\ {\lambda }_j^{{\widehat{W}}^1}={\widehat{w}}_{0,1}+{\widehat{w}}_{1,1}{\rho }_n^j+{\widehat{w}}_{-1,1}{\rho }_n^{-j},\\ {\lambda }_j^{{\widehat{W}}^2}={\widehat{w}}_{0,2},\\ {\lambda }_j^{{\widehat{W}}^{m-2}}={\widehat{w}}_{0,-2},\\ {\lambda }_j^{{\widehat{W}}^{m-1}}={\widehat{w}}_{0,-1}+{\widehat{w}}_{1,-1}{\rho }_n^j+{\widehat{w}}_{-1,-1}{\rho }_n^{-j},& (10)\end{array}$

for 0 ≤ j ≤ n−1. Let

${\Lambda }_{\ell }=\text{diag}\left[{\lambda }_0^{{\widehat{W}}^{\ell }},{\lambda }_1^{{\widehat{W}}^{\ell }},\cdots ,{\lambda }_{n-1}^{{\widehat{W}}^{\ell }}\right],\ell =0,1,m-1.$

The matrix Ŵ is a block circulant matrix, and each block Ŵ^j, 0 ≤ j ≤ m−1, is circulant. Let $P_m=\text{diag}(1,{\rho }_m,{\rho }_m^2,\dots ,{\rho }_m^{m-1})$ and ⊗ denote the Kronecker product of matrices. By Davis (2013), the eigenvalues of Ŵ are given by the diagonal entries of

$\begin{array}{l}I\text{ }\otimes \text{ }{\Lambda }_0+P_m\otimes {\Lambda }_1+P_m^2\otimes {\Lambda }_2+P_m^{m-2}\otimes {\Lambda }_{m-2}\\ +P_m^{m-1}\otimes {\Lambda }_{m-1}\\ =\text{bdiag}\left[D_0,D_1,D_2,0_m,\cdots ,0_m,D_{m-2},D_{m-1}\right],\end{array}$

where

$D_i={\displaystyle \sum _{\ell \in L}}{\rho }_m^{\ell }{\Lambda }_{\ell }$

and L = {0, 1, 2, m−2, m−1}.

Thus, the eigenvalues of Ŵ are

$\begin{array}{l}{\lambda }_{jk}^{\widehat{W}}={\lambda }_j^{{\widehat{W}}^0}+{\rho }_m^k{\lambda }_j^{{\widehat{W}}^1}\\ +{\rho }_m^{2k}{\lambda }_j^{{\widehat{W}}^2}+{\rho }_m^{-2k}{\lambda }_j^{{\widehat{W}}^{m-2}}+{\rho }_m^{-k}{\lambda }_j^{{\widehat{W}}^{m-1}},\\ 0\le j\le n-1,0\le k\le m-1,\end{array}$

which means the eigenvalues of J = −cI_N+Ŵ are

$\begin{array}{l}{\lambda }_{jk}^J=-c\%+{\lambda }_{jk}^{\widehat{W}}\\ +{\lambda }_j^{{\widehat{W}}^0}+{\rho }_m^k{\lambda }_j^{{\widehat{W}}^1}+{\rho }_m^{2k}{\lambda }_j^{{\widehat{W}}^2}+{\rho }_m^{-2k}{\lambda }_j^{{\widehat{W}}^{m-2}}+{\rho }_m^{-k}{\lambda }_j^{{\widehat{W}}^{m-1}},\\ 0\le j\le n-1,0\le k\le m-1.\end{array}$

By direct calculation, we find that the eigenvalues of J are

$\begin{array}{l}{\lambda }_{jk}^J\\ =-[w_{1,0}{H}^{\prime }({\psi }_h)(1-{\rho }_n^j)+w_{-1,0}{H}^{\prime }(-{\psi }_h)(1-{\rho }_n^{-j})\\ +w_{0,1}{H}^{\prime }({\psi }_v)(1-{\rho }_m^k)+w_{0,-1}{H}^{\prime }(-{\psi }_v)(1-{\rho }_m^{-k})\\ +w_{1,1}{H}^{\prime }({\psi }_h+{\psi }_v)(1-{\rho }_n^j{\rho }_m^k)\\ +w_{-1,1}{H}^{\prime }({\psi }_v-{\psi }_h)(1-{\rho }_n^{-j}{\rho }_m^k)\\ +w_{1,-1}{H}^{\prime }({\psi }_h-{\psi }_v)(1-{\rho }_n^j{\rho }_m^{-k})\\ +w_{-1,-1}{H}^{\prime }(-{\psi }_h-{\psi }_v)(1-{\rho }_n^{-j}{\rho }_m^{-k})\\ +w_{2,0}{H}^{\prime }(2{\psi }_h)(1-{\rho }_n^{2j})\\ +w_{-2,0}{H}^{\prime }(-2{\psi }_h)(1-{\rho }_n^{-2j})\\ +w_{0,2}{H}^{\prime }(2{\psi }_v)(1-{\rho }_m^{2k})+w_{0,-2}{H}^{\prime }(-2{\psi }_v)(1-{\rho }_m^{-2k})],\end{array}$

for 0 ≤ j ≤ n−1, 0 ≤ k ≤ m−1.

The real parts of the eigenvalues are

$\begin{array}{l}\Re \left({\lambda }_{jk}^J\right)\\ =-[w_{1,0}{H}^{\prime }\left({\psi }_h\right)\left(1-cos\left(\frac{2\pi j}{n}\right)\right)\\ +w_{-1,0}{H}^{\prime }\left(-{\psi }_h\right)\left(1-cos\left(\frac{2\pi j}{n}\right)\right)\\ +w_{0,1}{H}^{\prime }\left({\psi }_v\right)\left(1-cos\left(\frac{2\pi k}{m}\right)\right)+w_{0,-1}{H}^{\prime }\left(-{\psi }_v\right)\\ \left(1-cos\left(\frac{2\pi k}{m}\right)\right)\\ +w_{-1,1}{H}^{\prime }\left({\psi }_v-{\psi }_h\right)\left(1-cos\left(\frac{2\pi j}{n}\right)cos\left(\frac{2\pi k}{m}\right)\right)\\ +w_{1,1}{H}^{\prime }\left({\psi }_h+{\psi }_v\right)\left(1-cos\left(\frac{2\pi j}{n}\right)cos\left(\frac{2\pi k}{m}\right)\right)\\ +w_{1,-1}{H}^{\prime }\left({\psi }_h-{\psi }_v\right)\left(1-cos\left(\frac{2\pi j}{n}\right)cos\left(\frac{2\pi k}{m}\right)\right)\\ +w_{-1,-1}{H}^{\prime }\left(-{\psi }_h-{\psi }_v\right)\left(1-cos\left(\frac{2\pi j}{n}\right)cos\left(\frac{2\pi k}{m}\right)\right)\\ +w_{2,0}{H}^{\prime }\left(2{\psi }_h\right)\left(1-cos\left(\frac{4\pi j}{n}\right)\right)& (11)\end{array}$

$\begin{array}{l}+w_{-2,0}{H}^{\prime }(-2{\psi }_h)\left(1-\cos \left(\frac{4\pi j}{n}\right)\right)\\ +w_{0,2}{H}^{\prime }(2{\psi }_v)\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)\\ +w_{0,-2}{H}^{\prime }(-2{\psi }_v)\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)].& (12)\end{array}$

Solutions will be stable whenever the real parts of the eigenvalues are negative.

In the following, we consider solutions when the coupling is symmetric in each direction (horizontal, vertical, and diagonal). In this case, w_{±1, 0} = h₁, w_{0, ±1} = v₁, w_{±1, ±1} = d, w_{±2, 0} = h₂, and w_{0, ±2} = v₂.

Then the real parts of the eigenvalues are

$\begin{array}{l}\Re \left({\lambda }_{jk}^J\right)\\ =-2\{h_1{H^{\prime }}_{od{d}^{}}({\psi }_h)\left(1-\cos \left(\frac{2\pi j}{n}\right)\right)+v_1{H^{\prime }}_{odd}({\psi }_v)\left(1-\cos \left(\frac{2\pi k}{m}\right)\right)\\ +d\left[{H^{\prime }}_{odd}({\psi }_h+{\psi }_v)\left(1-\cos \left(\frac{2\pi k}{m}+\frac{2\pi j}{n}\right)\right)\right]\\ +d\left[{H^{\prime }}_{odd}({\psi }_h-{\psi }_v)\left(1-\cos \left(\frac{2\pi k}{m}-\frac{2\pi j}{n}\right)\right)\right]\\ +h_2{H^{\prime }}_{odd}(2{\psi }_h)\left(1-\cos \left(\frac{4\pi j}{n}\right)\right)+v_2{H^{\prime }}_{odd}(2{\psi }_v)\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)\},& (13)\end{array}$

for 0 ≤ j ≤ n−1, 0 ≤ k ≤ m−1.

In summary, equations (12) and (13) give the eigenvalues for any cluster solution as described in Proposition 1 and hence can be used to determine their stability. Recall that such solutions are referred to as (p_v, p_h) p-cluster solutions, and have constant phase differences ψ_h and ψ_v between adjacent cells in the horizontal and vertical directions, respectively.

We now describe some special cases, which give rise to particular patterns of activity.

Case 1 (Synchronous solution). For any values of n and m, the system admits the solution with ψ_h = ψ_v = 0. This corresponds to a 1-cluster solution (p_h = p_v = 1); hence all neurons are firing in phase. The eigenvalues in this case have

$\begin{array}{l}\Re \left({\lambda }_{jk}^J\right)=-2{H^{\prime }}_{odd}(0)[h_1\left(1-\cos \left(\frac{2\pi j}{n}\right)\right)+v_1\left(1-\cos \left(\frac{2\pi k}{m}\right)\right)\\ +2d\left(1-\cos \left(\frac{2\pi k}{m}\right)\cos \left(\frac{2\pi j}{n}\right)\right)\\ +h_2\left(1-\cos \left(\frac{4\pi j}{n}\right)\right)+v_2\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)].& (14)\end{array}$

This solution will be stable if H′(0) > 0.

Case 2 (Diagonal Stripes). If n and m have a common factor p, then it is possible to have a solution where ${\psi }_h={\psi }_v=\psi =\frac{2\ell \pi }{p}$, where l < p and gcd(l, p) = 1. This corresponds to a p-cluster solution (with p_h = p_v = p) where nearest neighbors (vertically and horizontally) are out of phase by ψ. In the 2D network this appears as a diagonal stripe pattern, where neurons are in-phase with their diagonal neighbors above-right and below-left. The eigenvalues in this case have

$\begin{array}{l}\Re \left({\lambda }_{jk}^J\right)\\ =-2\{{H^{\prime }}_{odd}(\psi )\left[h_1\left(1-\cos \left(\frac{2\pi j}{n}\right)\right)+v_1\left(1-\cos \left(\frac{2\pi k}{m}\right)\right)\right]\\ +d\left[{H^{\prime }}_{odd}(2\psi )\left(1-\cos \left(\frac{2\pi k}{m}+\frac{2\pi j}{n}\right)\right)+{H^{\prime }}_{odd}(0)\left(1-\cos \left(\frac{2\pi k}{m}-\frac{2\pi j}{n}\right)\right)\right]\\ {H^{\prime }}_{odd}(2\psi )\left[h_2\left(1-\cos \left(\frac{4\pi j}{n}\right)\right)+v_2\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)\right]\}.& (15)\end{array}$

There is also a p-cluster solution with ${\psi }_h=\psi =\frac{2\ell \pi }{p}$ and ψ_v = 2π−ψ. In the 2D network this appears as a diagonal stripe pattern, where neurons are in-phase with their diagonal neighbors above-left and below-right. The eigenvalues in this case are

$\begin{array}{l}\Re \left({\lambda }_{jk}^J\right)=-2\{{H^{\prime }}_{odd}(\psi )\left[h_1\left(1-\cos \left(\frac{2\pi j}{n}\right)\right)+v_1\left(1-\cos \left(\frac{2\pi k}{m}\right)\right)\right]\\ +d\left[{H^{\prime }}_{odd}(0)\left(1-\cos \left(\frac{2\pi k}{m}+\frac{2\pi j}{n}\right)\right)+{H^{\prime }}_{odd}(2\psi )\left(1-\cos \left(\frac{2\pi k}{m}-\frac{2\pi j}{n}\right)\right)\right]\\ {H^{\prime }}_{odd}(2\psi )\left[h_2\left(1-\cos \left(\frac{4\pi j}{n}\right)\right)+v_2\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)\right]\}.& (16)\end{array}$

Case 3 (Horizontal Stripes). If p is a factor of m, then the system admits a solution with ψ_h = 0 and ${\psi }_v=\frac{2\ell \pi }{p}$, where ℓ < p and gcd(ℓ, p) = 1. This is a p-cluster solution (with p_h = 1, p_v = p) where neurons are synchronized with all their horizontal neighbors and have phase difference ψ_v with their nearest vertical neighbors. Thus, these appear as horizontal stripes. The eigenvalues in this case have

$\begin{array}{l}\Re \left({\lambda }_{jk}^J\right)\\ =-2\{{H^{\prime }}_{odd}({\psi }_v)\left[v_1\left(1-\cos \left(\frac{2\pi k}{m}\right)\right)+2d\left(1-\cos \left(\frac{2\pi k}{m}\right)\cos \left(\frac{2\pi j}{n}\right)\right)\right]\\ +{H^{\prime }}_{odd}(2{\psi }_v)v_2\left(1-\cos \left(\frac{4\pi k}{m}\right)\right)\\ +{H^{\prime }}_{odd}(0)\left[h_1\left(1-\cos \left(\frac{2\pi j}{n}\right)\right)+h_2\left(1-\cos \left(\frac{4\pi j}{n}\right)\right)\right]\}.& (17)\end{array}$

Case 4 (Vertical Stripes). Similarly, if p is a factor of n, then the system admits a solution with ψ_v = 0 and ${\psi }_h=\frac{2\ell \pi }{p}$, where ℓ < p and gcd(ℓ, p) = 1. The eigenvalues will be the same as in Case 3 with the roles of v₁ and h₁ swapped and those of v₂ and h₂ swapped. This is a p-cluster solution (with p_h = p, p_v = 1) where neurons are synchronized with all their vertical neighbors and have phase difference ψ_h with their nearest horizontal neighbors. Thus, these appear as vertical stripes.

Figure 2 illustrates some of the possible cluster solutions that exist in a 6 × 6 network. In this figure, circles of the same color correspond to neurons that spike synchronously (in-phase), and thus belong to the same cluster. Figure 2A shows the synchronous solution. Figure 2B shows a 2-cluster diagonal stripe solution of the first type, i.e., with ψ_h = ψ_v = π. Figure 2C shows a 3-cluster horizontal stripe solution. Figure 2D shows a 2-cluster vertical stripe solution. Figure 2E shows a 3-cluster diagonal stripe solution of the second type, i.e., with ${\psi }_h=\frac{2\pi }{3}$, ${\psi }_v=\frac{4\pi }{3}$. Figure 2F shows a (2, 3) 6-cluster solution, corresponding to ${\psi }_h=\pi ,{\psi }_v=\frac{2\pi }{3}$. Each row splits into two clusters and each column into three, giving six clusters overall in the network. Similarly Figure 2G shows a (2, 6) 6-cluster solution, corresponding to ${\psi }_h=\pi ,{\psi }_v=\frac{\pi }{3}$ and Figure 2H shows a (3, 6) 6-cluster solution, corresponding to ${\psi }_h=\frac{2\pi }{3},{\psi }_v=\frac{\pi }{3}$.

FIGURE 2

Figure 2. Illustration of cluster solutions in a 6 × 6 network. Circles that are the same color represent neurons that are in the same cluster. (A) synchronous solution. (B) 2-cluster diagonal stripe solution. (C) 3-cluster horizontal stripe solution. (D) 2-cluster vertical stripe solution. (E) 3-cluster diagonal stripe solution. (F) (2, 3) 6-cluster solution. (G) (2, 6) 6-cluster solution. (H) (3, 6) 6-cluster solution.

4. Applications and numerical results

In this section, we consider existence and stability of cluster solutions in the 2D, biophysical inhibitory network of Wang-Buzsaki neurons described in Section 2.2.

4.1. Homogeneous first nearest neighbor coupling—effect of diagonal coupling

We first consider the case of first nearest neighbor coupling only, so h₂ = v₂ = 0. To study the effect of diagonal coupling with otherwise homogeneous coupling strengths, we choose either h₁ = v₁>0, d = 0 or h₁ = v₁ = d>0. Note that these cases correspond to each neuron being coupled to four (left, right, up, down) or eight (left, right, up, down, four diagonal) neighboring neurons, respectively. In our model, $H_{odd}^{\prime }(\varphi )<0$ for $\varphi \in (0,\frac{17\pi }{32})\cup (\frac{47\pi }{32},2\pi )$ and $H_{odd}^{\prime }(\varphi )>0$ for $\varphi \in \left[\frac{17\pi }{32},\frac{47\pi }{32}\right]$ (Figure 3). In the following, we consider the stability of the solutions described in Section 3 above.

FIGURE 3

Figure 3. (A) The interaction function for the phase model, H(ϕ). (B) The odd part of the interaction function, H_odd(ϕ). Both are numerically computed for the biophysical inhibitory network of Wang-Buzsaki neurons (calculated using XPPAUT). (C,D) The derivatives of these functions, calculated using finite differences. For computation of stability conditions, we note that $H_{odd}^{\prime }(0)\approx -0.11,H_{odd}^{\prime }(\frac{\pi }{3})\approx -1.14$, $H_{odd}^{\prime }(\frac{\pi }{2})\approx -0.18,H_{odd}^{\prime }(\frac{2\pi }{3})\approx 0.78$, $H_{odd}^{\prime }(\pi )\approx 1.67$.

4.1.1. Diagonal stripes

As shown in Section 3, if n and m are divisible by p, the system admits p-cluster diagonal stripe solutions with phase differences ψ_h = ψ_v = ψ and ψ_h = ψ, ψ_v = 2π−ψ. Consideration of Equations (15), (16) with h₂ = v₂ = 0 and d = 0 shows that these diagonal stripe solutions will be stable if $H_{odd}^{\prime }(\psi )>0$ and unstable if $H_{odd}^{\prime }(\psi )<0$. For example, the 2- and 3-cluster diagonal stripe solutions (corresponding to ψ = π and $\psi =\frac{2\pi }{3}$ or $\frac{4\pi }{3}$, respectively) will be stable if they exist in a given network, while the 4- and 6-cluster solutions ($\psi =\frac{\pi }{2}$ or $\frac{3\pi }{2}$ and $\psi =\frac{\pi }{3}$ or $\frac{5\pi }{3}$, respectively) will be unstable.

Setting d > 0 in Equations (15), (16) shows that the addition of diagonal coupling can destabilize a stable diagonal stripe solution if $min(h_1,v_1)H_{odd}^{\prime }(\psi )+d\left(H_{odd}^{\prime }(2\psi )+H_{odd}^{\prime }(0)\right)<0$. For example, the 2-cluster diagonal stripe solution can be destabilized in our network model if

$\begin{array}{ll}d>-\frac{H_{odd}(\pi )}{2H_{odd}^{\prime }(0)}min(h_1,v_1)\approx 7.59min(h_1,v_1),& (18)\end{array}$

but the 3-cluster diagonal stripe solutions cannot be destabilized as $H_{odd}^{\prime }(\frac{4\pi }{3})+H_{odd}^{\prime }(0)>0.$ Diagonal cluster solutions that are unstable with d = 0 need to be considered individually. Analysis of the eigenvalues (15)-(16) shows that with homogeneous coupling (h₁ = v₁ = d), the 4-cluster and 6-cluster diagonal stripe solutions are unstable in the 4 × 4 and 6 × 6 networks, respectively.

4.1.2. Horizontal/vertical stripes

As shown in Section 3, if m is divisible by p, the system admits a p-cluster horizontal stripe solution with phase differences ψ_h = 0 and ψ_v≠0. Since $H_{odd}^{\prime }(0)<0$ for our model, if v₂ = h₂ = 0 and d = 0 in Equation (17) the eigenvalues ${\lambda }_{j0}^J$ will have positive real part. So without diagonal coupling, we expect the horizontal stripe solutions (and similarly the vertical stripe solutions) to be unstable. However, if d>0, then these eigenvalues have real part

$\Re \left({\lambda }_{j0}^J\right)=-2\left[2d{H}_{odd}^{\prime }\left({\psi }_v\right)+h_1{H}_{odd}^{\prime }\left(0\right)\right]\left(1-cos\left(\frac{2\pi j}{n}\right)\right).$

Thus, if

$\begin{array}{ll}2d{H}_{odd}^{\prime }({\psi }_v)+h_1{H}_{odd}^{\prime }(0)>0,& (19)\end{array}$

the p-cluster horizontal stripe solutions may be stable. From Figure 3, it is clear that this will only be possible for a restricted range $~(\frac{\pi }{2},\frac{3\pi }{2})$ of values of ψ_v. In particular, diagonal coupling may be able to stabilize the 2-cluster and 3-cluster horizontal stripe solutions, but not the 4-cluster or 6-cluster ones. Detailed analysis of the eigenvalues shows that for homogeneous coupling, h₁ = v₁ = d, the 2-cluster solution is stable in the 4 × 4 network and both the 2- and 3-cluster solutions are stable in the 6 × 6 network. A similar analysis applies to the vertical stripe solutions.

4.1.3. Other cluster solutions

As described in Proposition 3, if p_h divides n and p_v divides m, then the system admits a (p_h, p_v) p-cluster solution with p = lcm(p_h, p_v). Here, we focus on the cases with p_h≠p_v and min{p_h, p_v}>1. If a network is symmetric (n = m) we may assume that p_h<p_v, as similar stability results for p-cluster solutions with p_h > p_v can be obtained by interchanging the horizontal and vertical directions.

In the case of nearest neighbor coupling only (i.e., v₂ = h₂ = d = 0), we know from Equation (14) that this p-cluster is stable if $H_{odd}^{\prime }({\psi }_h)>0$ and $H_{odd}^{\prime }({\psi }_v)>0$, and this solution is unstable if $H_{odd}^{\prime }({\psi }_h)<0$ or $H_{odd}^{\prime }({\psi }_v)<0$. For instance, suppose that p_h = 2 (for which n has to be even), then ψ_h = π and $H_{odd}^{\prime }(\pi )>0$. In 2 × 4 or 4 × 4 networks, p_v = 1, 2 or 4. Since we assume that p_h ≠ p_v and p_v > 1, p_v = 4. In this case, we have a (2, 4) 4-cluster solution (since lcm(2, 4) = 4). However, this solution will be unstable as $H_{odd}^{\prime }({\psi }_v)<0$ for all possible ψ_v. Similarly, in 2 × 6, 4 × 6 or 6 × 6 networks, we must have p_v = 3 or 6. If p_v = 3, it leads to a (2, 3) 6-cluster solution, and this solution is stable, since $H_{odd}^{\prime }({\psi }_v)>0$ for all the possible ψ_v. If p_v = 6, it leads to a (2, 6) 6-cluster solution. However, this solution will be unstable as $H_{odd}^{\prime }({\psi }_v)<0$ for all possible ψ_v. In 6 × 6 networks we can also have a (3, 6) 6-cluster solution which is similarly unstable.

With the diagonal coupling, d > 0, a stable (or unstable) solution can be maintained if $H_{odd}^{\prime }({\psi }_h\pm {\psi }_v)>0$ (or $H_{odd}^{\prime }({\psi }_h\pm {\psi }_v)<0$). Otherwise adding diagonal coupling may change the stability. In 4 × 4 networks, regarding the (2, 4) 4-cluster solution discussed before, adding diagonal neighbors could not reverse the instability as $H_{odd}^{\prime }({\psi }_h\pm {\psi }_v)<0$. In contrast, in 6 × 6 networks, for the (2, 3) 6-cluster solution that is stable without diagonal coupling, adding diagonal neighbors tends to destabilize the solution, and the stability could be lost if d ≫ max{h₁, v₁}. For the (2, 6) 6-cluster solutions, which are unstable without diagonal coupling, including diagonal neighbors tends to stabilize the solution, and may change the solution from unstable to stable if d ≫ v₁. For the (3, 6) 6-cluster solutions adding diagonal neighbors brings in both stabilizing and destabilizing terms. Further analysis of the eigenvalues with homogeneous coupling (h₁ = v₁ = d) shows that in the 6 × 6 network the (2, 3) solution is stabilized, the (2, 6) is destabilized and the (3, 6) solution remains unstable.

We illustrate the stability of these multiple cluster solutions in numerical solutions of inhibitory networks of Wang-Buzsaki neurons. In these simulations, synaptic strength is homogeneous at the maximum conductance level h₁ = v₁ = g_syn, and d = g_syn when it is non-zero.

4.1.4. Numerical simulations in 4 × 4 networks

For a 4 × 4 network with h₁ = v₁ = d = g_syn, 4 different cluster solutions are stable: 2-cluster horizontal, vertical, and diagonal stripe solutions, and a 4-cluster solution in which groups of 4 cells fire synchronously with a $\frac{\pi }{2}$ phase difference between successively firing groups. We show the numbering and coloring of our 4 × 4 and 6 × 6 networks in Figure 4. When the cells are numbered consecutively across rows of the 4 × 4 network lattice, the clusters in this 4-cluster solution consist of cells {1, 3, 9, 11}, {2, 4, 10, 12}, {6, 8, 14, 16}, and {5, 7, 13, 15}. This solution differs from those discussed so far since the horizontal and vertical phase differences ψ_h and ψ_v are not the same for all cells but instead alternate between $\frac{\pi }{2}$ and $-\frac{\pi }{2}$ between neighboring cells. The stability of these 4 solutions in the same network is demonstrated in Figure 5. In each simulation, the network is initialized in one of the stable solutions (A: 2-cluster horizontal stripe; B: 2-cluster diagonal stripe) and a transient perturbation, in the form of increased external current applied to a subset of cells, switches the network to a different stable solution (A: 2-cluster vertical stripe; B: 4-cluster solution). The corresponding 2D clusters are shown above each raster plot. If the diagonal coupling is removed (d = 0) other simulations (not shown) indicate that the 2-cluster horizontal stripe solution becomes unstable, while the 2-cluster diagonal stripe solution remains stable.

FIGURE 4

Figure 4. Illustration of the 2D configuration of cells for numerical simulations: 4 × 4 network (left) and 6 × 6 network (right). In both, the cells are numbered consecutively across rows of the network lattice, and odd numbered cells are a lighter shade compared to even numbered cells. In the figures that follow, the raster plots use the same colors.

FIGURE 5

Figure 5. Multiple stable cluster solutions in 4 × 4 inhibitory neural networks of Wang-Buzsaki model neurons. In the network cells are coupled to their first nearest neighbors in the horizontal, vertical, and diagonal directions (h₁ = v₁ = d = g_syn, h₂ = v₂ = 0) (A) Network was initiated in the 2-cluster horizontal stripe solution and a transient perturbation (additional applied current to a subset of cells) was given over the time interval t ∈ [1500, 1800] (between dashed vertical lines) which switches the solution to the 2-cluster vertical stripe solution. (B) Network was initiated in the 2-cluster diagonal stripe solution and a similar transient perturbation to a different subset of cells switched the solution to the stable 4-cluster solution. Above each raster plot, we show the clusters in 2D at the beginning and end. Colors inside the circle correspond to the colors in the raster plot. Circles with the same background color represent neurons that are in the same cluster.

4.1.5. Numerical simulations in 6 × 6 networks

In the 6 × 6 network, 2-, 3-, and 6-cluster solutions are stable. For example, without diagonal coupling in the network (d = 0), both 2- and 3-cluster diagonal stripe solutions are stable and transient perturbations can switch the network between these stable solutions (Figure 6A). In this simulation, the network is initialized in the 3-cluster diagonal stripe solution and at time t∈[1500, 1800] ms the applied current is transiently increased to a subset of cells. This perturbation switches the solution to the stable 2-cluster diagonal stripe solution. A perturbation of the same magnitude to a different subset of cells switches the 3-cluster diagonal stripe solution to a (2, 3) 6-cluster solution such that ψ_h = π and ${\psi }_v=\frac{2\pi }{3}$ (Figure 6B). The corresponding 2D clusters are shown above each raster plot. In the (2, 3) 6-cluster solution, cells across rows break into 2 clusters while cells along columns separate into 3 clusters, as illustrated in Figure 2F. In this simulation, the evolution to the stable (2, 3) 6-cluster solution took longer following the transient perturbation, compared to the evolution to the 2-cluster diagonal stripe solution. Generally, we expect that the evolution time to a new stable solution will depend on the specific perturbation given to induce the switch of solutions.

FIGURE 6

Figure 6. Multiple stable cluster solutions in 6 × 6 inhibitory neural network. Cells are coupled to their first nearest neighbors in the horizontal and vertical directions only (h₁ = v₁ = g_syn, d = h₂ = v₂ = 0) (A) Network is initialized in a stable 3-cluster diagonal stripe solution. A transient perturbation to a subset of cells, given at t ∈ [1500, 1800] ms, switches the solution to a stable 2-cluster diagonal stripe solution. (B) Network is initialized in a stable 3-cluster diagonal stripe solution and a transient perturbation to a different subset of cells switches the solution to a stable (2, 3) 6-cluster solution. Above each raster plot, we show the clusters in 2D at the beginning and end. Colors inside the circle correspond to the colors in the raster plot. Circles with the same background color represent neurons that are in the same cluster.

In the 6 × 6 network with the addition of diagonal first nearest neighbor coupling (i.e., d = g_syn), both 2- and 3-cluster horizontal and vertical stripe solutions are stable, and transient perturbations can switch the network between these solutions (Figure 7). For example, for a network in the 3-cluster vertical stripe solution, a transient perturbation to a subset of cells switches the solution to the stable 2-cluster horizontal stripe solution (Figure 7A). In the same network initialized in the 3-cluster horizontal stripe solution, a transient perturbation to a different subset of cells switches the network to a stable (2, 6) 6-cluster solution with ψ_h = π and ${\psi }_v=\frac{\pi }{3}$ (Figure 7B). In this solution, cells across rows break into 2 clusters while cells along columns separate into 6 clusters, as illustrated in Figure 2G.

FIGURE 7

Figure 7. Multiple stable cluster solutions in the 6 × 6 network. Cells are coupled to their first nearest neighbors in the horizontal, vertical, and diagonal directions (h₁ = v₁ = d = g_syn, h₂ = v₂ = 0). (A) Network is initialized in a stable 3-cluster vertical stripe solution. A transient perturbation to a subset of cells, given at t ∈ [1500, 1800] ms, switches the solution to a stable 2-cluster horizontal stripe solution. (B) Network is initialized in a stable 3-cluster horizontal stripe solution and a transient perturbation to a different subset of cells switches the solution to a stable (2, 6) 6-cluster solution.

4.2. Effect of second nearest neighbor homogeneous coupling

Now we study the effect on the stability of the cluster solutions of adding second nearest coupling with homogeneous coupling strengths, i.e., h₂ = v₂ = d = h₁ = v₁>0. For notational simplicity, we denote the first nearest neighbor coupling considered in the last section as 8 nearest neighbor coupling.

4.2.1. Diagonal stripes

As discussed in Section 4.1, if both m and n are even, the system admits a 2-cluster diagonal stripe solution with constant phase difference ψ_h = ψ_v = π. Adding the second nearest neighbor coupling tends to destabilize the diagonal stripes since $H_{odd}^{\prime }(2{\psi }_h)=H_{odd}^{\prime }(2{\psi }_v)=H_{odd}^{\prime }(0)<0$.

If m and n are divisible by 3, there is a 3-cluster diagonal stripe solution. In our model, this solution is stable with 8 nearest neighbors coupling. Adding second nearest neighbor coupling will increase the stability as $H_{odd}^{\prime }(2{\psi }_h)=H_{odd}^{\prime }(2{\psi }_v)=H_{odd}^{\prime }(2\pi /3)>0$. Thus, the 3-cluster diagonal solution will be stable with any strengths of second nearest neighbor coupling.

If n and m are divisible by 6, there is a 6-cluster diagonal stripe solution. In our model, this solution is unstable with 8 nearest neighbors coupling. Adding second nearest neighbor coupling may increase (${\psi }_h=\frac{\pi }{3}$) or decrease (${\psi }_h=\frac{5\pi }{3}$) the stability.

Further analysis of the eigenvalues in Equations (15), (16) for the 6 × 6 network with homogeneous coupling (h₂ = v₂ = d = h₁ = v₁>0), shows that the 2-cluster diagonal stripe solutions remain stable, the 3-cluster solutions remain stable and the 6-cluster solutions remain unstable.

4.2.2. Horizontal/vertical stripes

If m is divisible by p, the system admits a p-cluster horizontal stripe solution with phase difference ψ_h = 0 and ψ_v ≠ 0. Including the second nearest neighbor coupling may or may not change the stability of horizontal stripe solutions. For example, in 6 × 6 networks, there is a 6-cluster horizontal stripe solution with ${\psi }_v=\frac{\pi }{3}$ or $\frac{5\pi }{3}$. This solution is unstable regardless of the presence of the second nearest neighbor coupling since there is always an eigenvalue with positive real part (From Equation 17, $\Re ({\lambda }_{0k}^J)>0,k=1,2,\dots ,5$). The same network also admits 2-cluster and 3-cluster horizontal stripe solutions, which are stable with 8 nearest neighbor coupling (see Section 4.1). Evaluation of the eigenvalues in Equation (17) shows that homogeneous second nearest neighbor coupling (h₂ = v₂ = h₁ = v₁ = d) cannot destabilize either of these solutions. However, we shall see in the next section that heterogeneous coupling could destabilize these solutions, e.g., if h₂ ≫ v₂. Similar results are expected for vertical stripe solutions by swapping the horizontal and vertical directions.

4.2.3. Other cluster solutions

If a (p_h, p_v) p-cluster solution is stable with 8 nearest neighbors coupling, then adding the second nearest neighbors can strengthen the stability if $H_{odd}^{\prime }(2{\psi }_h)>0$ and $H_{odd}^{\prime }(2{\psi }_v)>0$ but tends to destabilize the solution if $H_{odd}^{\prime }(2{\psi }_h)<0$ or $H_{odd}^{\prime }(2{\psi }_v)<0$. On the other hand, if a p-cluster solution with 8 nearest neighbors coupling is unstable, including the second nearest neighbors may change the stability. For instance, consider the (2, 3) 6-cluster solution in 6 × 6 networks. This solution could be stable in the absence of second nearest neighbor coupling if d≪min{h₁, v₁} but could be destabilized by adding second nearest neighbors if v₂ ≫ max{v₁, h₁}. Evaluation of the eigenvalues in Equation (13) shows that homogeneous second nearest neighbor coupling (h₂ = v₂ = h₁ = v₁ = d) does not change the stability of any of the (p_h, p_v) 6-cluster solutions in the 6 × 6 network.

4.2.4. Numerical simulations in 6 × 6 networks

When the coupling to the second nearest neighbors in horizontal and vertical directions are added to the 6 × 6 network (h₂ = v₂ = g_syn), $\text{2}\overline{/}$, $\text{3}\overline{/}$, and 4-cluster solutions are stable and transient perturbations can switch the network between these stable solutions (Figure 8). In Figure 8A, the network is initialized in a stable 2-cluster horizontal stripe solution. A transient perturbation to a certain subset of cells switches the 2-cluster horizontal stripe solution to a stable 4-cluster solution in which groups of 9 cells fire synchronously with a $\frac{\pi }{2}$ phase difference between successively firing groups. Similar to the 4-cluster solution shown in Figure 5B, this solution differs from those considered in Section 3 since the horizontal and vertical phase differences ψ_h and ψ_v alternate between $\frac{\pi }{2}$ and $-\frac{\pi }{2}$ between neighboring cells. In Figure 8B, the same network is initialized in the 3-cluster vertical stripe solution and a transient perturbation to a subset of cells switches this solution to a stable 2-cluster solution with non-uniform phase differences between neighboring cells. Thus, some cells have constant ψ_h and ψ_v, but other cells do not.

FIGURE 8

Figure 8. Stable cluster solutions in the 6 × 6 network. Cells are coupled to their first nearest neighbors in the horizontal, vertical, and diagonal directions, and to their second nearest neighbors in the horizontal and vertical directions (h₁ = v₁ = d = h₂ = v₂ = g_syn). (A) Network is initialized in a stable 2-cluster horizontal stripe solution. A transient perturbation to a subset of cells, given at t ∈ [1500, 1800] ms, switches the solution to a stable 4-cluster solution. (B) Transition from a stable 3-cluster vertical stripe solution to a stable 2-cluster solution after the similar transient perturbation is applied.

When the diagonal coupling is set to zero (d = 0) while keeping all other couplings non-zero in the 6 × 6 network, numerical simulation confirms that the 2-cluster horizontal stripe solution shown in Figure 8A becomes unstable, as predicted by evaluation of the eigenvalues in Equation (15) with n = m = 6, ψ_h = 0 and ψ_v = π.

4.3. Heterogeneous coupling strengths

The numerical simulations discussed so far have considered h₁ = v₁ = d = h₂ = v₂ for simplicity. However, our derived expressions for the eigenvalues e.g., Equation (13), do not depend on this assumption. We will consider the implications of this more general expression in this section.

First, we note that if the coupling strengths are slightly perturbed from this homogeneous case, then we don't expect that the stability should change. To see this, note that the eigenvalues in Equation 13 are linear in all the coupling strengths. So if the coupling strengths are perturbed from the homogeneous case, i.e., $h_1=A+B{ĥ}_1,v_1=A+B{\widehat{v}}_1,d=A+B\widehat{d},h_2=A+B{ĥ}_2,v_2=A+B{\widehat{v}}_2$, then the eigenvalues satisfy

$\Re \left({\lambda }_{jk}^J\right){|}_{h_1,v_1,\dots }=\Re \left({\lambda }_{jk}^J\right){|}_{A,A,\dots }+B\Re \left({\lambda }_{jk}^J\right){|}_{{ĥ}_1,\widehat{v_1},\dots }$

So if B is small enough the change in the eigenvalue will be small enough that the real part should not change sign.

If the coupling strengths are changed significantly, of course, the stability of solutions may change. In the rest of this section we use the stability conditions in Section 3 to study how introducing strongly heterogeneous coupling strengths may destabilize cluster solutions or give rise to new stable cluster solutions.

It was shown in Section 4.1 that with first nearest neighbor coupling only (h₁ = v₁ = d, h₂ = v₂ = 0), the 2-cluster diagonal stripe solution in the 6 × 6 network is stable, but with sufficiently strong diagonal coupling (d > ≈ 7.59h₁) this solution can be destabilized. See Equation (18). Numerical simulations of the network (not shown) confirm that with h₁ = v₁ = 1, h₂ = v₂ = 0, a value of d > 7.6 is required to destabilize this solution. Further evaluation of the eigenvalues given by Equation (15) for the 6 × 6 network shows that sufficiently strong diagonal and second nearest neighbor coupling (d = h₂ = v₂>h₁ = v₁) should also be able to destabilize the 2-cluster diagonal stripe solution. Numerical simulations of the Wang-Buzsaki inhibitory network confirming this are shown in Figure 9. In Figure 9A with homogeneous coupling strengths (h₁ = v₁ = d = h₂ = v₂ = 1), the network is initialized on the 2-cluster diagonal stripe solution and then returns to this solution after a transient perturbation to a certain subset of cells is applied at t ∈ [1500, 1600] in the network. This demonstrates the stability of this solution. In Figure 9B, heterogeneous coupling with stronger second nearest neighbor and diagonal coupling is introduced (h₁ = v₁ = 1, d = h₂ = v₂ = 4). Here, the same perturbation applied to the 2-cluster diagonal solution switches the network to a stable 8-cluster solution.

FIGURE 9

Figure 9. Stable cluster solutions in the 6 × 6 network. Cells are coupled to their first nearest neighbors in the horizontal, vertical, and diagonal directions, and to their second nearest neighbors in the horizontal and vertical directions. (A,B) The effect of heterogeneous coupling strengths on the stability of the 2-cluster diagonal stripe solution; (A) With homogeneous coupling strengths (h₁ = v₁ = d = h₂ = v₂ = 1) the network returns to the 2-cluster diagonal stripe solution after a sufficiently small perturbation is applied for t ∈ [1500, 1600] ms. (B) With stronger coupling to the first nearest neighbors in the diagonal direction and to the second nearest neighbors in horizontal and vertical directions (h₁ = v₁ = 1, d = h₂ = v₂ = 4), the same perturbation switches the network to a stable 8-cluster solution. (C) Stability of (3, 6) 6-cluster solution with weaker coupling to the first nearest neighbors in the vertical direction (h₁ = d = h₂ = v₂ = 1, v₁ = 0.4). The network starts on the solution and then returns to it after a transient perturbation applied for t ∈ [1500, 1800] ms. The corresponding 2D clusters for each stable solution after transient perturbations are shown next to each raster plot.

With homogeneous coupling strengths in a 6 × 6 network, (p_h, p_v) 6-cluster solutions with p_v = 6 are unstable (Section 4.1). However, if the horizontal coupling to the first nearest neighbors becomes stronger than the vertical coupling (e.g., h₁ = d = h₂ = v₂ = 1 but v₁ = 0.4), this heterogeneity gives rise to a stable (3, 6) 6-cluster solution such that ${\psi }_h=\frac{2\pi }{3}$ and ${\psi }_v=\frac{\pi }{3}$ (Figure 9C). In this solution, cells across rows break into 3 clusters while cells along columns separate into 6 clusters, as illustrated in Figure 2H.

We can further break the coupling symmetry in our networks if we relax the assumption that the connectivity matrix W is symmetric, that is, if we let w_{i, j} ≠ w_{j, i}. This returns us to the general stability condition (12). Note that this condition depends on H, not H_odd and, from Figures 3B,D, there is asymmetry in H(ϕ). This could lead to interesting results if W is not symmetric.

5. Discussion

We have shown the existence of synchronized cluster solutions and found their stability conditions in networks of intrinsically oscillating neurons on a general 2D torus. The stability conditions are derived generally for potential coupling to nearest neighbors in the vertical and horizontal directions, nearest neighbors in the diagonal directions, and second nearest neighbors in the vertical and horizontal directions. We considered specific solutions arising from various connectivity configurations and different coupling strengths. We also considered heterogeneous coupling strengths. We showed that cluster solutions could change stability if we add or remove connections between nearest neighbors in the diagonal directions or second nearest neighbors in the vertical and horizontal directions. We then used simulations of biophysical inhibitory neural networks (4 × 4 and 6 × 6) to verify our analytical predictions about which solutions would be stable under various configurations. See Tables 2, 3 for a summary of the results with homogeneous coupling. For some solutions, the clusters in the 2D network can be thought of as generalizations of cluster solutions that occur in 1D networks, where isolating a ring in the horizontal or vertical direction results in a solution with constant phase differences between neighboring cells. Examples of this would be the diagonal, horizontal, and vertical stripe solutions. However, the 2D nature of the coupling introduces other solutions that are not simple generalizations of 1D solutions. For example we found p-cluster solutions with p_h clusters in the horizontal direction and p_v clusters in the vertical direction. Further, in our numerical simulations we observed stable solutions where the phase differences between neighboring cells are not uniform across the network but do follow a regular pattern, such as in Figures 5B, 8A, and stable solutions where the number of neurons in different clusters is different, such as in Figure 9B. While the stability of some of these observed solutions makes sense intuitively, we have not yet been able to find general stability conditions for these types of solution.

TABLE 2

Table 2. Cluster solutions for the 4 × 4 network and the stability as predicted by the phase model analysis for various homogeneous coupling arrangements.

TABLE 3

Table 3. Cluster solutions for the 6 × 6 network and the stability as predicted by the phase model analysis for various homogeneous coupling arrangements.

Our existence and stability analyses apply to general m×n 2D networks. While our numerical simulations focused on square (m = n) 2D networks, similar types of solutions will be exhibited in rectangular (m ≠ n) as long as the existence and stability conditions are met. For example, 2-cluster horizontal (vertical) stripe solutions will exist in any size network that has an even number of rows (columns). Similarly, 2-cluster (3-cluster) diagonal stripe solutions will exist in any size network where the number of rows and columns is even (divisible by 3). Additionally, the types of (p_h, p_v) p-cluster solutions that we observed in square networks will exist in rectangular networks as long as the number of columns is divisible by p_h and the number of rows is divisible by p_v. The stability of any of these solutions can be determined using our expressions for the eigenvalues.

Diagonal coupling seems to be particularly important in determining the stable cluster solutions and differentiates the 2D network from the 1D network. Without diagonal coupling, if either the p_h- or p_v-cluster solution is unstable for the 1D network, then the (p_h, p_v) p-cluster solution is unstable for the 2D network. However, adding diagonal coupling increases the number of predicted stable solutions. For example, in the 6 × 6 network, the number of stable solutions increases from 9 to 15 (see Table 3). We can also find new types of p-cluster solutions for the 2D network. Even if p-cluster solutions cannot exist in 1D networks with n and m neurons, there can exist p-cluster solutions in the 2D m×n networks if there are p_h, p_v (which divide m, n, respectively) such that p = lcm(p_h, p_v).

While our analysis results can be applied to any type of network, we chose to focus on networks of inhibitory neurons. Intuitively, one might expect that the most stable configurations in such networks would be when neurons that are connected don't fire synchronously, equivalently, that the clusters would contain neurons that are not directly coupled. Indeed we found a very stable 4-cluster solution in both the 4 × 4 and 6 × 6 networks where neurons in the same cluster were minimally connected (see Figures 5B, 8A). However, we also found stable solutions where neurons in the same cluster are directly coupled. The stability of many of these solutions was enabled by the 2D structure of the network and specifically the diagonal coupling discussed above. The relative phases of neurons in different clusters is an important factor in determining stability. In phase model analysis the effect of the relative phase on stability is determined by the slope of the interaction function, at the relative phase (see Figure 3). Not surprisingly, we found that solutions with relative phases of π are the most strongly stable; however solutions with any relative phases $\varphi \in [\frac{17\pi }{32},\pi ]$ are also stable. Interestingly, the most strongly destabilizing relative phase was not 0, but around $\frac{3\pi }{16}$. This makes sense if we keep in mind that coupling between cells is not instantaneous in our model network but through simulated chemical synapses which introduce a small delay. Such a small delay changes the phase when the inhibition is felt by the post-synaptic neuron. Stability of a solution in a network with a particular coupling structure occurs when the coupling strength is stronger between neurons with relative phases that are stabilizing than between neurons whose relative phases are destabilizing.

Cluster solutions in model networks have been linked to neural assemblies in the brain as these are solutions where the neurons form groups, with neurons in the same group firing synchronously and neurons in different groups firing phase-locked with non-zero phase difference. Our analysis results show how such assemblies can arise spontaneously, generated only by the connectivity of the network and the dynamics of the neurons. In our model of 2D networks of intrinsically oscillatory neurons, a vast range of cluster solutions can exist, with multiple stable solutions occurring for the same parameter values. For the small 6 × 6 inhibitory networks we studied in detail, with certain coupling configurations we found four different stable cluster solutions in our numerical simulations. Larger networks would be expected to have more such solutions. In these simulations, we were able to switch the network between solutions by adding transient input to a subset of the neurons. This gives one mechanism for the formation of neural assemblies, where the different stable cluster solutions represent different assemblies and transient external input to the network switches the system from one assembly to another. We note that the number of clusters in the solution affects the network level behavior. For example, a 2-cluster solution will give a network firing rate that is twice the intrinsic firing rate of the neurons. Thus, the co-existence of stable cluster solutions gives a mechanism for the network to respond to different transient inputs with different network firing rates. We also showed how changing the coupling configuration could switch the stable cluster solutions that occur. For example, adding diagonal coupling could increase the number of stable solutions and changing the relative strengths of different coupling connections could change the stability of cluster solutions. Thus, changes in network connectivity gives another mechanism for changing the number and composition of the neural assemblies that a network can exhibit. In our example simulations, neurons formed connections with between 4 and 12 of their neighbors. Thus, the connectivity ranged from sparse (~10% for 6 × 6 networks with 4 connections) to dense (~75% connectivity for 4 × 4 networks with 12 connections). Our results indicate that the structure of connectivity is more important for forming clusters than density of connectivity.

Neural assembly firing has been identified in the striatum, a sparsely coupled, inhibitory network that is part of the basal ganglia circuit (Carrillo-Reid et al., 2008; Miller et al., 2008; Adler et al., 2013; Barbera et al., 2016; Klaus et al., 2017). Based on an extensive review of anatomical and experimental evidence, Burke et al. (2017) proposed that lateral inhibition between spiny neurons in the striatum organizes the cells into assemblies, with different assemblies corresponding to different behaviors. Our work supports this proposal. In our model, the external input that switches between assemblies would correspond to excitatory input to the striatum from the cortex or thalamus. Further, the strength of synapses in the striatum can be affected by changes in dopamine level (Lemos et al., 2016; Dobbs et al., 2017), which would correspond to changing the connections strengths in our network. Of course our model is much simpler than the striatum, which includes multiple different populations of spiny neurons as well as several populations of interneurons (Burke et al., 2017). Future work could include extending our model to include more populations.

Recent work on memory formation suggests that some brain regions have pre-configured neural ensembles, allowing memories to be encoded quickly (Grosmark and Buzsáki, 2016; Farooq et al., 2019; Miyawaki and Mizuseki, 2022). The internal generation of multiple assemblies we find in our study gives a mechanism for this. Although the model networks in our simulations were inhibitory, the analysis does not depend on the type of coupling. A topic for future work would be to apply our results to excitatory or excitatory/inhibitory networks.

It would be interesting to see how our work extends to more biologically realistic coupling regimes. One obvious extension would be to take the continuum limit (i.e., m, n → ∞), resulting in a 2D partial differential equation for the phases of the neurons, ϕ(x, y). While there have been several studies of 1D continuum models (Ermentrout, 1992; Strogatz, 2000; Wiley et al., 2006; Kazanci and Ermentrout, 2007; Sethia et al., 2011; Girnyk et al., 2012; Omel'chenko et al., 2014; Heitmann and Ermentrout, 2015), we know of only one study of a 2D continuum model (Heitmann et al., 2012). Analysis of the existence and stability of wave solutions (the continuum analog of our phase-locked cluster solutions) has only been done in 1D models. This work has primarily focused on either homogeneous coupling or symmetric distance-dependent coupling, and either excitatory or center-surround coupling, where neurons receive excitation from near neighbors and inhibition from further neighbors.The extension to purely inhibitory 2D networks can in principle be done, although the analysis may be challenging. Our work on discrete networks and numerical studies of discrete and continuum models (Heitmann et al., 2012; Spreizer et al., 2017) indicates that it would be of particular interest to consider non-monotonic coupling and heterogeneous coupling. Further, some interesting behavior in our models results from the non-square matrices, which could be incorporated by choosing different length scales in a continuum model. Another extension of our work would be to relax or change the structure in the connectivity matrix. Other biologically relevant possibilities include small-world networks and networks with hubs (Bullmore and Sporns, 2009; Bassett and Bullmore, 2017).

Data availability statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: ModelDB, https://modeldb.yale.edu/267329.

Author contributions

SC, JM, and XW conducted the phase model analyses. VB and HR implemented and conducted numerical simulations. All authors contributed to conceptualization, methodology, writing the original draft, reviewing and editing the final manuscript. All authors contributed to the article and approved the submitted version.

Funding

All authors received support from the American Institute of Mathematics Structured Quartet Research Ensembles program. SC was supported by the Natural Sciences and Engineering Research Council of Canada.

Acknowledgments

We thank the American Institute of Mathematics for providing support through its Structured Quartet Research Ensembles program, from which this study was first initiated. SC acknowledges the support of the Natural Sciences and Engineering Research Council of Canada.

Conflict of interest

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

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Achuthan, S., and Canavier, C. C. (2009). Phase-resetting curves determine synchronization, phase locking, and clustering in networks of neural oscillators. J. Neurosci. 29, 5218–5233. doi: 10.1523/JNEUROSCI.0426-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Adler, A., Finkes, I., Katabi, S., Prut, Y., and Bergman, H. (2013). Encoding by synchronization in the primate striatum. J. Neurosci. 33, 4854–4866. doi: 10.1523/JNEUROSCI.4791-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Angulo-Garcia, D., Berke, J. D., and Torcini, A. (2016). Cell assembly dynamics of sparsely-connected inhibitory networks: a simple model for the collective activity of striatal projection neurons. PLoS Comput. Biol. 12, e1004778. doi: 10.1371/journal.pcbi.1004778

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashwin, P., and Swift, J. W. (1992). The dynamics of n weakly coupled identical oscillators. J. Nonlinear Sci. 2, 69–108. doi: 10.1007/BF02429852

CrossRef Full Text | Google Scholar

Barbera, G., Liang, B., Zhang, L., Gerfen, C. R., Culurciello, E., Chen, R., et al. (2016). Spatially compact neural clusters in the dorsal striatum encode locomotion relevant information. Neuron 92, 202–213. doi: 10.1016/j.neuron.2016.08.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, D. S., and Bullmore, E. T. (2017). Small-world brain networks revisited. Neuroscientist 23, 499–516. doi: 10.1177/1073858416667720

PubMed Abstract | 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

Burke, D. A., Rotstein, H. G., and Alvarez, V. A. (2017). Striatal local circuitry: a new framework for lateral inhibition. Neuron 96, 267–284. doi: 10.1016/j.neuron.2017.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, S. A., and Wang, Z. (2018). Phase models and clustering in networks of oscillators with delayed coupling. Phys. D Nonlinear Phenomena 363, 44–55. doi: 10.1016/j.physd.2017.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Canavier, C. C., Kazanci, F. G., and Prinz, A. A. (2009). Phase resetting curves allow for simple and accurate prediction of robust n: 1 phase locking for strongly coupled neural oscillators. Biophys. J. 97, 59–73. doi: 10.1016/j.bpj.2009.04.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrillo-Reid, L., Tecuapetla, F., Tapia, D., Hernández-Cruz, A., Galarraga, E., Drucker-Colin, R., et al. (2008). Encoding network states by striatal cell assemblies. J. Neurophysiol. 99, 1435–1450. doi: 10.1152/jn.01131.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Culp, J. M. (2021). Cluster solutions in networks of weakly coupled oscillators on a 2D square torus. Math. Appl. Sci. Eng. Vol. 2, 149–218. doi: 10.5206/mase/14147

CrossRef Full Text | Google Scholar

Davis, P. J. (2013). Circulant Matrices. Providence, RI: American Mathematical Society.

Google Scholar

Dobbs, L., Lemos, J. C., and Alvarez, V. A. (2017). Restructuring of basal ganglia circuitry and associated behaviors triggered by low striatal d2 receptor expression: implications for substance use disorders. Genes Brain Behav. 16, 56–70. doi: 10.1111/gbb.12361

PubMed Abstract | CrossRef Full Text | Google Scholar

Dragoi, G., and Buzsáki, G. (2006). Temporal encoding of place sequences by hippocampal cell assemblies. Neuron 50, 145–157. doi: 10.1016/j.neuron.2006.02.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., Fries, P., and Singer, W. (2001). Dynamic predictions: oscillations and synchrony in top-down processing. Nat. Rev. Neurosci. 2, 704–716. doi: 10.1038/35094565

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G., and Terman, D. (2010). Mathematical Foundations of Neuroscience. New York, NY: Springer. doi: 10.1007/978-0-387-87708-2

CrossRef Full Text | Google Scholar

Ermentrout, G. B. (1992). Stable periodic solutions to discrete and continuum arrays of weakly coupled nonlinear oscillators. SIAM J. Appl. Math. 52, 1665–1687. doi: 10.1137/0152096

CrossRef Full Text | Google Scholar

Farooq, U., Sibille, J., Liu, K., and Dragoi, G. (2019). Strengthened temporal coordination within pre-existing sequential cell assemblies supports trajectory replay. Neuron 103, 719–733. doi: 10.1016/j.neuron.2019.05.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Galán, R. F., Ermentrout, G. B., and Urban, N. N. (2006). Predicting synchronized neural assemblies from experimentally estimated phase-resetting curves. Neurocomputing 69, 1112–1115. doi: 10.1016/j.neucom.2005.12.055

CrossRef Full Text | Google Scholar

Girnyk, T., Hasler, M., and Maistrenko, Y. (2012). Multistability of twisted states in non-locally coupled kuramoto-type models. Chaos Interdiscipl. J. Nonlinear Sci. 22, 013114. doi: 10.1063/1.3677365

PubMed Abstract | CrossRef Full Text | Google Scholar

Golomb, D., Hansel, D., Shraiman, B., and Sompolinsky, H. (1992). Clustering in globally coupled phase oscillators. Phys. Rev. A 45, 3516. doi: 10.1103/PhysRevA.45.3516

PubMed Abstract | CrossRef Full Text | Google Scholar

Golomb, D., and Rinzel, J. (1994). Clustering in globally coupled inhibitory neurons. Phys. D Nonlinear Phenomena 72, 259–282. doi: 10.1016/0167-2789(94)90214-3

CrossRef Full Text | Google Scholar

Gray, C. M., Koening, P., Engel, G., and Singer, W. (1989). Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature 338, 334–337. doi: 10.1038/338334a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Grosmark, A. D., and Buzsáki, G. (2016). Diversity in neural firing dynamics supports both rigid and learned hippocampal sequences. Science 351, 1440–1443. doi: 10.1126/science.aad1935

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansel, D., Mato, G., and Meunier, C. (1993). Clustering and slow switching in globally coupled phase oscillators. Phys. Rev. E 48, 3470. doi: 10.1103/PhysRevE.48.3470

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, K., Csicsvari, J., Hirase, H., Dragoi, G., and Buzsáki, G. (2003). Organization of cell assemblies in the hippocampus. Nature 424, 552–556. doi: 10.1038/nature01834

PubMed Abstract | CrossRef Full Text | Google Scholar

Heitmann, S., and Ermentrout, G. B. (2015). Synchrony, waves and ripple in spatially coupled kuramoto oscillators with Mexican hat connectivity. Biol. Cybernet. 109, 333–347. doi: 10.1007/s00422-015-0646-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Heitmann, S., Gong, P., and Breakspear, M. (2012). A computational role for bistability and traveling waves in motor cortex. Front. Comput. Neurosci. 6, 67. doi: 10.3389/fncom.2012.00067

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoppensteadt, F., and Izhikevich, E. (1997). Weakly Connected Neural Networks. New York, NY: Springer-Verlag. doi: 10.1007/978-1-4612-1828-9

CrossRef Full Text | Google Scholar

Kazanci, F. G., and Ermentrout, B. (2007). Pattern formation in an array of oscillators with electrical and chemical coupling. SIAM J. Appl. Math. 67, 512–529. doi: 10.1137/060661041

CrossRef Full Text | Google Scholar

Kazanci, F. G., and Ermentrout, B. (2008). Wave formation through the interactions between clustered states and local coupling in arrays of neural oscillators. SIAM J. Appl. Dyn. Syst. 7, 491–509. doi: 10.1137/070699147

CrossRef Full Text | Google Scholar

Kilpatrick, Z. P., and Ermentrout, B. (2011). Sparse gamma rhythms arising through clustering in adapting neuronal networks. PLoS Comput. Biol. 7, e1002281. doi: 10.1371/journal.pcbi.1002281

PubMed Abstract | CrossRef Full Text | Google Scholar

Klaus, A., Martins, G. J., Paixao, V. B., Zhou, P., Paninski, L., and Costa, R. M. (2017). The spatiotemporal organization of the striatum encodes action space. Neuron 95, 1171–1180. doi: 10.1016/j.neuron.2017.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopell, N., and Ermentrout, G. (2002). “Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators,” in Handbook of Dynamical Systems, Vol 2: Toward Applications, ed B. Fiedler (Amsterdam: Elsevier), 3–54. doi: 10.1016/S1874-575X(02)80022-4

CrossRef Full Text | Google Scholar

Laurent, G., and Davidowitz, H. (1994). Encoding of olfactory information with oscillating neural assemblies. Science 265, 1872–1875. doi: 10.1126/science.265.5180.1872

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemos, J. C., Friend, D. M., Kaplan, A. R., Shin, J. H., Rubinstein, M., Kravitz, A. V., et al. (2016). Enhanced gaba transmission drives bradykinesia following loss of dopamine d2 receptor signaling. Neuron 90, 824–838. doi: 10.1016/j.neuron.2016.04.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y.-X., Wang, Y.-Q., and Miura, R. (2003). Clustering in small networks of excitatory neurons with heterogeneous coupling strengths. J. Comput. Neurosci. 14, 139–159. doi: 10.1023/A:1021902717424

PubMed Abstract | CrossRef Full Text | Google Scholar

Mancilla, J. G., Lewis, T. J., Pinto, D. J., Rinzel, J., and Connors, B. W. (2007). Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex. J. Neurosci. 27, 2058–2073. doi: 10.1523/JNEUROSCI.2715-06.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, B., Walker, A., Shah, A., Barton, S., and Rebec, G. (2008). Dysregulated information processing by medium spiny neurons in striatum of freely behaving mouse models of Huntington's disease. J. Neurophysiol. 100, 2205–2216. doi: 10.1152/jn.90606.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, J., Ryu, H., Teymuroglu, Z., Wang, X., Booth, V., and Campbell, S. A. (2015). “Clustering in inhibitory neural networks with nearest neighbor coupling,” in Applications of Dynamical Systems in Biology and Medicine, eds T. Jackson and A. Radunskaya (New York, NY: Springer), 99–121. doi: 10.1007/978-1-4939-2782-1_5

CrossRef Full Text | Google Scholar

Miyawaki, H., and Mizuseki, K. (2022). De novo inter-regional coactivations of preconfigured local ensembles support memory. Nat. Commun. 13, 1–21. doi: 10.1038/s41467-022-28929-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Muldoon, S. F., Soltesz, I., and Cossart, R. (2013). Spatially clustered neuronal assemblies comprise the microstructure of synchrony in chronically epileptic networks. Proc. Natl. Acad. Sci. U.S.A. 110, 3567–3572. doi: 10.1073/pnas.1216958110

PubMed Abstract | CrossRef Full Text | Google Scholar

Okuda, K. (1993). Variety and generality of clustering in globally coupled oscillators. Phys. D 63, 424–436. doi: 10.1016/0167-2789(93)90121-G

CrossRef Full Text | Google Scholar

Omel'chenko, O. E., Wolfrum, M., and Laing, C. R. (2014). Partially coherent twisted states in arrays of coupled phase oscillators. Chaos Interdiscipl. J. Nonlinear Sci. 24, 023102. doi: 10.1063/1.4870259

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponzi, A., and Wickens, J. (2010). Sequentially switching cell assemblies in random inhibitory networks of spiking neurons in the striatum. J. Neurosci. 30, 5894–5911. doi: 10.1523/JNEUROSCI.5540-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryu, H., Miller, J., Teymuroglu, Z., Wang, X., Booth, V., and Campbell, S. A. (2021). Spatially localized cluster solutions in inhibitory neural networks. Math. Biosci. 336, 108591. doi: 10.1016/j.mbs.2021.108591

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakaguchi, H., and Kuramoto, Y. (1986). A soluble active rotater model showing phase transitions via mutual entertainment. Prog. Theoret. Phys. 76, 576–581. doi: 10.1143/PTP.76.576

CrossRef Full Text | Google Scholar

Saraga, F., Ng, L., and Skinner, F. K. (2006). Distal gap junctions and active dendrites can tune network dynamics. J. Neurophysiol. 95, 1669–1682. doi: 10.1152/jn.00662.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwemmer, M., and Lewis, T. (2012). “The theory of weakly coupled oscillators,” in Phase Response Curves in Neuroscience, eds N. Schultheiss, A. Prinz, and R. Butera (New York, NY: Springer), 3–31. doi: 10.1007/978-1-4614-0739-3_1

CrossRef Full Text | Google Scholar

Sethia, G. C., Sen, A., and Atay, F. M. (2011). Phase-locked solutions and their stability in the presence of propagation delays. Pramana 77, 905–915. doi: 10.1007/s12043-011-0184-2

CrossRef Full Text | Google Scholar

Spreizer, S., Angelhuber, M., Bahuguna, J., Aertsen, A., and Kumar, A. (2017). Activity dynamics and signal representation in a striatal network model with distance-dependent connectivity. eNeuro 4. doi: 10.1523/ENEURO.0348-16.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Wang, X.-J., and Buzsáki, G. (1996). Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J. Neurosci. 16, 6402–6413. doi: 10.1523/JNEUROSCI.16-20-06402.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiley, D. A., Strogatz, S. H., and Girvan, M. (2006). The size of the sync basin. Chaos Interdiscipl. J. Nonlinear Sci. 16, 015103. doi: 10.1063/1.2165594

PubMed Abstract | CrossRef Full Text | Google Scholar

Zanette, D. H. (2000). Propagating structures in globally coupled systems with time delays. Phys. Rev. E 62, 3167. doi: 10.1103/PhysRevE.62.3167

PubMed Abstract | CrossRef Full Text | Google Scholar