Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 2, Number 3, September 2021, pp.194-218 https://doi.org/10.5206/mase/14147 CLUSTER SOLUTIONS IN NETWORKS OF WEAKLY COUPLED OSCILLATORS ON A 2D SQUARE TORUS JORDAN MICHAEL CULP Abstract. We consider a model for an N × N lattice network of weakly coupled neural oscillators with periodic boundary conditions (2D square torus), where the coupling between neurons is assumed to be within a von Neumann neighborhood of size r, denoted as von Neumann r-neighborhood. Using the phase model reduction technique, we study the existence of cluster solutions with constant phase differences (Ψh, Ψv) between adjacent oscillators along the horizontal and vertical directions in our network, where Ψh and Ψv are not necessarily to be identical. Applying the Kronecker product representation and the circulant matrix theory, we develop a novel approach to analyze the stability of cluster solutions with constant phase difference (i.e., Ψh, Ψv are equal). We begin our analysis by deriving the precise conditions for stability of such cluster solutions with von Neumann 1-neighborhood and 2 neighborhood couplings, and then we generalize our result to von Neumann r-neighborhood coupling for arbitrary neighborhood size r ≥ 1. This developed approach for the stability analysis indeed can be extended to an arbitrary coupling in our network. Finally, numerical simulations are used to validate the above analytical results for various values of N and r by considering an inhibitory network of Morris-Lecar neurons. 1. Introduction An average neuron forms and receives one to ten thousand synaptic connections for sending and receiving information. Since there are at least 1011 neurons in the human brain, there are thus 1014 ∼ 1015 synaptic connections that are formed in the brain. The summation of this input at the cellular level combine to allow neurons to perform complicated information processing and when considered as a whole brain, or neural region, allow for the complex cognitive task that we use to live our lives to be performed. It is assumed that it is the structure and details inherent in these connections that hold the secret to how these tasks manifest. Thus, understanding our brains ability to organize and create coherent patterns out of the collection of electrical activity from billions of coupled individual neurons is of much research interest. The modeling of coupled oscillators have been applied to study a number of biological and physical systems, for example [38, 28, 42, 34, 30, 49, 46, 23, 18]. The synchronization and cluster solutions, as defined below, in networks of large populations of neurons play an important role in various brain functions [43, 9, 41, 22, 16, 6, 27, 48, 14, 47, 44, 40, 45]. The theory of weakly coupled oscillators [13, 25, 31] is a classical tool for using dynamical systems theory to study oscillations in neural networks, where the phase reduction method has been utilized to reduce a model of a weakly coupled neural network to a phase model. This phase model is used to study the existence and stability of cluster solutions in the network of oscillators. These cluster solutions are Received by the editors 1 August 2021; accepted 13 September 2021; published online September 19 2021. 2010 Mathematics Subject Classification. 92B20,92B25,34D20. Key words and phrases. phase models, coupled oscillators, synchronization, phase-locking, clustering solutions, stability. 194 CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 195 phase locked solutions to the phase model where oscillators separate into subgroups. The oscillators within each subgroup are synchronized with zero phase difference, while those in different subgroups are phase locked. The existence and stability of cluster solutions in networks of weakly coupled oscillators has been studied within the context of various network topologies and coupling schemes, see, e.g., [2, 24, 32, 38, 28] and the references therein. In [13], oscillators are modeled as a nonlinear system with stable limit cycle with bidirectional nearest neighbor coupling, where that the existence and necessary conditions for the stability of phase locking behavior in the network are established. In [26] oscillators are also modeled as a nonlinear system with stable limit cycle with bidirectional nearest neighbor coupling on a ring of oscillators. A Hopf bifurcation on the ring network topology is also studied to obtain different types of stable oscillators. A review on the stability of cluster solutions of oscillators can be found in see [2, 5]. General networks of identical oscillators are considered in [15] where network symmetries are examined in order to obtain different types of phase locking behavior. The existence and stability of cluster solutions in linear array (chains) and ring networks with uni- or bi-directional coupling have been shown in [29, 2, 5, 11, 13]. Examples of analysis of phase models with time delayed coupling can be found in [7, 8]. Moreover, the study of cluster solutions in models with all to all coupling can be found in [2, 32, 37, 17, 21]. The known results for the analytical study of cluster solutions on 1D network topologies are more extensive than the results for the 2D network topologies. The study of two dimensional lattices has been studied in the case of finite square lattice [35, 33]. In [33], the existence and stability of rotating wave solutions exist on finite square lattices with “no flux” boundary conditions and nearest neighbor coupling is studied. In [35], phaselocked behavior is studied, as the two-dimensional (and higher dimensional) arrays, with nearest neighbor coupling, can be decomposed into two one-dimensional problems, if some conditions on the intrinsic frequencies are met. In [3], rotating wave solutions on infinite lattices are considered. In [4], sufficient conditions for local asymptotic stability of phase-locked solutions in coupled phase models on infinite lattices are considered. In this work, we aim to extend the known analytical results concerning the existence and stability of cluster solutions in a two dimensional network topologies to include two dimensional finite lattices with periodic boundary conditions. Our work differs mainly from the previously mentioned results in that we derive general conditions for the existence and stability of cluster solutions on our network topology without specificity to a particular type of pattern such as stable rotating solutions (spiral waves), and for larger (greater than first) nearest neighbor coupling. Furthermore, we are able to find a convenient representation of our Jacobian matrix by utilizing results from circulant matrix theory and the Kronecker product that allows for a novel approach in the analysis of the stability of the cluster solutions on a 2D network topology. The results in this work are shown for a general N×N lattice of weakly coupled oscillators with periodic boundary conditions, with a von Neumann 1 and 2 neighborhood coupling size. We allow for further generalization of this model by considering situations where the horizontal and vertical phase differences are not equal, and where the neighborhood size can be extended to a general size r. 2. Notations In this section we provide some preliminary mathematical theory and notation to be used in this work. 196 J. CULP 2.1. Matrices. Define a circulant matrix generated by constants c1,c2,c3, · · · ,cN−1,cN as circ(c1,c2,c3, · · · ,cN−1,cN ) =   c1 c2 c3 · · · cN−1 cN cN c1 c2 · · · cN−2 cN−1 cN−1 cN c1 . . . cN−3 cN−2 ... ... . . . . . . . . . ... c3 c4 · · · cN c1 c2 c2 c3 · · · cN−1 cN c1   Likewise we will define an m×m block circulant matrix generated by the n1 ×n2 matrices, M1,M2, M3, · · · ,Mm−1,Mm as bcirc(M1,M2, · · · ,Mm−1,Mm) =   M1 M2 M3 · · · Mm−1 Mm Mm M1 M2 · · · Mm−2 Mm−1 Mm−1 Mm M1 . . . Mm−3 Mm−2 ... ... . . . . . . . . . ... M3 M4 · · · Mm M1 M2 M2 M3 · · · Mm−1 Mm M1   The Kronecker product of n1 ×n2 matrix M1 and m1 ×m2 matrix M2 is denoted as M1 ⊗M2 =   m11M2 m12M2 · · · m1,n2M2 ... ... ... mn1,1M2 mn1,2M2 · · · mn1,n2M2   For more information on known results and properties of circulant matrices and the Kronecker product used thorough out this work, please refer to [10]. Now let us introduce the notation for some special N ×N matrices used throughout the paper. Let A1 = circ(0, 1, 0, · · · , 0), A−1 = circ(0, 0, 0, · · · , 1). It is easy to see that A−1 = A N−1 1 = (A1) −1. Moreover, as A1 is a circulant matrix, its spectrum is known explicitly. Specifically, the eigenvalues of A1 are {λ0,λ1, · · · ,λN−1} and the corresponding eigenvectors are {u0,u1, · · · ,uN−1}, where λj = ω j, uj = 1 √ N (1,ωj,ω2j, · · · ,ω(N−1)j)T , 0 ≤ j ≤ N − 1 with ω = exp(2π √ −1/N) is the primitive Nth root of the unity. Since A−1 = A N−1 1 , A−1uj = A N−1 1 uj = λ N−1 j uj. Note that λj = ω j and ωN = 1. We have λN−1j = ω −j = λ−1j and hence A−1uj = λ −1 j uj. It is straightforward to see that A −k 1 uj = λ −k j uj for k ∈ N. Lemma 2.1. Let p,q ∈{−(N − 1), · · · ,−2,−1, 0, 1, 2, · · · ,N − 1} be given. Then (A p 1 ⊗A q 1)(ui ⊗uj) = ω pi+qjui ⊗uj, 0 ≤ i,j ≤ N − 1. Proof. By the commutative property of the Kronecker product, (A p 1 ⊗A q 1)(ui ⊗uj) = (A p 1ui) ⊗ (A q 1uj) = (λ p i ui) ⊗ (λ q juj) = (λ p i λ q j)ui ⊗uj. for all 0 ≤ i,j ≤ N − 1. Since λi = ωi for all 0 ≤ i ≤ N − 1, λ p i λ q j = ω pi+qj. � CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 197 This result shows that, for any p,q ∈{−(N−1), · · · ,−2,−1, 0, 1, 2, · · · ,N−1}, Ap1⊗A q 1 has the same eigenvectors {ui⊗uj : 0 ≤ i,j ≤ N−1} and the corresponding eigenvalues are {ωpi+qj : 0 ≤ i,j ≤ N−1}. In particular, (IN ⊗A1)(ui ⊗uj) = ωjui ⊗uj, (IN ⊗A−1)(ui ⊗uj) = ω−jui ⊗uj, (A1 ⊗ IN )(ui ⊗uj) = ωiui ⊗uj, (A−1 ⊗ IN )(ui ⊗uj) = ω−iui ⊗uj. (2.1) 2.2. Mapping the 2D Lattice into a 1D Array. Figure 1 shows how we map the indices of a two dimensional N ×N lattice into an one dimensional array from 1, . . . ,N2. Here, for 0 ≤ i,j ≤ N − 1, we use the mapping f : Z+ × Z+ → N defined by f(i,j) = jN + (i + 1). Clearly f is a one-to-one and onto mapping from {0, 1, 2, · · · ,N − 1}×{0, 1, 2, · · · ,N − 1} to {1, 2, . . . ,N2}. As we will see later, this mapping is used in the construction of our circulant coupling matrix, and allows us to conveniently represent a two dimensional network structure in one dimension. (0, 0) (N − 1, 0) (N − 1, 1)(0, 1) (0, 2) (0,N − 1) (N − 1,N − 1) f−→ 1 N 2NN + 1 2N + 1 (N − 1)N + 1 N2 Figure 1. Mapping the 2D Lattice into a 1D Array 2.3. von Neumann Neighborhoods. In this work, the von Neumann r-neighborhood on a 2D square lattice is defined as the set of rth adjacent nodes of a central node. As compared to the classical definition, the center is the common node is removed for convenience. Note that a von Neumann 1- neighborhood is all the nodes at a Manhattan distance of 1, and that an extension consists of taking the set of points at a Manhattan distance of r > 1. In Figure 2, an example of von Neumann 2-neighborhood around node P is plotted. The number of nodes in a von Neumann r-neighborhood is r2 + (r + 1)2 −1. In our example, r = 2 and we have 22 + (2 + 1)2 − 1 = 13 − 1 = 12, with 4 cells at distance 1 and 8 cells at distance 2. 2.4. Phase Reduction Method. In this subsection, we review the phase reduction method by which a general network model of identical weakly coupled oscillators can be reduced to a phase model. Consider a two dimensional network model of identical, weakly coupled oscillators on an N × N lattice with periodic boundary conditions. dyij dt = F(yij) + � N−1∑ ĩ,j̃=0 c(i,j),(̃i,j̃)G(yij,yĩj̃), 0 ≤ i,j ≤ N − 1. (2.2) Here yij ∈ Rm, 0 < � � 1 is the coupling strength, F is the vector field of the isolated oscillator, G is a function describing the coupling between oscillators and (c(i,j),(̃i,j̃)) is the coupling matrix, where as c(i,j),(̃i,j̃) ≥ 0 denotes the coupling between the oscillators yij and yĩj̃. Furthermore, we assume that for an isolated oscillator, the solution to dyij dt = F(yij) exhibits an exponentially asymptotically T−periodic 198 J. CULP 1 2 2 1 2 2 1 2 2 1 2 2 P Figure 2. The von Neumann 2-neighborhood of node P stable solution, denoted by X̂ij(t) on the domain 0 ≤ t ≤ T. Here, the frequency of the ijth isolated oscillator with period T is given by Ω = 2π T . Under the assumption of sufficiently weak coupling between oscillators, the theory of weakly coupled oscillators allows us to reduce the number of equations representing the dynamics of each oscillator to a single differential equation representing the change of the phase of each oscillator along the corresponding T-periodic limit cycle. This phase model will take the form, dθij dt = Ω + � N−1∑ ĩ,j̃=0 c(i,j),(̃i,j̃)H(θĩj̃ −θij), 0 ≤ i,j ≤ N − 1. where H is the interaction function with period 2π that satisfies H(θ) = 1 T ∫ T 0 Z(t) ·G ( X̂(t),X̂(t + θ/Ω) ) dt. Here Z is known as the phase response curve of the isolated oscillator, which is the unique periodic solution of the linearized adjoint system: dZ dt = −{DF(X̂)}TZ with the normalization 1 T ∫ T 0 Z(t) ·F(X̂)dt = 1 where DF(X̂) is the Jacobian of F with respect to X evaluated at X = X̂. The phase solution of the above phase model is θij(t) = Ωt + φij(t). (2.3) Here φij(t) represents the relative (initial) phase of the ij th oscillator. Under the assumption of weak coupling, the dynamics of the relative phase of each oscillator satisfies the differential equation (to the first order of �): dφij dt = � N−1∑ ĩ,j̃=0 c(i,j),(̃i,j̃)H(θĩj̃ −θij), (2.4) Note now that for two identical oscillators, θĩj̃(t),θij(t), by using the definition in (2.3) we would have that θĩj̃(t) −θij(t) = φĩj̃(t) −φij(t). Hence, in our work we will frequently refer to the phase model in CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 199 terms of dφij dt . A more rigorous mathematical discussion on the theory of weakly coupled oscillators can be found at [1, 19, 20, 39]. 3. Phase Model Analysis In this section, we study clustering dynamics of a neural network model on a square torus with von Neumann r-neighborhood coupling. Particularly, we analyze the existence of cluster solutions with constant phase difference and derive the precise condition for the stability of these solutions. Note that Cr := {(k,l) ∈ Z+ × Z+ : 0 < |k| + |l| ≤ r} is the set of indices in the von Neumann r-neighborhood on a square lattice, for r ∈ N. Throughout the paper, all the 2D indices are defined under modulo N unless otherwise stated. Consider a network model consisting of N × N identical weakly coupled neural oscillators on a torus with von Neumann r-neighborhood coupling: dθi,j dt = Ω + � ∑ 0<|k|+|l|≤r k,l∈Z wk,lH(θi+k,j+l −θi,j), 0 ≤ i,j ≤ N − 1 (3.1) where the coupling is assumed to be relative-location dependent, and wk,l = c(i+k,j+l),(i,j) is the coupling between the ijth oscillator and (i + k,j + l)th oscillator for all (i,j). 3.1. Existence of Cluster Solutions. Assume that the phase difference only depends on the relative location; that is, the phase difference between (i + k,j + l)th oscillator and the ijth oscillator is Ψk,l = θi+k,j+l −θi,j, ∀ 0 ≤ i,j ≤ N − 1. (3.2) By this assumption, there are two types of phase differences Ψ1,0 and Ψ0,1, where the first one is in the horizontal direction and the second one is in the vertical direction. For simplicity, we define Ψh = Ψ1,0 and Ψv = Ψ0,1. So in general Ψk,l = kΨh + lΨv. In view of (3.1) and (3.2), dΨk,l/dt = 0 for all (k,l). On the other hand, it follows from the periodic boundary condition that the horizontal and vertical phase differences should satisfy NΨh = NΨv = 0 mod 2π. Then NΨh = 2khπ, NΨv = 2kvπ for some kh,kv ∈ {0, 1, 2, · · · ,N − 1}. If Ψh = Ψv = 0 (i.e., kh = kv = 0), it gives a synchronization or one-cluster solution. If Ψh = 0, Ψv > 0 (i.e., kh = 0, kv ∈ N), it gives a solution of horizontal stripes. If Ψh > 0, Ψv = 0 (i.e., kh ∈ N, kv = 0) , it gives a solution of vertical stripes. Otherwise, suppose that Ψh > 0, Ψv > 0 (i.e., kh, kv ∈ N). Define ph = gcd(N,kh) and pv = gcd(N,kv). Let nd = N/pd and md = kd/pd with d = h,v. Then the solution with phase difference (Ψh, Ψv) gives us a cluster solution that is of period nh and nv along the horizontal and vertical directions, respectively, where Ψh = 2mhπ/nh and Ψv = 2mvπ/nv. It is clear that the period of this cluster solution is n := lcm(nh,nv). This solution defines an n-cluster solution with phase difference (Ψh, Ψv). Let NN = {2, 3, · · · ,N}. It leads to the following existence result. Theorem 3.1. The system (3.1) admits four types of cluster solutions with constant horizontal and vertical phase differences between adjacent oscillators. (i) There exists a synchronization solution with phase difference (Ψh, Ψv) = (0, 0). (ii) Assume that nv is a factor of N with nv ∈ NN . There exists an nv-cluster solution of horizontal stripes with phase difference (Ψh, Ψv) = (0, Ψv), where Ψv = 2mvπ/nv for some mv ∈ N satisfying mv < nv and gcd(mv,nv) = 1. (iii) Assume that nh is a factor with N with nh ∈ NN . There exists an nh-cluster solution of vertical stripes with phase difference (Ψh, Ψv) = (Ψh, 0), where Ψh = 2mhπ/nh for some mh ∈ N satisfying mh < nh and gcd(mh,nh) = 1. 200 J. CULP (iv) Assume that both nv and nh are factors of N with nh,nv ∈ NN . Define n = lcm(nh,nv). There exists an n-cluster solution with phase difference (Ψh, Ψv), where Ψd = 2mdπ/nd for some md ∈ N satisfying md < nd and gcd(md,nd) = 1 for d = h,v. . Furthermore, if Ψh = Ψv = Ψ, cluster solutions will have constant phase difference Ψ between adjacent oscillators in the horizontal and vertical directions. We denote such a solution as a cluster solution with constant phase difference Ψ. the following result. Corollary 3.2. The system (3.1) admits two types of cluster solutions with constant phase difference Ψ. (i) There exists a synchronization solution with Ψ = 0. (ii) Assume that n is a factor of N with n ∈ NN . There exists an n-cluster solution with constant phase difference Ψ, where Ψ = 2mπ/n for some m ∈ N satisfying m < n and gcd(m,n) = 1. 3.2. Stability of Cluster Solutions with Constant Phase Difference. In this subsection, we study the stability of cluster solutions and our investigation will be focused on cluster solutions with constant phase difference Ψ (between adjacent cells in the horizontal and vertical directions). 3.3. von Neumann 1-neighborhood Coupling: r = 1. First, let’s consider the case where the coupling is in the von Neumann 1-neighborhood. This is the case of the nearest neighbor coupling. Accordingly, model (3.1) becomes dθij dt = Ω + � ∑ |k|+|l|=1 k,l∈Z wk,lH(θi+k,j+l −θij) = Ω + � [ h−1H(θi−1,j −θij) + h1H(θi+1,j −θij) + v1H(θi,j+1 −θij) + v−1H(θi,j−1 −θij) ] , (3.3) where h±1 := w±1,0 and v±1 := w0,±1 for 0 ≤ i,j ≤ N − 1. Suppose that N ≥ 3 and n divides N. Consider an n-cluster solution with constant phase difference Ψ θ̄ (1) ij (t) = (Ω + �Ω (1) c )t + (i + j)Ψ, 0 ≤ i,j,≤ N − 1 (3.4) for some Ωc > 0. Substituting (3.4) into (3.3), we find that Ω(1)c = (h−1 + v−1)H(−Ψ) + (h1 + v1)H(Ψ). Let θij(t) = θ̄ (1) ij (t) + yij(t). Linearizing (3.3) at θ̄ (1) ij (t), we have dyij dt = � { h−1 [ H′(θi−1,j −θij)yi−1,j −H′(θi−1,j −θij)yij ] + h1 [ H′(θi+1,j −θij)yi+1,j −H′(θi+1,j −θij)yij ] + v1 [ H′(θi,j+1 −θij)yi,j+1 −H′(θi,j+1 −θij)yi,j ] + v−1 [ H′(θi,j−1 −θij)yi,j−1 −H′(θi,j−1 −θij)yi,j ]} . (3.5) In what follows and throughout the rest of the paper, we assume that indices k,l ∈ Z and we will omit this assumption over all the related summations. Since θi+k,j+l −θij = (k + l)Ψ, we can rewrite (3.5) CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 201 by using the notation wk,l for coupling strength as follows: dyij dt = � ∑ |k|+|l|=1 wk,l [ H′((k + l)Ψ)yi+k,j+l −H′((k + l)Ψ)yij ] = �   ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)yi+k,j+l − ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)yij   . (3.6) In view of the map f defined in Section 2.2, for any p ∈ {1, 2, · · · ,N2}, there exists a unique (i,j) ∈ {0, 1, 2, · · · ,N−1}×{0, 1, 2, · · · ,N−1} such that p = f(i,j). Hence we define xp = yij with p = f(i,j). Let x = (x1,x2, · · · ,xN2 )T and a(1) = ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ) = (h−1 + v−1)H ′(−Ψ) + (h1 + v1)H′(Ψ). So we can express (3.6) in terms of 1D indices, which is given by dx dt = �(M(1) −a(1)IN2 )x. (3.7) Here M(1) is a block circulant matrix and M(1) = bcirc ( h1H ′(Ψ)A1 + h−1H ′(−Ψ)A−1,v1H′(Ψ)IN, 0N, · · · , 0N,v−1H′(−Ψ)IN ) , where A1 and A−1 are defined in Section 2.1 and Ip is the p×p identity matrix for p ∈ N. Using Kronecker products, one can easily see that M(1) =IN ⊗ ( h1H ′(Ψ)A1 + h1H ′(−Ψ)A−1 ) + ( v1H ′(Ψ)A1 + v−1H ′(−Ψ)A−1 ) ⊗ IN. Then M(1) =h1H ′(Ψ)IN ⊗A1 + h−1H′(−Ψ)IN ⊗A−1 + v1H′(Ψ)A1 ⊗ IN + v−1H′(−Ψ)A−1 ⊗ IN = ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)(A1) l ⊗ (A1)k where the last equality is obtained by using wk,l notations for the coupling strengths, and (A1) −1 = A−1. It follows from Lemma 2.1, particularly, equations in (2.1), that M(1)ui ⊗uj = ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)ωli+kjui ⊗uj. This shows that the eigenvalues of M are {∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)ωli+kj : 0 ≤ i,j ≤ N − 1 } . So the eigenvalues of M −a1IN2 are λ (1) ij = ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)ωli+kj −a(1) = ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)(ωli+kj − 1), 0 ≤ i,j ≤ N − 1. Then <(λ(1)ij ) = ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ)[cos(2(li + kj)π/N) − 1] = −2 ∑ |k|+|l|=1 wk,lH ′((k + l)Ψ) sin2((li + kj)π/N). 202 J. CULP Note that sin2((li + kj)π/N) > 0 for all (k,l) ∈C1. So the cluster solution is stable if [w1,0H ′(Ψ) + w−1,0H ′(−Ψ)] sin2(jπ/N) + [w0,1H′(Ψ) + w0,−1H′(−Ψ)] sin2(iπ/N) > 0 for (i,j) 6= (0, 0). This leads to the following result. Theorem 3.3. Let N ≥ 3 and n be a factor of N with n ∈ NN . For our N ×N torus network model with von Neumann 1-neighborhood coupling, the n-cluster solution with constant phase difference Ψ is stable if [w1,0H ′(Ψ) + w−1,0H ′(−Ψ)] sin2(jπ/N) + [w0,1H′(Ψ) + w0,−1H′(−Ψ)] sin2(iπ/N) > 0 (3.8) for (i,j) 6= (0, 0), 0 ≤ i,j ≤ N − 1. Let Hodd(θ) = H(θ)−H(−θ) 2 , which is the odd part of the interconnection function H. Hence with the nearest neighbor coupling, we have the following result: Remark 3.1. (i) If Ψ = 0, it is the case of 1-cluster solution and the syncrhonization solution stable if H′odd(0) > 0, as H′odd(0) = H ′(0). (ii) If Ψ = π, it leads to an antiphase solution (2-cluster solution). For all the possible (k,l) in the nearest neighbor we have either H′(π) or H′(−π). By H′(π) = H′(−π) and H′odd(π) = H ′(π), the antiphase solution is stable is H′odd(π) > 0. (iii) If wkl ≡ w0 is constant for all k,l, the coupling is homogeneous and (3.8) becomes 2w0H ′ odd(Ψ)[sin 2(jπ/N) + sin2(iπ/N)] > 0 for 0 ≤ i,j ≤ N−1 with (i,j) 6= (0, 0). In this case, this cluster solution is stable if H′odd(Ψ) > 0. This precise stability condition on our 2D torus is in consistent with the result on the 1D ring [29]. 3.4. von Neumann 2-neighborhood Coupling: r = 2. Next, we consider the coupling is in the von Neumann 2-neighborhood, and we analyze the stability of the cluster solution with constant phase difference Ψ. dθij dt = Ω + � ∑ 0<|k|+|l|≤2 k,l∈Z wk,lH(θi+k,j+l −θij), 0 ≤ i,j ≤ N − 1. (3.9) Suppose that N ≥ 5 and n divides N. Consider an n-cluster solution with constant phase difference Ψ θ̄ (2) ij (t) = (Ω + �Ω (2) c )t + (i + j)Ψ, 0 ≤ i,j,≤ N − 1. (3.10) Here Ω(2)c = ∑ 0<|k|+|l|≤2 wk,lH((k + l)Ψ) which can be obtained by plugging (3.10) into (3.9). Let θij(t) = θ̄ (2) ij (t) + yij(t). Linearizing (3.9) at θ̄ (2) ij (t) to the first order gives dyij dt = � ∑ 0<|k|+|l|≤2 wk,l [ H′((k + l)Ψ)yi+k,j+l −H′((k + l)Ψ)yij ] = �   ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)yi+k,j+l − ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)yij   . CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 203 Let x = (x1,x2, · · · ,xN2 )T . Convert the indices from 2D to 1D. The above system can be written as dx dt = �(M(2) −a(2)IN2 )x, where M(2) = bcirc(M0,M1,M2, 0N, · · · , 0N,MN−2,MN−1) and a(2) = ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ) with M0 = w1,0H ′(Ψ)A1 + w2,0H ′(2Ψ)A21 + w−2,0H ′(−2Ψ)A2−1 + w−1,0H ′(−Ψ)A−1 M1 = w0,1H ′(Ψ)IN + w1,1H ′(2Ψ)A1 + w−1,1H ′(0)A−1 M2 = w0,2H ′(2Ψ)IN M3, . . . ,MN−3 = 0N MN−2 = w0,−2H ′(−2Ψ)IN MN−1 = w0,−1H ′(−Ψ)IN + w1,−1H′(0)A1 + w−1,−1H′(−2Ψ)A−1. Here 0N is the zero matrix of size N. Since any block circulant matrix can be expressed in terms of Kronecker products, M(2) = IN ⊗M0 + A1 ⊗M1 + A21 ⊗M2 + A N−2 1 ⊗MN−2 + A N−1 1 ⊗MN−1 = IN ⊗M0 + A1 ⊗M1 + A21 ⊗M2 + (A−1) 2 ⊗MN−2 + A−1 ⊗MN−1 = ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)(Al1 ⊗A k 1). By Lemma 2.1, (Al1 ⊗Ak1)(ui ⊗uj) = ωli+kjui ⊗uj and hence M(2)ui ⊗uj = ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)ωli+kjui ⊗uj, 0 ≤ i,j ≤ N − 1. It shows that the spectrum of M(2) is given by  ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)ωli+kj : 0 ≤ i,j ≤ N − 1   . So, the eigenvalues of M −a1IN2 are λ (2) ij = ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)ωli+kj −a(2) = ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)(ωli+kj − 1), 0 ≤ i,j ≤ N − 1. Hence, <(λ(2)ij ) = ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ)[cos(2(li + kj)π/N) − 1] = −2 ∑ 0<|k|+|l|≤2 wk,lH ′((k + l)Ψ) sin2((li + kj)π/N). 204 J. CULP Note that sin2((li + kj)π/N) > 0 for all (k,l) ∈ C2, when (i,j) 6= (0, 0) . Thus, the cluster solution is stable if [w2,0H ′(2Ψ) + w−2,0H ′(−2Ψ)] sin2(2jπ/N) + [w0,2H′(2Ψ) + w0,−2H′(−2Ψ)] sin2(2iπ/N) +[w1,1H ′(2Ψ) + w−1,−1H ′(−2Ψ)] sin2((i + j)π/N) + [w1,−1 + w−1,1]H′(0) sin2((i− j)π/N) +[w1,0H ′(Ψ) + w−1,0H ′(−Ψ)] sin2(jπ/N) + [w0,1H′(Ψ) + w0,−1H′(−Ψ)] sin2(iπ/N) > 0 for (i,j) 6= (0, 0) and 0 ≤ i,j ≤ N − 1. This leads to the following result. Theorem 3.4. Let N ≥ 5 and n be a factor of N with n ∈ NN . For our N ×N torus network model with von Neumann 2-neighborhood coupling, the n-cluster solution with constant phase difference Ψ is stable if [w2,0H ′(2Ψ) + w−2,0H ′(−2Ψ)] sin2(2jπ/N) + [w0,2H′(2Ψ) + w0,−2H′(−2Ψ)] sin2(2iπ/N) +[w1,1H ′(2Ψ) + w−1,−1H ′(−2Ψ)] sin2((i + j)π/N) + [w1,−1 + w−1,1]H′(0) sin2((i− j)π/N) +[w1,0H ′(Ψ) + w−1,0H ′(−Ψ)] sin2(jπ/N) + [w0,1H′(Ψ) + w0,−1H′(−Ψ)] sin2(iπ/N) > 0 (3.11) for (i,j) 6= (0, 0), 0 ≤ i,j ≤ N − 1. Thus, for the von Neumann 2-neighborhood coupling in our network model, we have the following result. Again H′odd(0) = H ′(0) and H′odd(π) = H ′(π). Remark 3.2. (i) If Ψ = 0, it is the case of 1-cluster solution and the synchronization solution stable if H′odd(0) > 0. (ii) If Ψ = π, it leads to an antiphase solution (2-cluster solution). For all the possible (k,l) in the nearest neighbor we have H′(π), H′(−π), H′(2π) and H′(0). So a stronger sufficient condition to guarantee that the antiphase solution is stable is H′odd(π) > 0 and H ′ odd(0) > 0. (iii) If wkl ≡ w0 is constant for all k,l, the coupling is homogeneous and (3.11) is equivalent to H′odd(Ψ)[sin 2(iπ/N) + sin2(jπ/N)] +H′odd(2Ψ)[sin 2(2iπ/N) + sin2(2jπ/N) + sin2((i+j)π/N)] + H′odd(0) sin 2((i − j)π/N) > 0 for 0 ≤ i,j ≤ N − 1 with (i,j) 6= (0, 0). Hence, in this case, a stronger sufficient condition to guarantee that this cluster solution is stable is that H′odd(0) > 0, H′odd(Ψ) > 0 and H ′ odd(2Ψ) > 0. 3.5. von Neumann r-neighborhood Coupling: r ∈ N. In this subsection, we extend the stability result to the general von Neumann r-neighborhood coupling for any given r ∈ N. By the similar method as that used in Section 3.5, it is straightforward to establish the result as follows. Theorem 3.5. Let N ≥ 2r + 1 and n be a factor of N with n ∈ NN . For our N ×N torus network model with von Neumann r-neighborhood coupling, the n-cluster solution with constant phase difference Ψ is stable if∑ 0<|k|+|l|≤r k,l∈Z wk,lH ′((k + l)Ψ) sin2((li + kj)π/N) > 0, (i,j) 6= (0, 0), 0 ≤ i,j ≤ N − 1. (3.12) Remark 3.3. (i) If Ψ = 0, it is the case of 1-cluster solution and the synchronization solution stable if H′odd(0) > 0. (ii) A stronger sufficient condition to guarantee that this cluster solution is stable is H′((k+l)Ψ) > 0 for (k,l) ∈Cr. CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 205 (iii) If wkl ≡ w0 is constant for all k,l, the coupling is homogeneous and (3.12) is equivalent to∑ 0<|k|+|l|≤r k,l∈Z H′odd((|k + l|)Ψ) sin 2((li + kj)π/N) > 0 for 0 ≤ i,j ≤ N − 1 with (i,j) 6= (0, 0). Hence, in this case, a stronger sufficient condition to guarantee that this cluster solution is stable is that H′odd((|k + l|)Ψ) > 0 for (k,l) ∈Cr. (iv) The stability result can be extended to an arbiturary coupling, where the set of coupling is denoted as C. Suppose that N is sufficiently large such that C ⊂ {0, 1, 2, · · · ,N − 1} × {0, 1, 2, · · · ,N − 1}. Theorem 3.6. Let n be a factor of N with n ∈ NN . For our N ×N torus network model with coupling C, the n-cluster solution with constant phase difference Ψ is stable if∑ (k,l)∈C wk,lH ′((k + l)Ψ) sin2((li + kj)π/N) > 0, (i,j) 6= (0, 0), 0 ≤ i,j ≤ N − 1. (3.13) 4. Numerical Results In this section we apply our analytical results to a network of N2 identical Morris-Lecar oscillators. We will run simulations for various values of N and various values of our coupling radius r to validate our analytical results. 4.1. Model and Parameter Analysis. We are using the dimensionless formulation of the Morris- Lecar model, [36, 7, 8], as given by the following system, where i ∈ [1, . . . ,N2],  v′i = Iapp −gCam∞(vi)(vi −vCa) −gKwi(vi −vk) −gL(vi −vL) −gsyn� N2∑ j=1 wijsj(vi −vsyn) w′i = φλ(vi)(w∞(vi) −wi) s′i = αh(vi)(1 −si) − si τs m∞(vi) = 1 2 ( 1 + tanh ( vi −V1 V2 )) w∞(vi) = 1 2 ( 1 + tanh ( vi −V3 V4 )) λ(vi) = cosh ( vi −V3 2V4 ) h(vi) = 5 exp(−(vi −vpre)/.1) Here we use the parameter values from [8], with the exception of Iapp,�,gsyn, all of which are given in Table 1. In the absence of coupling, each oscillator in the network has an exponentially asymptotically stable limit cycle with period ≈ 11.93 and frequency ≈ 0.53. 4.2. Calculating H and Hodd. The interaction function H can be calculated from the uncoupled single cell version of our Morris-Lecar model. The calculation of the interaction function H and it’s corresponding odd part, Hodd, were done in XPPAUT with the parameters given in Table 1. Data was exported from XPPAUT and imported to Python. See the text [12] for more information and tutorials on using XPPAUT. In Python, an univariate spline was then calculated for the function H and Hodd from the given data. From the univariate spline we can approximate the derivative of H and Hodd over one period. 206 J. CULP Parameter Name Value vCa Calcium equilibrium potential 1 vk Potassium equilibrium potential -0.7 vL Leak equilibrium potential -0.5 vsyn Synaptic reversal potential -0.625 vpre Pre-Synaptic membrane potential -0.1 gK Potassium ionic conductance 2 gL Leak ionic conductance 0.5 gCa Calcium potential conductance 1 φ Potassium rate constant 1/3 V1 Calcium activation potential -0.01 V2 Calcium reciprocal slope 0.15 V3 Potassium activation potential 0.1 V4 Potassium reciprocal slope 0.145 Iapp Applied current 0.123 gsyn Synaptic conductance 0.025 � Coupling strength 1 2r(r+1) α Synaptic activation constant 1.0 τs Synaptic decay constant 1.0 Table 1. Parameter Values 4.3. Numerical Simulations. Here we begin by considering the case when N = 5 with nearest and second nearest neighbor coupling. In each case, we would expect the existence of a synchronous solution, Ψ = 0, along with the 5-clusters Ψ = 2πk 5 , k = 1, 2, 3, 4. The values of H′(Ψ),H′(−Ψ) and H′odd(Ψ) for Ψ = 2πk 5 , k = 0, 1, 2, 3, 4 are given in Table 2. Stability analysis using some these values for the nearest and second nearest neighbor case will be provided below. k Ψ H′(Ψ) H′(−Ψ) H′odd(Ψ) 0 0 -0.5494149289091872 -0.5494149289091872 -0.6448920358548498 1 2π 5 -0.29701606374660317 0.18635559734182014 -0.05424507664335855 2 4π 5 0.12050167876589174 0.3801261880992415 0.24992749301022554 3 6π 5 0.3801261880992415 0.12050167876589174 0.24993264907886742 4 8π 5 0.18635559734182014 -0.2970160637466032 -0.05424962016118515 Table 2. H′(Ψ),H′(−Ψ),H′odd for Ψ = 0, 2π 5 , 4π 5 , 6π 5 , 8π 5 4.3.1. von Neumann 1-neighborhood Coupling. In Table 3, we see that in the case of homogeneous coupling our model predicts unstable and stable 5-clusters for various Ψ with N = 5. The firing groups consist of the following clusters of cells C1 = {1, 10, 14, 18, 22},C2 = {2, 6, 15, 19, 23},C3 = {3, 7, 11, 20, 24}, C4 = {4, 8, 12, 16, 25},C5 = {5, 9, 13, 17, 21} CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 207 Ψ Predicted Stability Verified Numerically 0 unstable Yes 2π 5 unstable Yes 4π 5 stable Yes 6π 5 stable Yes 8π 5 unstable Yes Table 3. Prediction Under Homogeneous Coupling, N = 5,r = 1,w = 1 In the case of Ψ = 4π 5 , stability was verified for the firing orders (C1,C4,C2,C5,C3), (C2,C5,C3,C1,C4), (C3,C1,C4,C2,C5), (C4,C2,C5,C3,C1), (C5,C3,C1,C4,C2) Recall, the stability condition for the nearest neighbor case with N = 5 is given by( w1,0H ′(Ψ) + w−1,0H ′(−Ψ) ) sin2 ( πj 5 ) + ( w0,1H ′(Ψ)v1 + w0,−1H ′(−Ψ) ) sin2 ( πi 5 ) > 0 (4.1) for 0 ≤ i,j ≤ 4, where i,j are both not 0. For i,j = 1, 2, 3, 4, we have that sin2(π 5 ) = sin2( 4π 5 ) and sin2( 2π 5 ) = sin2( 3π 5 ), so that a sufficient condition for stability is that( w1,0H ′(Ψ) + w−1,0H ′(−Ψ) ) > 0 and ( w0,1H ′(Ψ)v1 + w0,−1H ′(−Ψ) ) > 0 (4.2) From Table 2, we can see that for k = 2, 3 both H′(Ψ),H′(−Ψ) > 0. Therefore our simplified stability conditions will be positive for any choice of w±1,0,w0,±1. On the other hand, for k = 1, 4, H ′(Ψ) and H′(−Ψ) have alternate signs. In the case of homogeneous coupling the conditionality of the inequalities in (4.2) are completely dependent on the value of H′odd(Ψ) in Table 2, which is negative in either case so that we would expect an unstable solution. We are interested in cases where stability of our cluster solution could change given a change in the coupling strength, particularly from unstable to stable, so we will focus on cases of heterogenous coupling with k = 1, 4. Let us define the function F(w1,w−1) = w1H ′(Ψ) +w−1H ′(−Ψ) where w1 takes the place of w1,0,w0,1 and w−1 takes the place of w−1,0,w0,−1. The plot of F(w1,w−1) yields planar graphs which easily shows that there are many choices of w1,w−1 that would satisfy our inequalities in either case. We will focus on the case when k = 4, Ψ = 8π 5 . Here, we consider the cases for when w1,0 = w0,1 = 1, w−1,0 = w0,−1 = α, and vary α. A summary of values used in simulations are provided in Table 4, along with the smallest and largest value of the stability condition in (4.1). In each case in Table 4, a stable 5-cluster could not be verified. Case Smallest Value Largest Value α 1 ≈ 0.0131199 ≈ 0.068697 1 2 2 ≈ 0.038730 ≈ 0.202793 1 4 3 ≈ 0.057937 ≈ 0.3033663 1 16 4 ≈ 0.063940 ≈ 0.3347953 1 256 Table 4. Smallest and Largest Value of the Stability Condition, (4.1), for the Four Cases with Nearest Neighbor Coupling, r = 1,N = 5, Ψ = 8π 5 208 J. CULP 4.3.2. von Neumann 2-neighborhood Coupling. A similar analysis as done in the nearest neighbor cou- pling case will also be applied in our second nearest neighbor case by using our stability condition in (3.11), for N = 5, with 0 ≤ i,j ≤ 4, both not zero. In Table 5, we have that in the case of homogeneous coupling our model predicts unstable 5-clusters for all possible Ψ. We will consider then the case of heterogeneous coupling in which we show a change in stability of our cluster solutions. Ψ Predicted Stability Verified Numerically 0 unstable Yes 2π 5 unstable Yes 4π 5 unstable Yes 6π 5 unstable Yes 8π 5 unstable Yes Table 5. Prediction Under Homogeneous Coupling, N = 5,r = 2,w = 1 We will consider the case when Ψ = 4π 5 , so that H′(0) ≈−0.5491, H′(Ψ) ≈ 0.1204, H′(−Ψ) ≈ 0.3798, H′(2Ψ) ≈ 0.1862, H′(−2Ψ) ≈−0.2965 (4.3) From (3.11) and (4.3), we can see which coupling weights are associated with negative and positive values of H′. We can thus choose the corresponding coupling weights so that our stability condition, (3.11), is satisfied for all i,j and our cluster solution is predicted to be stable. One such configuration is to assign the coupling weights associated with second nearest neighbor coupling (wk,l, |k| + |l| = 2) the same value, α1 and the coupling weights associated with nearest neighbor coupling, (wk,l, |k|+ |l| = 1), the same value, α2. One such assignment, α1 = 1 16 ,α2 = 1, will allow us to satisfy the conditions in (3.11). In Figure ??, we have diagrams representing the change in sign of the real part of the eigenvalues with our assigned values versus the homogeneous coupling case. The corresponding raster plots for the beginning and end of these simulations are given in Figure 7, where darker colors indicate larger values. In Case 1, we can see that with strictly negative real parts of our eigenvalues where (3.11) is satisfied, we have a stable 5-cluster, likewise, in Case 2, where the conditions of (3.11) are not satisfied, we have an unstable cluster solution, as predicted. 4.3.3. Third Nearest Neighbor Coupling. Here will consider an 8 × 8 lattice of oscillators with third (r = 3) nearest neighbor coupling. In Table 6, we have a summary of predicted cluster solutions. Notice that in the case of Ψ = π 2 ,π, 3π 2 we were unable to verify the predicted stability. In Figure 4, we have a diagram giving the real part of the eigenvalues and the beginning and end of a simulation in the case with Ψ = 3π 2 , the case for Ψ = π 2 is similar. Here we see that the simulation appears to produce a stable cluster as a result, in contrast to the unstable prediction. We also note that almost all the real parts of the eigenvalues are negative and much larger in magnitude than the two eigenvalues with positive real part, ≈ 0.1322. In Figure 5, we see the real part of the eigenvalues and the beginning and end of the simulation in the unverified unstable case with homogeneous coupling and Ψ = π. Here, the simulation appears to be stable although there is a mixture of positive and negative real parts of the eigenvalues, with the some large negative eigenvalues. In Figure 6, we consider the case with heterogenous coupling with values wk,l = 1/4 when |k| + |l| = 3; wk,l = 1 when |k| + |l| = 2; wk,l = 1/4 when |k| + |l| = 1. In this case, where most eigenvalues are positive and large in magnitude, we have an unstable solution, as would have been predicted. CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 209 (a) Fitness function (b) Anti-predation response level α Figure 3. The fitness function F1(α,v1) and F2(α,v2) and convergent dynamics of anti-predation response level α(t) on the two patches which are not connected by dis- persals. Ψ Predicted Stability Verified Numerically 0 unstable Yes π 4 unstable Yes π 2 unstable No 3π 4 unstable Yes π unstable No 5π 4 unstable Yes 3π 2 unstable No 7π 4 unstable Yes Table 6. Prediction Under Homogeneous Coupling, N = 8,r = 3,w = 1 In Figure 8, we show 5 firing diagrams across one period of our simulation with the coupling values given in wk,l = 1/16 when |k| + |l| = 3; wk,l = 1/8 when |k| + |l| = 2; wk,l = 1 when |k| + |l| = 1. The cell numberings are given within each cell in Figure 8. Cells that fire together share the same color in each diagram, and darker colors imply larger values. We expect 4 subgroups of 16 oscillators each, were the predicted grouping are given in (4.4). From Figure 8, there appears to be a traveling wave firing pattern which agrees with the firing order given in (4.4). C1 ={1, 5, 12, 16, 19, 23, 26, 30, 33, 37, 44, 48, 51, 55, 58, 62} C4 ={4, 8, 11, 15, 18, 22, 25, 29, 36, 40, 43, 47, 50, 54, 57, 61} C3 ={3, 7, 10, 14, 17, 21, 28, 32, 35, 39, 42, 46, 49, 53, 60, 64} C2 ={2, 6, 9, 13, 20, 24, 27, 31, 34, 38, 41, 45, 52, 56, 59, 63} (4.4) 4.3.4. Cluster Solutions in a Larger Network. Here we consider a 18 × 18 lattice of 324 oscillators. Under the assumption of homogeneous coupling weights, we will begin with the case of nearest neighbor 210 J. CULP (a) Real Part of the Eigenvalues (b) Beginning and End of Simulation Figure 4. Unverified Third Nearest Neighbor Homogenous Coupling Case with N = 8, Ψ = 3π 2 (a) Real Part of the Eigenvalues (b) Beginning and End of Simulation Figure 5. Unverified Third Nearest Neighbor Homogenous Coupling Case with N = 8, Ψ = π coupling and increase our neighborhood radius, r, and study the change in stability. Here we find that for the case of nearest neighbor coupling the only predicted stable cluster solutions are for Ψ = 5π 9 , 2π 3 , 7π 9 , 8π 9 ,π, 10π 9 , 11π 9 , 4π 3 , 13π 9 (4.5) When we move to the case with second nearest neighbor homogeneous coupling, we have that there are no predicted stable cluster solutions, but when we move to the third nearest neighbor case we have that Ψ = 4π 9 , 14π 9 are the only stable cluster solutions. Therefore, for homogeneous coupling weights, there are no cluster solutions that are predicted to stable for the nearest neighbor and remain stable as our neighborhood radius increases. CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 211 (a) Real Part of the Eigenvalues (b) Beginning and End of Simulation Figure 6. Unverified Third Nearest Neighbor Heterogeneous Coupling Case with N = 8, Ψ = π Under the values of Ψ in (4.5), with nearest neighbor coupling, we find that for second nearest neighbor coupling with the coupling weights; wk,l = 1/4 when |k| + |l| = 2; wk,l = 1 when |k| + |l| = 1, that only Ψ = 5π 9 , 2π 3 , 4π 3 , 13π 9 (4.6) are predicted to be stable. When we furthermore extend our analysis to third nearest neighbor coupling with the coupling weights; wk,l = 1/8 when |k| + |l| = 3; wk,l = 1/4 when |k| + |l| = 2; wk,l = 1 when |k| + |l| = 1 for (4.6), we find that stability is only predicted in the case for Ψ = 5π 9 , 13π 9 . If we change our coupling weights to be 1, 1 2 for 1, 1 4 in the case of second nearest neighbor coupling, then there are no predicted stable cluster solutions. Likewise, if we change the coupling weights in the third nearest neighbor case to 1, 1 2 , 1 4 for 1, 1 4 , 1 8 , then there are no predicted stable clusters. Given in Figure 9 are the plots of the stable cluster solutions, with Ψ = 5π 9 , for the third nearest neighbor heterogeneous coupling case. The plots for the first and second nearest neighbor coupling cases are similar and will be omitted. 4.3.5. Numerical Results for the Case with Ψh 6= Ψv. Here we examine the firing patterns for a 16×16 lattice of oscillators with various values of Ψh 6= Ψv and third nearest neighbor coupling. In Figure 10, we use the coupling weights; wk,l = 1/8 when |k|+|l| = 3; wk,l = 1/4 when |k|+|l| = 2; wk,l = 1 when |k| + |l| = 1, and kh = 4,kv = 12. Using these coupling weights and choice of kh,kv, we predict 4 stable clusters of 64 oscillators each, where that Ψh = π 2 , Ψv = 3π 2 . In Figure 11, we plot our stable cluster solutions over one period (17 frames, where frame 0 and 16 are equivalent) with the same coupling weights as above and kh = 9,kv = 8. Here we predict 16 stable clusters of 16 oscillators each, where that Ψh = 9π 8 , Ψv = π. In Figure 10, we have a traveling wave pattern that appears to move through from the upper left to bottom right of the plot. This pattern is opposed to the traveling wave pattern in Figure 8, in which the wave appears to travel from the upper right to bottom left of the plot. In Figure 11, there appears to be a traveling vertical checkerboard pattern. 212 J. CULP (a) Stable Solution (b) Unstable Solution Figure 7. Effect of Varying Coupling Values on Stability of Cluster Solutions in Sec- ond Nearest Neighbor Coupling Case 5. Discussion In this work we considered a 2D lattice of N ×N weakly coupled oscillators with period boundary conditions and von Neumann neighborhood r coupling. For our analytical results we derived conditions CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 213 Figure 8. Firing Sequence Over One Period, T for the existence and stability of cluster solutions in our network under various coupling situations, while taking advantage of known results, involving circulant matrices and the Kronecker product, so that our stability analysis can be easily performed. We would like to extend our analytical results to include existence and stability analysis for an arbitrary N1 × N2,N1 6= N2 lattice of oscillators with periodic boundary conditions. This extension would dramatically increase the generality of the model, although it is not clear to what extent our representation techniques could be applied. Additionally, in the case for non-square lattices, it seems likely from our work with Ψh 6= Ψv, that the existence of any cluster solutions would require a condition of the kind in which N1,N2 are multiples of each other. Expanding our 2D network topology results to 3D, where a von Neumann neighborhood is also defined, is also possible. There are additional extensions to our model that could also be investigated. It will be interesting to consider a modification of our model to include an explicit time delay. This modification has biology significance and been studied in [7] for a ring topology. There are also certainly a large diversity of interesting coupling structures possible in larger networks, such as a randomly connected network, that we could incorporate into our network topology. Similarly as above, it would be uncertain as to which our convenient representation technique would be limited under these coupling assumptions. In Table 6, we were not able to numerically verify the case with Ψ = π 2 ,π, 3π 2 . Recalling, Figure 4 in our numerical section, in the case of Ψ = π 2 , 3π 2 , we have a majority negative real part of the eigenvalues, which are large in magnitude compared to the two small positive eigenvalues. It is possible 214 J. CULP Figure 9. 18 × 18 Lattice of Coupled Oscillators with Third Nearest Neighbor Coupling Figure 10. Firing Sequence Over One Period, T, for N = 16,r = 3, Ψh = π 2 , Ψv = 3π 2 that a longer simulation time might be needed to show the correct stability. We should also note that since our positive eigenvalues are close to zero, it might be the case that there could be numerical errors in the approximation of our interaction function H or it’s derivative H′ that are resulting in CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 215 Figure 11. Firing Sequence Over One Period, for N = 16,r = 3, Ψh = 9π 8 , Ψv = π incorrectly giving the positive value on these small eigenvalues. Additional numerical errors could have also occurred in our numerical ODE solver. In the case with N = 5 and nearest neighbor coupling we were not successful in showing a change of stability from adjusting our coupling weights. It is possible that the inability to show a change in 216 J. CULP stability is due to the same numerical issues as previously mentioned for our work in the case with Ψ = π 2 ,π, 3π 2 . Additional investigation could still yield a change in stability, possibly by using more convenient values of Iapp and τs, but we were unable to find these values. In our numerical examples, we consider relatively small values of N. Therefore, additional simulations with large values of N should be run. Here the ability to produce more diverse firing patterns, such as the rotating spiral wave [33, 3], on the 2D torus with our periodic boundary conditions, could be studied. References [1] P. Ashwin, S. Coombes, and R. Nicks, Mathematical frameworks for oscillatory network dynamics in neuroscience, The Journal of Mathematical Neuroscience 6 (2016), no. 1, 1–92. [2] P. Ashwin and J.W. Swift, The dynamics of n weakly coupled identical oscillators, Journal of Nonlinear Science 2 (1992), no. 1, 69–108. [3] J. Bramburger, Rotating wave solutions to lattice dynamical systems i: the anti-continuum limit, Journal of Dynamics and Differential Equations 31 (2019), no. 1, 469–498. [4] , Stability of infinite systems of coupled oscillators via random walks on weighted graphs, Transactions of the American Mathematical Society 372 (2019), no. 2, 1159–1192. [5] E. Brown, P. Holmes, and J. Moehlis, Globally coupled oscillator networks, Perspectives and Problems in Nolinear Science, Springer, 2003, pp. 183–215. [6] G. Buzsáki and A. Draguhn, Neuronal oscillations in cortical networks, science 304 (2004), no. 5679, 1926–1929. [7] S.A. Campbell and I. Kobelevskiy, Phase models and oscillators with time delayed coupling, Discrete & Continuous Dynamical Systems-A 32 (2012), no. 8, 2653. [8] S.A. Campbell and Z. Wang, Phase models and clustering in networks of oscillators with delayed coupling, Physica D: Nonlinear Phenomena 363 (2018), 44–55. [9] S. Coombes and J.R. Terry, The dynamics of neurological disease: integrating computational, experimental and clinical neuroscience, 2012. [10] P.J. Davis, Circulant matrices, Wiley, 1979. [11] G.B. Ermentrout, The behavior of rings of coupled oscillators, Journal of mathematical biology 23 (1985), no. 1, 55–74. [12] , Simulating, analyzing, and animating dynamical systems: a guide to xppaut for researchers and students, SIAM, 2002. [13] G.B. Ermentrout and N. Kopell, Frequency plateaus in a chain of weakly coupled oscillators, i., SIAM journal on Mathematical Analysis 15 (1984), no. 2, 215–237. [14] J. Fell and N. Axmacher, The role of phase synchronization in memory processes, Nature reviews neuroscience 12 (2011), no. 2, 105–118. [15] M. Golubitsky, I. Stewart, and D.G. Schaeffer, Singularities and groups in bifurcation theory: Volume ii, vol. 69, Springer Science & Business Media, 2012. [16] H. Haken, J.A Kelso, and H. Bunz, A theoretical model of phase transitions in human hand movements, Biological cybernetics 51 (1985), no. 5, 347–356. [17] D. Hansel, G. Mato, and C. Meunier, Clustering and slow switching in globally coupled phase oscillators, Physical Review E 48 (1993), no. 5, 3470. [18] , Phase dynamics for weakly coupled hodgkin-huxley neurons, EPL (Europhysics Letters) 23 (1993), no. 5, 367. [19] F.C. Hoppensteadt and E.M. Izhikevich, Weakly connected neural networks, vol. 126, Springer Science & Business Media, 2012. [20] E.M. Izhikevich and Y. Kuramoto, Weakly coupled oscillators, Encyclopedia of mathematical physics 5 (2006), 448. [21] F.G. Kazanci and G.B. Ermentrout, Wave formation through the interactions between clustered states and local coupling in arrays of neural oscillators, SIAM Journal on Applied Dynamical Systems 7 (2008), no. 2, 491–509. [22] J.A. Kelso, Dynamic patterns: The self-organization of brain and behavior, MIT press, 1995. [23] N. Kopell and G.B. Ermentrout, Coupled oscillators and the design of central pattern generators, Mathematical biosciences 90 (1988), no. 1-2, 87–109. [24] , Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators, Handbook of dynamical systems 2 (2002), 3–54. [25] Y. Kuramoto, Chemical oscillations, waves, and turbulence, Courier Corporation, 2003. CLUSTER SOLUTIONS ON A 2D SQUARE TORUS 217 [26] C.R. Laing, Rotating waves in rings of coupled oscillators, Dynamics and stability of systems 13 (1998), no. 4, 305–318. [27] G. Laurent, M. Stopfer, R.W. Friedrich, M.I. Rabinovich, A. Volkovskii, and H.D. Abarbanel, Odor encoding as an active, dynamical process: experiments, computation, and theory, Annual review of neuroscience 24 (2001), no. 1, 263–297. [28] J.G. Mancilla, T.J. Lewis, D.J. Pinto, J Rinzel, and B.W. Connors, Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex, Journal of Neuroscience 27 (2007), no. 8, 2058–2073. [29] J. Miller, H. Ryu, Z. Teymuroglu, X. Wang, V. Booth, and S.A. Campbell, Clustering in inhibitory neural networks with nearest neighbor coupling, Applications of dynamical systems in biology and medicine, Springer, 2015, pp. 99– 121. [30] R.E. Mirollo and S. Strogatz, Synchronization of pulse-coupled biological oscillators, SIAM Journal on Applied Mathematics 50 (1990), no. 6, 1645–1662. [31] J.C. Neu, Coupled chemical oscillators, SIAM Journal on Applied Mathematics 37 (1979), no. 2, 307–315. [32] K. Okuda, Variety and generality of clustering in globally coupled oscillators, Physica D: Nonlinear Phenomena 63 (1993), no. 3-4, 424–436. [33] J.E. Paullet and G.B. Ermentrout, Stable rotating waves in two-dimensional discrete active media, SIAM Journal on Applied Mathematics 54 (1994), no. 6, 1720–1744. [34] C.S. Peskin and Courant Institute of Mathematical Sciences, Mathematical aspects of heart physiology, Courant Institute Lecture Notes, Courant Institute of Mathematical Sciences, New York University, 1975. [35] L. Ren and G.B. Ermentrout, Monotonicity of phaselocked solutions in chains and arrays of nearest-neighbor coupled oscillators, SIAM journal on mathematical analysis 29 (1998), no. 1, 208–234. [36] J. Rinzel and G.B. Ermentrout, Analysis of neural excitability and oscillations, Methods in neuronal modeling 2 (1998), 251–292. [37] H. Sakaguchi and Y. Kuramoto, A soluble active rotater model showing phase transitions via mutual entertainment, Progress of Theoretical Physics 76 (1986), no. 3, 576–581. [38] F. Saraga, L. Ng, and F.K. Skinner, Distal gap junctions and active dendrites can tune network dynamics, Journal of neurophysiology 95 (2006), no. 3, 1669–1682. [39] M.A. Schwemmer and T.J. Lewis, The theory of weakly coupled oscillators, Phase response curves in neuroscience, Springer, 2012, pp. 3–31. [40] W. Singer, Synchronization of cortical activity and its putative role in information processing and learning, Annual review of physiology 55 (1993), no. 1, 349–374. [41] P.S.G. Stein, D.G Stuart, S. Grillner, and A.I. Selverston, Neurons, networks, and motor behavior, MIT press, 1999. [42] A. Takamatsu, T. Fujii, and I. Endo, Time delay effect in a living coupled oscillator system with the plasmodium of physarum polycephalum, Physical Review Letters 85 (2000), no. 9, 2026. [43] M.S. Titcombe, R. Edwards, and A. Beuter, Mathematical modelling of parkinsonian tremor., Nonlinear Studies 11 (2004), no. 3, 363–384. [44] P. Uhlhaas, G. Pipa, B. Lima, L. Melloni, S. Neuenschwander, D. Nikolić, and W. Singer, Neural synchrony in cortical networks: history, concept and current status, Frontiers in integrative neuroscience 3 (2009), 17. [45] J.L.P. Velazquez, Brain research: a perspective from the coupled oscillators field, NeuroQuantology 4 (2006), no. 2, 155–165. [46] S.S. Wang and H.G. Winful, Dynamics of phase-locked semiconductor laser arrays, Applied physics letters 52 (1988), no. 21, 1774–1776. [47] X.J. Wang, Neurophysiological and computational principles of cortical rhythms in cognition, Physiological reviews 90 (2010), no. 3, 1195–1268. [48] M. Wehr and G. Laurent, Odour encoding by temporal sequences of firing in oscillating neural assemblies, Nature 384 (1996), no. 6605, 162–166. [49] H.G. Winful and S.S Wang, Stability of phase locking in coupled semiconductor laser arrays, Applied physics letters 53 (1988), no. 20, 1894–1896. 218 J. CULP Department of Cell Biology & Anatomy, University of Calgary, Calgary, Alberta, Canada Current address: same E-mail address: jordan.culp@ucalgary.ca