CLUSTER SOLUTIONS IN NETWORKS OF WEAKLY COUPLED OSCILLATORS ON A 2D SQUARE TORUS

. 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 diﬀerences (Ψ 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 diﬀerence (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.


Introduction
An average neuron forms and receives one to ten thousand synaptic connections for sending and receiving information. Since there are at least 10 11 neurons in the human brain, there are thus 10 14 ∼ 10 15 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 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.

Notations
In this section we provide some preliminary mathematical theory and notation to be used in this work.

2.1.
Matrices. Define a circulant matrix generated by constants c 1 , c 2 , c 3 , · · · , c N −1 , c N as Lemma 2.1. Let p, q ∈ {−(N − 1), · · · , −2, −1, 0, 1, 2, · · · , N − 1} be given. Then Proof. By the commutative property of the Kronecker product, This result shows that, for any p, q ∈ {−(N −1), · · · , −2, −1, 0, 1, 2, · · · , N −1}, A p 1 ⊗A q 1 has the same eigenvectors {u i ⊗u j : 0 ≤ i, j ≤ N −1} and the corresponding eigenvalues are {ω pi+qj : 0 ≤ i, j ≤ N −1}. In particular, 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, . . . , N 2 . Here, for 0 ≤ i, j ≤ N − 1, we use the mapping f : Clearly f is a one-to-one and onto mapping from {0, 1, 2, · · · , N − 1} × {0, 1, 2, · · · , N − 1} to {1, 2, . . . , N 2 }. 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. 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 r th 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 1neighborhood 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 r 2 + (r + 1) 2 − 1.

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.
Here y ij ∈ R m , 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),(ĩ,j) ) is the coupling matrix, where as c (i,j),(ĩ,j) ≥ 0 denotes the coupling between the oscillators y ij and yĩj. Furthermore, we assume that for an isolated oscillator, the solution to dyij dt = F (y ij ) exhibits an exponentially asymptotically T −periodic  Figure 2. The von Neumann 2-neighborhood of node P stable solution, denoted byX ij (t) on the domain 0 ≤ t ≤ T . Here, the frequency of the ij th 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, where H is the interaction function with period 2π that satisfies Here Z is known as the phase response curve of the isolated oscillator, which is the unique periodic solution of the linearized adjoint system: 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 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 ): 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 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].

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 C r := {(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: where the coupling is assumed to be relative-location dependent, and w k,l = c (i+k,j+l),(i,j) is the coupling between the ij th oscillator and (i + k, j + l) th oscillator for all (i, j).

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 ij th oscillator is 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 Then the solution with phase difference (Ψ h , Ψ v ) gives us a cluster solution that is of period n h and n v along the horizontal and vertical directions, respectively, where Ψ h = 2m h π/n h and Ψ v = 2m v π/n v . It is clear that the period of this cluster solution is n := lcm(n h , n v ). This solution defines an n-cluster solution with phase difference (Ψ h , Ψ v ). Let It leads to the following existence result.
(iv) Assume that both n v and n h are factors of N with n h , n v ∈ N N . Define n = lcm(n h , n v ). There exists an n-cluster solution with phase difference . 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. (i) There exists a synchronization solution with Ψ = 0.
(ii) Assume that n is a factor of N with n ∈ N N . 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.

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 Suppose that N ≥ 3 and n divides N . Consider an n-cluster solution with constant phase difference Ψθ 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) by using the notation w k,l for coupling strength as follows: In view of the map f defined in Section 2.2, for any So we can express (3.6) in terms of 1D indices, which is given by Here M (1) is a block circulant matrix and where A 1 and A −1 are defined in Section 2.1 and I p is the p × p identity matrix for p ∈ N.
Using Kronecker products, one can easily see that Then where the last equality is obtained by using w k,l notations for the coupling strengths, and (A 1 ) −1 = A −1 . It follows from Lemma 2.1, particularly, equations in (2.1), that This shows that the eigenvalues of M are , which is the odd part of the interconnection function H. Hence with the nearest neighbor coupling, we have the following result: 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 Ψ.
Suppose that N ≥ 5 and n divides N . Consider an n-cluster solution with constant phase difference Ψθ Here which can be obtained by plugging (3.10) into (3.9). Let θ ij (t) =θ (2) ij (t) + y ij (t). Linearizing (3.9) atθ Let x = (x 1 , x 2 , · · · , x N 2 ) T . Convert the indices from 2D to 1D. The above system can be written as and Here 0 N is the zero matrix of size N .
Since any block circulant matrix can be expressed in terms of Kronecker products, It shows that the spectrum of M (2) is given by 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 (π). 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.
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) ∈ C r . (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 ∈ N N . For our N × N torus network model with coupling C, the n-cluster solution with constant phase difference Ψ is stable if

Numerical Results
In this section we apply our analytical results to a network of N 2 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.

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, . . . , N 2 ], Here we use the parameter values from [8], with the exception of I app , , g syn , 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 H odd . 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, H odd , 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 H odd from the given data. From the univariate spline we can approximate the derivative of H and H odd over one period.   Table 2. Stability analysis using some these values for the nearest and second nearest neighbor case will be provided below.  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 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 (C 1 , C 4 , C 2 , C 5 , C 3 ), (C 2 , C 5 , C 3 , C 1 , C 4 ), (C 3 , C 1 , C 4 , C 2 , C 5 ), Recall, the stability condition for the nearest neighbor case with N = 5 is given by for 0 ≤ i, j ≤ 4, where i, j are both not 0. For i, j = 1, 2, 3, 4, we have that sin 2 ( π 5 ) = sin 2 ( 4π 5 ) and sin 2 ( 2π 5 ) = sin 2 ( 3π 5 ), so that a sufficient condition for stability is that 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 , w 0,±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 (w 1 , w −1 ) = w 1 H (Ψ) + w −1 H (−Ψ) where w 1 takes the place of w 1,0 , w 0,1 and w −1 takes the place of w −1,0 , w 0,−1 . The plot of F (w 1 , w −1 ) yields planar graphs which easily shows that there are many choices of w 1 , 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 w 1,0 = w 0,1 = 1, w −1,0 = w 0,−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.  A similar analysis as done in the nearest neighbor coupling 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. unstable Yes Table 5. Prediction Under Homogeneous Coupling, N = 5, r = 2, w = 1 We will consider the case when Ψ = 4π 5 , so that 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 (w k,l , |k| + |l| = 2) the same value, α 1 and the coupling weights associated with nearest neighbor coupling, (w k,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.

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 w k,l = 1/4 when |k| + |l| = 3; w k,l = 1 when |k| + |l| = 2; w k,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.   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 w k,l = 1/16 when |k| + |l| = 3; w k,l = 1/8 when |k| + |l| = 2; w k,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     Under the values of Ψ in (4.5), with nearest neighbor coupling, we find that for second nearest neighbor coupling with the coupling weights; w k,l = 1/4 when |k| + |l| = 2; w k,l = 1 when |k| + |l| = 1, that only are predicted to be stable. When we furthermore extend our analysis to third nearest neighbor coupling with the coupling weights; w k,l = 1/8 when |k| + |l| = 3; w k,l = 1/4 when |k| + |l| = 2; w k,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 = Ψ v . Here we examine the firing patterns for a 16 × 16 lattice of oscillators with various values of Ψ h = Ψ v and third nearest neighbor coupling.
In Figure 10, we use the coupling weights; w k,l = 1/8 when |k| + |l| = 3; w k,l = 1/4 when |k| + |l| = 2; w k,l = 1 when |k| + |l| = 1, and k h = 4, k v = 12. Using these coupling weights and choice of k h , k v , 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 k h = 9, k v = 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.

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 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 N 1 × N 2 , N 1 = N 2 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 = Ψ v , that the existence of any cluster solutions would require a condition of the kind in which N 1 , N 2 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   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 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 I app 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.
Department of Cell Biology & Anatomy, University of Calgary, Calgary, Alberta, Canada Current address: same E-mail address: jordan.culp@ucalgary.ca