# Plastic systemic inhibition controls amplitude while allowing phase pattern
in a stochastic neural field model

Conor L. Morrison  Department of Statistics, University of British Columbia,
Vancouver, BC Canada V6T 1Z4  Priscilla E. Greenwood  Department of
Mathematics, University of British Columbia, Vancouver, BC Canada V6T 1Z2
Lawrence M. Ward  Department of Psychology and Djavad Mowafaghian Centre for
Brain Health, 2136 West Mall, University of British Columbia, Vancouver, BC
Canada V6T 1Z4

###### Abstract

Oscillatory phase pattern formation and amplitude control for a linearized
stochastic neuron field model was investigated by simulating coupled
stochastic processes defined by stochastic differential equations. It was
found, for several choices of parameters, that pattern formation in the phases
of these processes occurred if and only if the amplitudes were allowed to grow
large. Stimulated by recent work on homeostatic inhibitory plasticity, we
introduced static and plastic (adaptive) systemic inhibitory mechanisms to
keep the amplitudes stochastically bounded in subsequent simulations. The
systems with static systemic inhibition exhibited bounded amplitudes but no
sustained phase patterns, whereas the systems with plastic systemic inhibition
exhibited both bounded amplitudes and sustained phase patterns. These results
demonstrate that plastic inhibitory mechanisms in neural field models can
stochastically control amplitudes while allowing patterns of phase
synchronization to develop. Similar mechanisms of plastic systemic inhibition
could play a role in regulating oscillatory functioning in the brain.

###### pacs:

##  I Introduction

Mathematical models that exhibit oscillations have played a key role in
modeling neural phenomena. Such models include models of individual neurons —
e.g., leaky integrate-and-fire neurons Sacerdote and Giraudo (2013), the
Izhikevich neuron Izhikevich (2003), or the Hodgkin-Huxley neuron Hodgkin and
Huxley (1952) — and models of neuron populations, such as the Wilson-Cowan
(W-C) model Wilson and Cowan (1972) or neural field models such as those
studied by Faugeras and Inglis Faugeras and Inglis (2015). Perhaps the most
influential model of neuron populations has been that of Wilson and Cowan
Wilson and Cowan (1972). They proposed a nonlinear rate model for two
interacting populations of neurons:

$\begin{split}&dV_{E}(t)=\bigg{[}-V_{E}(t)+g\Big{[}a_{E}\Big{(}S_{EE}V_{E}(t)\\\
&-S_{EI}V_{I}(t)-\theta_{E}+P_{E}(t)\Big{)}\Big{]}\bigg{]}dt,\end{split}$ (1)
$\begin{split}&dV_{I}(t)=\bigg{[}-V_{I}(t)+g\Big{[}a_{I}\Big{(}-S_{II}V_{I}(t)\\\
&+S_{IE}V_{E}(t)-\theta_{I}\Big{)}\Big{]}\bigg{]}dt,\end{split}$ (2)

where $V_{E}(t),V_{I}(t)$ are voltages of excitatory and inhibitory neuron
populations, respectively, $S_{EE},S_{II},S_{EI},S_{IE}$ are synaptic
efficacies that control their interactions, $g$ is a sigmoid threshold
function, $a_{E},a_{I},\theta_{E},\theta_{I}$ are constants, and $P_{E}(t)$ is
input current. In the W-C model (1), (2) and in most neural field models a
nonlinear transform $g$ is introduced in order to bound the amplitudes, where
$g$ is typically the logistic function or some other sigmoid function. This
allows the formation of limit cycles that resemble oscillations recorded from
neural systems.

Nonlinear models of neuron populations like the W-C model can be linearized
and centered at a relevant fixed point, resulting in systems like that studied
by Kang et al. Kang  _et al._ (2010) and by Greenwood et al. Greenwood  _et
al._ (2015). The latter system, which we will employ in what follows, is
governed by stochastic differential equations for mean voltages of excitatory
and inhibitory neuron populations — $E$ and $I$, respectively — given by,

$\begin{split}&\tau_{E}dE(t)=(-E(t)+S_{EE}E(t)\\\
&-S_{EI}I(t))dt+\sigma_{E}dW_{E}(t),\end{split}$ (3)

and,

$\begin{split}&\tau_{I}dI(t)=(-I(t)-S_{II}E(t)\\\
&+S_{IE}E(t))dt+\sigma_{I}dW_{I}(t),\end{split}$ (4)

or equivalently by the matrix equation,

$\begin{split}d\mathbb{V}=[E(t),I(t)]^{\top}=-\mathbb{A}\mathbb{V}(t)dt+\mathbb{N}d\mathbb{W},\end{split}$
(5)

where,

$\begin{split}-\mathbb{A}=\begin{pmatrix}(-1+S_{EE})\tau_{E}^{-1}&-S_{EI}\tau_{E}^{-1}\\\
S_{IE}\tau_{I}^{-1}&(-1-S_{II})\tau_{I}^{-1}\end{pmatrix},\end{split}$ (6)

where $d\mathbb{W}=(dW^{E},dW^{I})^{\top}$ are independent standard Brownian
motions, and where,

$\begin{split}\mathbb{N}=\begin{pmatrix}\sigma_{E}\tau_{E}^{-1}&0\\\
0&\sigma_{I}\tau_{I}^{-1}\end{pmatrix}.\end{split}$ (7)

Here $\tau_{E}$ and $\tau_{I}$ represent time constants, $S_{II}$, $S_{IE}$,
$S_{EI}$, and $S_{EE}$ represent mean synaptic connection efficacies, and
$\sigma_{I}$ and $\sigma_{E}$ are amplitudes. We will refer to such a system
as an EI-pair system.

The model (5) is a form of W-C model but with the addition of time constants
and stochastic terms and the omission of several other constants, the input
current, and, most importantly for our purposes, the amplitude-bounding
sigmoid function $g$ (as in (1), (2)). Without the amplitude-bounding
function, the model (5) does not exhibit a noisy limit cycle, as a stochastic
form of the W-C model does in Wallace  _et al._ (2011). In the absence of a
limit cycle, however, oscillations can be sustained owing to the stochasticity
in the system. These oscillations are referred to as ‘quasi-cycles’ Bressloff
(2010) and occur when the matrix $-\mathbb{A}$, as in (5), has complex
eigenvalues $-\lambda\pm\omega i$, with $0<\lambda$, and the matrix
$\mathbb{N}\neq 0$. With $\mathbb{N}=0$ the system damps to the fixed point
because the real part of the eigenvalues, $-\lambda$, is negative. Adding
moderate noise sustains noisy oscillations at the ‘natural’ frequency
$\omega$. Quasi-cycles in Wilson-Cowan-type population models were first
described by Bressloff Bressloff (2010), and those in systems such as (5) were
investigated mathematically by Baxendale and Greenwood Baxendale and Greenwood
(2011).

Recent work has investigated _coupled_ systems of quasi-cycle oscillatory EI-
pairs, a single pair of which whose stochasticity has been shown to reproduce
qualitative behaviour observed in neural systems, such as gamma bursts
Greenwood  _et al._ (2015). In Greenwood  _et al._ (2016), a Kuramoto-type
coupling scheme was applied to the phases
($\theta(t)=\arctan(\frac{I(t)}{E(t)})$) and amplitudes
($Z(t)=||\mathbb{V}(t)||$) of EI-pairs under an approximation given by
Baxendale and Greenwood in Baxendale and Greenwood (2011). A linear type
coupling scheme was used, with coupling strengths given by a coupling matrix
$\mathbb{C}=(\mathbb{C}_{k,j})$, where $\mathbb{C}_{k,j}$ represents the
coupling effect of neuron $j$ to neuron $k$. Phase synchronization was
demonstrated to occur at a critical value of $||\mathbb{C}||$ – in line with
the work done by Kuramoto Kuramoto (1984). Similar work was done in Greenwood
and Ward (2019) on a lattice of discrete-Mexican-hat coupled EI-pairs,
demonstrating patterns of ordered phases in both space and time that appeared
quickly and then evolved slowly.

In some papers, such as Greenwood  _et al._ (2016) and Greenwood and Ward
(2019), the amplitudes of the coupled quasi-cycle oscillators (which
correspond to their peak-to-trough voltages) appeared to grow without bound
for some parameter regimes in the absence of a bounding function $g$. Although
amplitudes of limit cycles were not discussed by Kuramoto in Kuramoto (1984),
it is important to consider the behaviour of the amplitudes when appraising
phase pattern formation in these models. For example, dependence of phase
coherence between W-C models of two oscillating neural masses on the ratio of
their amplitudes was demonstrated in Daffertshofer and van Wijk (2011). And in
Greenwood and Ward (2019), amplitude patterns were dissociated from phase
patterns at weak coupling strengths but both appeared together with strong
coupling. In the latter case, the simulations were run for just long enough to
reveal spatial and temporal patterns, but not long enough for amplitudes to
grow to extremely large levels. The justification for this procedure was that
oscillatory patterns in the brain are seldom stable for very long, usually
only for periods of a few hundred milliseconds. Greenwood and Ward Greenwood
and Ward (2019) suggested that there could be a mechanism that limited
oscillatory amplitudes in these cases. In particular they noted that global
inhibition or changes in local coupling strengths, among other mechanisms,
would promote the development of transient spatial and temporal amplitude and
phase patterns among neural oscillations.

In this paper we demonstrate, via numerical simulations of a Mexican-hat-
coupled model similar to that of Greenwood and Ward (2019), that sustained
phase patterns occur only when the amplitude processes are unbounded. The
introduction of ‘plastic systemic inhibition,’ however, can control the
amplitudes while still permitting sustained phase pattern formation to occur.
By ‘plastic systemic inhibition,’ we mean a system-dependent regulatory
mechanism that targets the intrinsic damping parameter, $\lambda$, of each
coupled EI-pair, instead of bounding each voltage by a sigmoid. We also tested
a ‘static’ inhibitory mechanism, which is analogous to the plastic mechanism
except that the effect on $\lambda$ is fixed instead of adaptive. As it turns
out, only with the use of the adaptive plastic systemic inhibitory mechanism
do we see bounded amplitudes combined with sustained phase patterns; the
static mechanism bounds amplitudes but does not allow for pattern formation in
the system. Additionally, we study how the eigenvalues of the total coupled
system influence amplitude growth and phase pattern formation, and show that
pattern formation is apparently stochastic when using static systemic
inhibition.

We conjecture that our plastic systemic inhibitory mechanism may be closely
related to homeostatic inhibitory plasticity Vogels  _et al._ (2011).
Inhibitory plasticity has come into recent interest as a mechanism that can
induce both stability and rich dynamics in neural networks. The role of
inhibitory plasticity in maintaining the excitation-inhibition balance, the
stabilization of recurrent network dynamics, and sensory-response de-
correlation is discussed in Sprekeler (2017). Synaptic plasticity in
inhibitory synapses is explored in Vogels  _et al._ (2011), which explains
sparse firing patterns observed in response to natural stimuli, as well as
providing a homeostatic mechanism that generates asynchronous and irregular
network states. In Hellyer  _et al._ (2016), a local homeostatic inhibitory
plasticity scheme is shown to regulate network activity and cause rich and
spontaneous dynamics to emerge over a large range of brain configurations,
which otherwise have a limited range of dynamic regimes.

In Section II.1 we summarize the derivation of the phase and amplitude
processes with a generic linear-type coupling scheme in the case where every
EI-pair has the same parameters. Section II.2 describes the matrix equations
that summarize the dynamics of the coupled system, analogous to (5). Section
II.3 describes the systemic inhibitory mechanisms. The choice of coupling
coefficients, the parameters for the EI-pairs, and the simulation parameters
for the numerical results are provided in Section II.4. The results of our
simulations are presented in Section III, with a more general discussion in
Section IV, and the conclusion in Section V.

##  II Model Development

###  II.1 Itô Transformation to Phase and Amplitude Processes

We consider a system of $N$ EI-pairs,
$(E_{1},I_{1}),...,(E_{k},I_{k}),...,(E_{N},I_{N})$, each of which, in the
absence of coupling, obeys (5), i.e.,

$d\mathbb{V}_{k}(t)=[E_{k}(t),I_{k}(t)]^{\top}=-\mathbb{A}\mathbb{V}_{k}(t)dt+\mathbb{N}d\mathbb{W}_{k}(t),$
(8)

where for each $k$, $d\mathbb{W}_{k}(t)=(dW_{k}^{E},dW_{k}^{I})^{\top}$, where
$\\{W_{k}^{E},W_{k}^{I},k=1,...,N\\}$ are independent standard Brownian
motions. Note that the matrices $\mathbb{A},\mathbb{N}$ are identical for each
EI-pair. We introduce a coupling term, $\mathbb{M}_{k}(t;\mathbb{C})$, into
(8), given by,

$\mathbb{M}_{k}(t;\mathbb{C})=\sum_{j}\mathbb{C}_{k,j}\mathbb{V}_{j}(t)$ (9)

where $\mathbb{C}=(\mathbb{C}_{k,j})$ is a coupling matrix whose diagonal
entries are zero and whose off-diagonal entries, $\\{\mathbb{C}_{k,j}\\}$,
represent the post-synaptic connection of EI-pair $j$ to EI-pair $k$. The
introduction of the coupling term modifies (8) into,

$d\mathbb{V}_{k}(t)=\big{(}-\mathbb{A}\mathbb{V}_{k}(t)+\mathbb{M}_{k}(t;\mathbb{C})\big{)}dt+\mathbb{N}d\mathbb{W}_{k}(t),$
(10)

Next we transform $-A$ into “normal form”, which is a change of basis that
aims to replace the matrix $-A$ in (10) with a matrix that is written in terms
of the intrinsic damping $\lambda$ and frequency $\omega$ of the system (5).
The intrinsic damping and frequency are the real and imaginary parts,
respectively, of the matrix $-A$; in other words $-A$ has eigenvalues given by
$-\lambda\pm\omega i$. These can be calculated as,

$\begin{split}\lambda=-\frac{tr(\mathbb{-A})}{2}=\frac{1-S_{EE}}{2\tau_{E}}+\frac{1+S_{II}}{2\tau_{I}},\end{split}$
(11)

and,

$\begin{split}&\omega=\frac{1}{2}\sqrt{tr(-\mathbb{A})^{2}-4\det(-\mathbb{A})}\\\
&=\sqrt{\lambda^{2}-\frac{S_{IE}S_{EI}+(1-S_{EE})(1+S_{II})}{\tau_{E}\tau_{I}}}.\end{split}$
(12)

We perform the change of basis
$\\{\mathbb{Y}_{k}=\mathbb{Q}_{k}^{-1}\mathbb{V}_{k}=(u_{k},v_{k})^{\top}\\}$
where the matrix $\mathbb{Q}$ is given by,

$\begin{split}\mathbb{Q}=\begin{pmatrix}-\omega&\lambda+(-1+S_{EE})\tau_{E}^{-1}\\\
0&S_{IE}\tau_{I}^{-1}\end{pmatrix},\end{split}$ (13)

to obtain the matrix $\mathbb{B}$ as,

$\begin{split}\mathbb{B}:=\mathbb{Q}^{-1}(-\mathbb{A})\mathbb{Q}=\begin{pmatrix}-\lambda&\omega\\\
-\omega&-\lambda\end{pmatrix}.\end{split}$ (14)

The coupling term is stable under the transformation, i.e.
$\mathbb{Q}^{-1}\mathbb{M}_{k}(t;\mathbb{C})\mathbb{Q}=\mathbb{M}_{k}(t;\mathbb{C})$.
Thus the system in the new basis is given by,

$d\mathbb{Y}_{k}(t)=(\mathbb{B}\mathbb{Y}_{k}(t)+\mathbb{M}_{k}(t;\mathbb{C}))dt+\mathbb{E}d\mathbb{W}_{k}(t),$
(15)

where $\mathbb{E}=\mathbb{Q}^{-1}\mathbb{N}$.

The derivation of the phase and amplitude processes corresponding to (15) are
given in Appendix A. In the case where every EI-pair has identical parameters,
these processes are given respectively by,

$dZ_{k}=\bigg{(}\frac{a^{2}+b^{2}+c^{2}-(a^{2}+b^{2})\cos(\theta_{k}(t))^{2}}{2Z_{k}}\frac{-c^{2}\sin(\theta_{k}(t))^{2}-bc\sin(2\theta_{k}(t))}{2Z_{k}}-\lambda
Z_{k}+\sum_{j}\mathbb{C}_{k,j}Z_{j}\cos(\theta_{k}-\theta_{j})\bigg{)}dt+dR_{k}.$
(16)

and,

$d\theta_{k}=\bigg{(}bc\frac{1-2\cos(\theta_{k}(t))^{2}}{Z_{k}^{2}}-\omega+\sum_{j}\mathbb{C}_{k,j}\frac{Z_{j}}{Z_{k}}\sin(\theta_{j}-\theta_{k})\bigg{)}dt+dS_{k},$
(17)

where,

$a=\frac{-\sigma_{E}}{\omega\tau_{E}}$ (18) $c=\frac{\sigma_{I}}{S_{IE}},$
(19) $b=\frac{-1+S_{EE}+\lambda\tau_{E}}{\omega\tau_{E}}c,$ (20)

and where the noise terms are given by,

$\displaystyle\begin{split}&dR_{k}=a\cos(\theta_{k}(t))dW_{k}^{E}\\\
&+(b\cos(\theta_{k}(t))+c\sin(\theta_{k}(t)))dW_{k}^{I},\end{split}$ (21)

and,

$\displaystyle\begin{split}&dS_{k}(t)=\\\
&\frac{-a\sin(\theta_{k}(t))dW_{k}^{E}+(c\cos(\theta_{k}(t))-b\sin(\theta_{k}(t)))dW_{k}^{I}}{Z_{k}}.\end{split}$
(22)

###  II.2 Matrix Form

We write the system of $N$ coupled EI-pairs as a matrix equation. We let
$\mathcal{Y}=(\mathbb{Y}_{1}^{\top},...,\mathbb{Y}_{N}^{\top})^{\top}$ and
$\mathcal{W}=(\mathbb{W}_{1}^{\top},...,\mathbb{W}_{N}^{\top})^{\top}$. The
block matrix $\mathcal{D}$, given by,

$\begin{split}\mathcal{D}=\begin{pmatrix}\mathbb{B}&0&\ldots&0&0\\\
0&\mathbb{B}&0&\ldots&0\\\ \vdots&0&\ddots&&\vdots\\\
0&\vdots&&\mathbb{B}&0\\\ 0&0&\ldots&0&\mathbb{B}\end{pmatrix},\end{split}$
(23)

captures the individual uncoupled deterministic dynamics of each EI-pair as in
(8) (without the noise term). Next, the block matrix
$\mathcal{K}(\mathbb{C})$, given by,

$\begin{split}\leavevmode\resizebox{381.5877pt}{}{$\mathcal{K}(\mathbb{C})=\begin{pmatrix}0&\mathbb{C}_{1,2}\mathbb{I}_{2}&\mathbb{C}_{1,3}\mathbb{I}_{2}&\ldots&\mathbb{C}_{1,N}\mathbb{I}_{2}\\\
\mathbb{C}_{2,1}\mathbb{I}_{2}&0&\mathbb{C}_{2,3}\mathbb{I}_{2}&\ldots&\mathbb{C}_{2,N}\mathbb{I}_{2}\\\
\mathbb{C}_{3,1}\mathbb{I}_{2}&\mathbb{C}_{3,2}\mathbb{I}_{2}&0&\ldots&\mathbb{C}_{3,N}\mathbb{I}_{2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\mathbb{C}_{N,1}\mathbb{I}_{2}&\mathbb{C}_{N,2}\mathbb{I}_{2}&\mathbb{C}_{N,3}\mathbb{I}_{2}&\ldots&0\end{pmatrix},$}\end{split}$
(24)

captures the coupling aspect of the system in (10), where $\mathbb{I}_{2}$ is
the $2\times 2$ identity matrix. The matrix $\mathcal{E}$, given by,

$\begin{split}\mathcal{E}=\begin{pmatrix}\mathbb{E}&0&\ldots&0&0\\\
0&\mathbb{E}&0&\ldots&0\\\ \vdots&0&\ddots&&\vdots\\\
0&\vdots&&\mathbb{E}&0\\\ 0&0&\ldots&0&\mathbb{E}\end{pmatrix},\end{split}$
(25)

provides the correct transformation of the i.i.d. Brownian motions in
$\mathcal{W}$. Thus the system of $N$ EI-pairs governed by (10) can be written
as the matrix equation,

$\begin{split}\leavevmode\resizebox{381.5877pt}{}{$d\mathcal{Y}(t)=\big{(}\mathcal{D}\mathcal{Y}(t)+\mathcal{K}(\mathbb{C})\mathcal{Y}(t)\big{)}dt+\mathcal{E}d\mathcal{W}(t).$}\end{split}$
(26)

Note that the matrix $\mathcal{D}$ is a constant matrix, defined by the choice
of parameters in (8), whereas the matrix $\mathcal{K}(\mathbb{C})$ depends
both on the parameters in (8), and on the choice of coupling matrix
$\mathbb{C}$. In our numerical simulations and analyses that follow, we are
concerned with the deterministic temporal eigenvalues of our coupled system,
which are the eigenvalues of the block matrix $\mathcal{L}(\mathbb{C})$, given
by,

$\displaystyle\mathcal{L}(\mathbb{C})=\mathcal{D}+\mathcal{K}(\mathbb{C})$
$\displaystyle=\begin{pmatrix}\mathbb{B}&\mathbb{C}_{1,2}\mathbb{I}_{2}&\mathbb{C}_{1,3}\mathbb{I}_{2}&\ldots&\mathbb{C}_{1,N}\mathbb{I}_{2}\\\
\mathbb{C}_{2,1}\mathbb{I}_{2}&\mathbb{B}&\mathbb{C}_{2,3}\mathbb{I}_{2}&\ldots&\mathbb{C}_{2,N}\mathbb{I}_{2}\\\
\mathbb{C}_{3,1}\mathbb{I}_{2}&\mathbb{C}_{3,2}\mathbb{I}_{2}&\mathbb{B}&\ldots&\mathbb{C}_{3,N}\mathbb{I}_{2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\mathbb{C}_{N,1}\mathbb{I}_{2}&\mathbb{C}_{N,2}\mathbb{I}_{2}&\mathbb{C}_{N,3}\mathbb{I}_{2}&\ldots&\mathbb{B}\end{pmatrix},$
(27)

###  II.3 Systemic Inhibitory Schemes

Now we introduce the systemic inhibitory mechanisms investigated in this
paper. The first is a ‘static’ systemic inhibitory mechanism where we subtract
$\delta\mathbb{I}_{2N\times 2N}$, for a systemic inhibition parameter
$\delta>0$, from (26) leading to the matrix equation for the statically
inhibited system given by,

$d\mathcal{Y}(t)=\bigg{(}\big{(}\mathcal{L}(\mathbb{C})-\delta\mathbb{I}_{2N\times
2N}\big{)}\mathcal{Y}(t)\bigg{)}dt+\mathcal{E}d\mathcal{W}(t).$ (28)

We also explore two ‘plastic’ systemic inhibitory mechanisms, where we instead
subtract $\delta\xi(t)$, for a systemic inhibition parameter $\delta>0$, from
(26) resulting in the matrix equation,

$d\mathcal{Y}(t)=\bigg{(}\big{(}\mathcal{L}(\mathbb{C})-\delta\xi(t)\big{)}\mathcal{Y}(t)\bigg{)}dt+\mathcal{E}d\mathcal{W}(t).$
(29)

The first plastic systemic inhibition mechanism is what we will call ‘binary-
type’ plastic systemic inhibition, where $\xi$ is given by,

$\displaystyle\xi(t)=\begin{pmatrix}T_{1}(t)\mathbb{I}_{2}&0&\ldots&0&0\\\
0&T_{2}(t)\mathbb{I}_{2}&0&\ldots&0\\\ \vdots&0&\ddots&&\vdots\\\
0&\vdots&&T_{N-1}(t)\mathbb{I}_{2}&0\\\
0&0&\ldots&0&T_{N}(t)\mathbb{I}_{2}\end{pmatrix},$ (30)

where $\mathbb{I}_{2}$ is the $2\times 2$ identity matrix and where,

$\displaystyle
T_{i}(t)=1\\{Z_{i}(t)>z^{*}\\}=\begin{cases}1&\text{if}\,\,Z_{i}(t)>z^{*}\\\
0&\text{otherwise}\\\ \end{cases},$

for a plastic amplitude threshold $z^{*}$. The second plastic mechanism we
explored is what we will call ‘saturation-type’ plastic systemic inhibition,
where $\xi(t)$ is given by,

$\displaystyle\xi(t)=\begin{pmatrix}U_{1}(t)\mathbb{I}_{2}&0&\ldots&0&0\\\
0&U_{2}(t)\mathbb{I}_{2}&0&\ldots&0\\\ \vdots&0&\ddots&&\vdots\\\
0&\vdots&&U_{N-1}(t)\mathbb{I}_{2}&0\\\
0&0&\ldots&0&U_{N}(t)\mathbb{I}_{2}\end{pmatrix},$ (31)

where

$\displaystyle U_{i}(t)=\frac{1}{1+\max\\{0,z^{*}-Z_{i}(t)\\}}.$

As we will see in Section III, both plastic inhibitory mechanisms have the
desired effect of keeping amplitudes stochastically bounded while allowing
pattern formation to occur in the phases. The saturation-type inhibition
mechanism is more biologically plausible, however, as the amplitude threshold
$z^{*}$ in that mechanism corresponds to a maximum, or ‘saturated,’ rate of
firing of a neuron or of a neural population.

When choosing a systemic inhibition parameter $\delta$ so that the resulting
system has a specific maximum eigenvalue real-part when systemic inhibition is
applied, it is practical to instead first pick a target maximal eigenvalue
real-part for the matrix $\mathcal{L}(\mathbb{C})-\delta\mathbb{I}_{2N\times
2N}$ in (28), or $\mathcal{L}(\mathbb{C})-\delta\xi(t)$ in (29), and then
choose $\delta$ accordingly. We did this because the eigenvalue real parts
provide a heuristic as to whether the system remains close to or deviates far
from the origin ($Z=0$). For example, one would expect that if the eigenvalues
of the matrix $\mathcal{L}(\mathbb{C})-\delta\mathbb{I}_{2N\times 2N}$ all
have negative real part, then the system will remain bounded – which is
exactly what we are trying to achieve by raising the systemic inhibition when
the amplitude is too large. Conversely, we expect amplitudes to be unbounded
when there is at least one eigenvalue with a positive real part. In the case
where $\mathcal{L}(\mathbb{C})$ is diagonalizable, subtracting
$\delta\mathbb{I}_{2N\times 2N}$ will subtract $\delta$ from all of the
eigenvalues. Hence when $\mathcal{L}(\mathbb{C})$ is diagonalizable, we can
numerically find the maximal eigenvalue real part, say $\lambda^{*}$, and pick
a target maximal eigenvalue real part, say $\Lambda^{*}$, and set
$\delta=\lambda^{*}-\Lambda^{*}$. Then
$\mathcal{L}(\mathbb{C})-\delta\mathbb{I}_{2N\times 2N}$ will have maximal
eigenvalue real-part $\Lambda^{*}$. It is not as clear whether this method
would work in general when the matrix $\mathcal{L}(\mathbb{C})$ is not
diagonalizable, but in our simulations that follow we re-compute the maximal
eigenvalue real part to ensure this works as intended. We will hereafter refer
to $\Lambda^{*}$ as the ‘adjusted eigenvalue bound.’ See Section IV for a
discussion on the biological interpretation of these systemic inhibitory
mechanisms.

###  II.4 Parameter Specifications and Mexican Hat Coupling

We consider our EI-pairs regarded as a one dimensional ring-lattice with
uniform spacing. That is, EI-pair $k$ has the integer position $k$ on the
ring-lattice $\\{1,...,N\\}\mod N$. Every EI-pair is set to have parameters
$\tau_{E}=0.003,\tau_{I}=0.006,S_{EE}=1.5,S_{II}=0.1,S_{IE}=4,$ and
$S_{EI}=1$, whereas the damping rates (real parts of eigenvalues of
$\mathbb{A}$) and intrinsic frequencies (imaginary parts of eigenvalues of
$\mathbb{A}$) were adjusted to fit the other parameters according to (11) and
(LABEL:omega). Population size is $N=100$. Although the specifications of our
simulations closely follow that of Greenwood and Ward (2019), that work did
not specify the constants $\sigma_{E}$ and $\sigma_{I}$, and so we used
$\sigma_{E}=\sigma_{I}=12$ as in Greenwood  _et al._ (2016). As mentioned in
Greenwood  _et al._ (2016), these sets of parameters yield a narrow-band
quasi-cycle oscillation at approximately 70 Hz, which belongs to the gamma
frequency range of 30-80 Hz which has been shown to be important in neural
oscillations Greenwood  _et al._ (2015). The Euler-Maruyama method was used to
numerically solve the stochastic differential equations in (16) and (17) with
time step $\Delta t=5\cdot 10^{-5}$

In this paper we use the ‘Mexican hat’ function, which is a difference of
Gaussian density functions, as a kernel in our coupling scheme. It is given
by,

$\begin{split}&m(x)=b_{1}\exp\big{[}-(\frac{x}{d_{1}})\big{]}-b_{2}\exp\big{[}-(\frac{x}{d_{2}})\big{]},\\\
&b_{2}>b_{1},\hskip 2.84544ptd_{2}>d_{1},\end{split}$ (32)

where $b_{1}$ and $b_{2}$ are the amplitudes and $d_{1}$ and $d_{2}$ the
diffusive parameters for two Gaussians. That is, we define the coupling matrix
$\mathbb{C}$ by,

$\begin{split}\mathbb{C}_{k,j}=Cm(k-j),\end{split}$ (33)

where $C>0$ is an adjustable parameter. Note that the difference $k-j$ in
$m(k-j)$ is actually the minimum distance between points $k$ and $j$
respecting the ring structure of the lattice. The coupling term,
$\mathbb{M}_{k}(t;\mathbb{C})$, in (10) is thus given by,

$\begin{split}\mathbb{M}_{k}(t;\mathbb{C})=\sum_{j=1}^{N}Cm(k-j)\mathbb{V}_{j}(t).\end{split}$
(34)

The spacing of the EI-pairs can be defined implicitly by a scaling of the
parameters $d_{1}$ and $d_{2}$ in (LABEL:mexhat), and for this reason we let
the EI-pairs be spaced one unit apart in (33). In Greenwood and Ward (2019),
$b_{2}$ and $d_{1}$ were fixed to be 1 (in order to make the analysis more
tractable), and then $b_{1}=1.3$ and $d_{2}=1.5$ were determined to be
suitable parameters for achieving pattern formation. Unlike in Greenwood and
Ward (2019), we employ unit distances between EI-pairs instead of a separation
distance of 0.2 per EI-pair, and consider diffusion parameters roughly five
times greater in order to account for this difference. In Greenwood and Ward
(2019), patterns were found to form with $C=1$ and more so at $C=8$, and so
here we chose $C=8$.

In order to measure the degree of spatial pattern formation in the simulated
phases, we plot the spatial sample entropy as a function of time. The version
of sample entropy that we use in this paper is defined in Richman and Moorman
(2000) – which provides background motivation as well as technical details.
This sample entropy is an unbiased estimator of the ‘Rényi entropy of order
2.’ Details are given in the ‘Sample Entropy’ section in the Supplementary
Material. In our context, high values of entropy suggest absence of a pattern,
whereas low values of entropy suggest the presence of a pattern. We computed
sample entropies using an algorithm given in Ramdani  _et al._ (2009) via an
implementation provided by Martínez-Cagigal (2018) (wherein we set tolerance
parameter $r=1$, embedding dimension parameter $m=1$, and use Chebychev
distance measure). Sample entropy provides a useful heuristic to aid visual
inspection and to quantify the degree of pattern formation.

For the simulations that follow in this paper we consider Mexican hat coupling
kernels for 7 sets of parameters $(b_{1},b_{2},d_{1},d_{2})$, but with the
values of $b_{2}=1,d_{1}=5$ held constant throughout. We enumerate the
subsequent couplers (i.e. the matrix $\mathbb{C}$ for each system) as couplers
A through G, as displayed in Table 1. The parameters control the shape of the
Mexican Hat function (see Figs. 1 and 2) and thus the visual appearance of the
phase patterns; larger values of $d_{2}$ relative to $d_{1}$ imply a wider
reach of the Mexican hat, and different ratios $b_{1}/b_{2}$ determine the
relative amount of inhibition around the excitatory centre. Our focus here is
on the formation of phase patterns rather than on their visual appearance, so
we won’t comment further on the character of the phase patterns although their
differences will be apparent and can be related to the parameters of the
respective Mexican hats.

The purpose of choosing a diverse range of couplers was to investigate whether
results regarding the presence or absence of pattern formation vary for
different values of coupler parameters that had a variety of max eigenvalue
real-parts. Some couplers in Table 1 are more realistic than others; e.g.
coupler A has a balance of excitation and inhibition, whereas coupler C has
almost no inhibition. Nevertheless, we were interested in testing for pattern
formation across a variety of dynamic regimes, and so we also simulated the
more unrealistic couplers. As it turns out the differences in the absence or
presence of pattern formation were limited. Couplers A through D behaved
essentially the same with regard to pattern formation, but differently than
couplers E through G – the three of which behaved essentially the same. In
Section 3 we present only the simulations for couplers A and F, relegating the
others to the Supplementary Material.

Pair | $b_{1}$ | $d_{2}$ | Max Eigenvalue Real-part
---|---|---|---
A | 2.6 | 19.1 | 115.8
B | 4.1 | 9.1 | 137.8
C | 3.6 | 3.6 | 175.1
D | 4.1 | 19.1 | 189.0
E | 1.1 | 5.1 | -3.462
F | 1.1 | 6.1175 | -0.0004060
G | 1.1 | 6.12 | 0.01516

Table 1: Parameter values in the coupling kernel (LABEL:mexhat) for each of
the 7 couplers. Figure 1: Mexican hat coupling kernels, as in (LABEL:mexhat),
with couplers: (A) $(b_{1},d_{2})=(2.6,19.1)$ with max eigenvalue real-part
115.8 (dotted curve), (B) $(b_{1},d_{2})=(4.1,9.1)$ with max eigenvalue real-
part 137.8 (solid curve), (C) $(b_{1},d_{2})=(3.6,3.6)$ with max eigenvalue
real-part 175.1 (dotted and dashed curve), (D) $(b_{1},d_{2})=(4.1,19.1)$ with
max eigenvalue real-part 189 (dashed curve). The parameters $b_{2}$ and
$d_{1}$ are set to 1 and 5, respectively, for each coupler. The line $y=0$ is
dashed for reference. It is relevant to note that the more realistic couplers
in this figure are coupler A (dotted curve) and D (dashed curve) as both have
a balance of excitation and inhibition. By contrast, coupler B (solid curve)
has very little inhibition compared to its excitation and coupler C (dotted
and dashed curve) has almost no inhibition. Figure 2: Mexican hat coupling
kernels, as in (LABEL:mexhat), with couplers: (E) $(b_{1},d_{2})=(1.1,5.1)$
with max eigenvalue real-part -3.461 (dotted curve), (F)
$(b_{1},d_{2})=(1.1,6.1175)$ with max eigenvalue real-part -0.04708 (solid
curve), (G) $(b_{1},d_{2})=(1.1,6.12)$ with max eigenvalue real-part 0.01516
(dashed curve). The parameters $b_{2}$ and $d_{1}$ are set to 1 and 5,
respectively, for each coupler. The line $y=0$ is dashed for reference. It is
relevant to note that coupler E (dotted curve) is unrealistic as it has
essentially no inhibition, whereas couplers F (solid curve) and G (dashed
curve) both have some balance of excitation and inhibition.

##  III Creating Sustained Phase Patterns with Bounded Amplitudes

We simulated (16) and (17) with the Mexican hat coupling parameters indicated
in Table 1 for a range of time periods depending on the parameters and the
results of the simulations. Time periods were generally shorter when
amplitudes grew without bound, and longer when amplitudes remained bounded, in
order to demonstrate the outcome convincingly. In what follows we display only
a subset of the simulation plots that illustrate our overall results. Results
not displayed but described can be found in the Supplementary Material.

###  III.1 Investigating Phase Patterns without Systemic Inhibition

Couplers A through D all induce positive maximum eigenvalue real-parts with an
order of magnitude of $10^{2}$ for their respective matrices
$\mathcal{L}(\mathbb{C})$ (II.2). The initial conditions for simulations
without systemic inhibition using these pairs have amplitudes taken from a
uniform distribution on $(0,1]$, and phases taken from a uniform distribution
on $[0,2\pi]$. We ran these simulations in order to illustrate clearly how the
model behaves when maximum eigenvalue real-parts are large and positive. Under
these conditions amplitudes grow very quickly and apparently without bound. A
representative example of these results is displayed in Fig. 3 for coupler A
(similarly for couplers B through D; see Figs. S1-S3). A clear pattern
develops over a short time in the phases as well as in the amplitudes. Sample
entropy of the phases is high initially but quickly decreases to a low level
as the phase pattern develops and evolves. In all of these four simulations,
the amplitudes grow to extreme orders of magnitude in a short period of time.
As expected, the systems with larger maximal eigenvalue real-parts tend to
have amplitudes with greater orders of magnitude when compared on equal time
scales; e.g. the system with coupler A has maximal eigenvalue real part 115.8
and amplitudes that reach a magnitude of $10^{4}$ at $t=0.06$ (Fig. 3) whereas
the system with coupler D has maximal eigenvalue real part of 189.0 and
amplitudes that reach magnitudes of $10^{5}$ at $t=0.06$ (Fig. S3).

Figure 3: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters was 115.8. Initial amplitudes were selected from
a uniform distribution on $(0,1]$ and initial phases from a uniform
distribution on $[0,2\pi]$.

We also simulated (16) and (17) for couplers E through G without systemic
inhibition. Fig. 4 displays the results using coupler F (similar results for E
and G in Figs. S4-S5, respectively), which has a very small negative real part
of the eigenvalues, and thus some damping. In this case the amplitudes do not
grow without bound but remain in the vicinity of $2\cdot 10^{2}$, even though
we ran the simulation for 1 sec. Unlike the results for large positive real
eigenvalues, no clear pattern develops in the phases or in the amplitudes, and
the sample entropy remains high. Even when the real part of the eigenvalues is
positive, but very small, this result is the same (Fig. S5).

Figure 4: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler F; $(b_{1},b_{2},d_{1},d_{2})=(1.1,1,5,6.1175)$.
The maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$),
given by (II.2), with these parameters was -0.00040595. Initial amplitudes
were selected from a uniform distribution on $(0,1]$ and initial phases from a
uniform distribution on $[0,2\pi]$.

Additional simulations were run for couplers E through G with the same
parameters except that the initial amplitudes were sampled from a uniform
distribution on (0,1]+5000, and the simulations were run for a longer time
period: 2-3 sec. Figure 5 displays the results for coupler F (similar results
for couplers E and G in Figs. S6-S7, respectively). When the amplitudes are
large, at the beginning of the run, clear spatial patterns similar to those
observed in Fig. 3 (and Figs. S1-S3) develop in the phases. These patterns
decay over time, however, as the amplitudes damp to near zero. The sample
entropies of the phases decrease initially as the spatial pattern appears but
then increases again as the pattern decays.

Figure 5: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler F; $(b_{1},b_{2},d_{1},d_{2})=(1.1,1,5,6.1175)$.
The sample entropy of the spatial phase distribution is also plotted for every
time point with tolerance $\tau=1$ and embedding dimension $m=1$. The maximum
eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by (II.2),
with these parameters is -0.00040595. Initial amplitudes are selected from a
uniform distribution on $(0,1]+5000$ and initial phases from a uniform
distribution on $[0,2\pi]$.

The results of the simulations illustrated in Fig. 5 (and in Figs. S6-S7)
imply that the lack of pattern formation observed in Fig. 4 (and in Figs.
S4-S5) is not because the simulations were run for an insufficient time period
for large amplitudes to develop. We have seen that phase pattern formation
occurs when the amplitudes appear to increase without bound. Therefore, if
patterns were to form at some point in time, we would also expect large
amplitudes to develop at or before that time point. But when we set the
amplitudes to large initial values, but with small positive or negative real
parts of the eigenvalues, the amplitudes rapidly decrease from those large
initial values, and the phase patterns that are initially present disappear.
It thus seems that the appearance of phase patterns in this model is dependent
on the sustaining of large amplitudes.

###  III.2 Static Systemic Inhibition Prevents and Destroys Phase Patterns

To study the effect of static systemic inhibition to bound amplitudes, we
simulated the model for couplers A through D, as specified in Table 1, with
initial amplitudes taken from a uniform distribution on $(0,1]$, and phases
taken from a uniform distribution on $[0,2\pi]$. The adjusted eigenvalue bound
$\Lambda^{*}=-10^{-3}$ (defined in Section II.4) was used for each simulation.
This bound was chosen so that the maximal eigenvalue real parts would all be
negative – ensuring that the amplitudes remain bounded – but close to zero to
minimize the magnitude of the inhibitory effect (i.e. to minimize $\delta$
which is defined in Section II.4). Fig. 6 displays a representative example of
these results for coupler A (similar results for couplers B through D; see
Figs. S8-S10). The results are qualitatively similar to Fig. 5 (and to Figs.
S6-S7), however the phase patterns seem to disappear faster as the amplitude
damps for these simulations with static systemic inhibition.

Figure 6: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters is 115.8. Initial amplitudes are selected from a
uniform distribution on $(0,1]+5000$ and initial phases from a uniform
distribution on $[0,2\pi]$. Static systemic inhibition (28) was used with
$\delta=115.8505$ corresponding to an adjusted eigenvalue bound of
$\Lambda^{*}=-10^{-3}$.

We also tried using a positive adjusted eigenvalue bound (still less than the
maximum eigenvalue real-part). This bound still leads to the application of
static systemic inhibition, but the resulting maximal eigenvalue real-part of
the inhibited system in (28) is positive, and not negative as with a negative
adjusted eigenvalue bound. The problem is to determine whether there is an
adjusted eigenvalue bound $\Lambda^{*}$ (or equivalently a systemic inhibition
parameter $\delta$) for which we will observe sustained phase patterns. As it
turns out, for certain choices of $\Lambda^{*}$, i.i.d. simulations do
sometimes, but not always, display sustained phase patterns. This is
demonstrated in Figs. 7 and 8. Both simulations used coupler A, independently
sampled initial phases from a uniform distribution on $[0,2\pi]$,
independently sampled initial amplitudes from a uniform distribution on
$(0,1]+1000$, used static systemic inhibition with an adjusted eigenvalue
bound of $\Lambda^{*}=10$, and were run independently. We see phase pattern
formation, and decreased sample entropy, in Fig 7 but not in Fig. 8. Note that
the amplitudes in Fig. 7 grow large and seemingly without bound similar to
previous simulations with pattern formation, whereas the amplitudes in Fig. 8
damp to low values similar to previous simulations without pattern formation.

Figure 7: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters is 115.8. Initial amplitudes are selected from a
uniform distribution on $(0,1]+1000$ and initial phases from a uniform
distribution on $[0,2\pi]$. The systemic inhibition parameter $\delta=105.9$
corresponding to an adjusted eigenvalue bound of $\Lambda^{*}=10$.

Figure 8: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters is 115.8. Initial amplitudes are selected from a
uniform distribution on $(0,1]+1000$ and initial phases from a uniform
distribution on $[0,2\pi]$. Static systemic inhibition (28) was used with
$\delta=105.9$ corresponding to an adjusted eigenvalue bound of
$\Lambda^{*}=10$.

The results in Figs. 7 and 8 suggest that both amplitude increase and
sustained phase pattern formation occur stochastically, i.e. that pattern
formation only occurs with some probability under static systemic inhibition
with a small positive real part of the eigenvalues. To explore this further,
we performed nine i.i.d. simulations using the same specifications as Figs. 7
and 8 for a duration of $T=2$, for a range of adjusted eigenvalue bounds
$\Lambda^{*}=8,9,10,11$, and $12$. Box plots of the amplitude distributions at
time $t=2$ are given for each simulation in Fig. 9. For each of the adjusted
eigenvalue bounds, the amplitude processes grow to be very large only
sometimes. Moreover, it is worthwhile to note that only one of nine
simulations have amplitude growth for adjusted eigenvalue bounds of
$\Lambda^{*}=8$ or $\Lambda^{*}=9$ (when inhibition is set to $\delta=107.8$
and $\delta=106.8$, respectively). By comparison, this occurred for three of
nine simulations with adjusted eigenvalue bounds of 10 or 11 ($\delta=105.8$
and $\delta=104.8$, respectively), and for four of 10 simulations with
adjusted eigenvalue bounds of 12 ($\delta=103.8$). This suggests that the
probability of amplitudes growing large (which has been demonstrated to be
associated with pattern formation in the phases) decreases as static systemic
inhibition increases.

Figure 9: Independent samples of terminal distributions of simulated amplitude
processes run for $T=2$ with Mexican Hat coupling (LABEL:mexhat) with coupler
A; $(b_{1},b_{2},d_{1},d_{2})=(1.1,1,5,6.12)$, for various adjusted eigenvalue
bounds. The maximum eigenvalue real part of the matrix
$\mathcal{L}(\mathbb{C}$), given by (II.2), with these parameters is 0.01516.
Initial amplitudes are selected from a uniform distribution on $(0,1]+1000$
and initial phases from a uniform distribution on $[0,2\pi]$. Note that each
of the five simulations corresponding to the five box plots in each subplot
are Independent as are the simulations between each subplot. Red pluses
represent outliers from the boxes.

###  III.3 Plastic Systemic Inhibition Permits Phase Patterns with Bounded
Amplitudes

To study the effects of plastic systemic inhibition, we ran simulations with
coupler A and binary-type and saturation-type plastic systemic inhibition (see
Section II.3). The results are shown in Figs. 10 and 11 for $z^{*}=100$ and in
Figs. 12 and 13 for $z^{*}=300$, for binary-type and saturation-type
inhibition, respectively. The result for binary-type plastic inhibition is
shown for $z^{*}=200$ in Fig. S11 – in which the results are intermediate
between Figs. 10 and 12. There is only a weak phase pattern apparent in Figs
10 and 11 ($z^{*}=100$), with relatively high sample entropy, but Figs. 12 and
13 ($z^{*}=300$) display clear phase patterns, with relatively small sample
entropies. There is also a clear difference between the amplitudes in these
simulations: there are no apparent amplitude patterns, but average amplitudes
range from about 100 when $z^{*}=100$, to 300 when $z^{*}=300$, as would be
expected given how the plastic systemic inhibitory mechanisms were designed.
Thus the appearance of phase patterns under these conditions depends on the
magnitude of the amplitudes, and does not depend on whether the amplitudes
display spatial patterns. Moreover, under these conditions phase patterns can
appear in the absence of very large oscillatory amplitudes, making this
situation much more similar to actual neural activity in brains. Thus plastic
systemic inhibition of the sort implemented here can limit amplitudes to more
biological ranges while allowing spatial patterns to form in the phases.

Figure 10: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters was 115.8. Initial amplitudes were selected from
a uniform distribution on $(0,1]$ and initial phases from a uniform
distribution on $[0,2\pi]$. Binary-type plastic systemic inhibition (29) was
used with $\delta=115.9$ – corresponding to an adjusted eigenvalue bound of
$\Lambda^{*}=-10^{-3}$ – with a plastic amplitude threshold of $z^{*}=100$.

Figure 11: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters was 115.8. Initial amplitudes were selected from
a uniform distribution on $(0,1]$ and initial phases from a uniform
distribution on $[0,2\pi]$. Saturation-type plastic systemic inhibition (28)
was used with $\delta=115.9$ – corresponding to an adjusted eigenvalue bound
of $\Lambda^{*}=-10^{-3}$ – with a plastic amplitude threshold of $z^{*}=100$.

Figure 12: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters was 115.8. Initial amplitudes were selected from
a uniform distribution on $(0,1]$ and initial phases from a uniform
distribution on $[0,2\pi]$. Binary-type plastic systemic inhibition (29) was
used with $\delta=115.9$ – corresponding to an adjusted eigenvalue bound of
$\Lambda^{*}=-10^{-3}$ – with a plastic amplitude threshold of $z^{*}=300$.

Figure 13: Phase and amplitude process simulation with Mexican Hat coupling
(LABEL:mexhat) with coupler A; $(b_{1},b_{2},d_{1},d_{2})=(2.6,1,5,19.1)$. The
maximum eigenvalue real part of the matrix $\mathcal{L}(\mathbb{C}$), given by
(II.2), with these parameters was 115.8. Initial amplitudes were selected from
a uniform distribution on $(0,1]$ and initial phases from a uniform
distribution on $[0,2\pi]$. Saturation-type plastic systemic inhibition (28)
was used with $\delta=115.9$ – corresponding to an adjusted eigenvalue bound
of $\Lambda^{*}=-10^{-3}$ – with a plastic amplitude threshold of $z^{*}=300$.

##  IV Discussion

###  IV.1 Comparing methods of systemic control

Beginning with the paper of Wilson and Cowan Wilson and Cowan (1972) most
studies involving neural population models have used a sigmoid function to
bound the values of the variables, which should remain in a dynamic range in
order to be biologically meaningful. Such a device has also been extended to
neural fields with coupling Faugeras and Inglis (2015). However, we were
interested in devising a mechanism that achieves the same goal without
introducing non-linearity in this way. By finding such a mechanism of bounding
variable values, we might preserve the connections that our system has to the
literature surrounding quasi-cyles, such as that in Baxendale and Greenwood
(2011), Greenwood  _et al._ (2015), Greenwood  _et al._ (2016), and Greenwood
and Ward (2019). By doing so, it is possible that important questions about
neural dynamics could be answered in future research on quasi-cycles.

Moreover, the saturation-type mechanism we proposed in this paper has a closer
connection to the underlying biological system than do the sigmoids in
previous works. Both mechanisms operate under the principle that the firing
rates of neuron populations reach a point of saturation, however our mechanism
can be interpreted as adaptively modifying the system parameters in (5) in
response to the overall magnitude of the variables. Sigmoids, on the other
hand, simply ‘correct’ the variables, by mapping them to a desired range of
values. Bounding variables by a modification of the parameters in (5) is not
new; in works such as Abeysuriya  _et al._ (2018), systemic inhibition was
implemented by modifying the parameter $S_{EI}$. However, this approach only
affects the oscillations of each EI-pair, by (LABEL:omega), and not the
damping $-\lambda$ of the EI-pairs, as per (11). By contrast, our systemic
inhibitory mechanisms subtract $\delta$ (in the static mechanism) or the
diagonals of $\delta\xi(t)$ (in the plastic mechanism) from the diagonals of
the matrix $\mathcal{L}(\mathbb{C})$, which are copies of the diagonals of
$\mathbb{B}$ containing the intrinsic damping parameter $-\lambda$. Thus, our
systemic inhibitory mechanisms are akin to selectively increasing the damping
parameter $\lambda$ for each EI-pair individually, while holding all other
parameters constant. In particular the saturation-type mechanism subtracts
$\delta(1/(1+\max(0,z^{*}-Z_{k}))$ from the damping term of each process, thus
increasing the damping parameter more the nearer the amplitude is to the
saturation threshold, $z^{*}$, which we conceptualize as the maximum firing
rate of a neural population.

For each individual, non-normalized system we have
$\lambda=((1-S_{EE})/\tau_{E}+(1+S_{II})/\tau_{I})/2$ ((11)). Hence, in terms
of the system parameters $(S_{EE},S_{IE},S_{EI},S_{II})$, changes in $\lambda$
would correspond to changes in self-excitation ($S_{EE}$) and/or in self-
inhibition ($S_{II}$). Either a decrease in self-excitation (decrease
$S_{EE}$) or an increase in self-inhibition (increase $S_{II}$) would result
in an increase in damping by making $-\lambda$ more negative. In our
normalized systems it was possible to change $\lambda$ without changing
$\omega$ by simply changing the $\lambda$ in the matrix $\mathbb{B}$ and
leaving $\omega$ unchanged (via either (28) or (29)). But, in the non-
normalized original system, (8), in order to maintain the oscillation
frequency at the same value when changing $S_{EE}$ or $S_{II}$, a change in
one or both of the interaction efficacies ($S_{EI},S_{IE}$) in (LABEL:omega)
would be necessary. It would make sense for changes in synapses of the same
neuron population to occur at the same time. For example, a decrease in
$S_{EE}$ could be accompanied by a corresponding decrease in $S_{EI}$,
increasing damping and maintaining the same oscillation frequency. A similar
mechanism could be implemented by the inhibitory population, with increases in
both $S_{II}$ and $S_{IE}$. These changes in synaptic efficacy would occur
rapidly, in an activity-dependent way, similar to the saturation-type plastic
inhibition we have studied. Such rapid changes in synaptic efficacies are
observed in real neural systems Kasai  _et al._ (2010), so this does
constitute a biologically plausible mechanism for bounding oscillatory
amplitudes.

###  IV.2 Parsing the Sources of Stability

It is natural to think of our system (10) as a coupling of the otherwise
i.i.d. Ornstein-Uhlenbeck (OU) processes satisfying (8). A further aim, after
our present results is to find behaviors that generalize over large families
of couplings, which are highly structured but at the same time may be regarded
as ‘partially random.’ Some randomness might allow us to take into account
that the brain structures of individuals develop according to a combination of
tight rules together with some stochasticity.

In order to study separately the nature of neural systems with randomness of
connectivity structure, some authors including Ipsen and Peterson (2020)
Wainrib and Touboul (2013) Rajan and Abbott (2006) have studied models of the
form

$\begin{split}\frac{dx_{i}}{dt}=-\frac{x_{i}}{\tau}+\sum_{j=1}^{n}J_{i,j}S(x_{j}),\hskip
7.11317pti=1,...,n.\end{split}$ (35)

Here $S$ is a sigmoid confining the values of $x_{i}$ to $[-1,1]$, $J$ is a
(one-way) connectivity matrix, $\frac{1}{\tau}$ is a decay rate that we can
compare to the $\lambda$ ((11)) in the present paper. The aim of Ipsen and
Peterson (2020) is to construct $J$, which is stochastic with strong
structural restrictions in terms of, e.g., excitatory-inhibitory balance, in
various forms, such that the stability of the system in terms of the
distribution of eigenvalues of $J$ is controlled with increasing system size,
$n$. They find that a critical value of $\frac{1}{\tau}$, at which the
stability of the system changes as the parameter $\tau$ increases, is
($\sigma_{\text{eff}}\sqrt{n}$). Here $\sigma_{{\text{eff}}}$ is a variance
related to $J$.

In comparing (35) with (10) we are representing our OU-type system, (10), by
just the small, negative part of its eigenvalue pair. The question of pattern
formation is absent, but the question of boundedness of the process is present
and handled by the sigmoid, $S$, in (35). The question arises: suppose we
replace $S$ with a “plastic inhibition”, increasing $\frac{1}{\tau}$ when ever
$|x_{i}|$ exceeds a certain value. Do results similar to those of Ipsen and
Peterson (2020) hold?

##  V Conclusions

In Greenwood and Ward (2019), where we studied synchronization patterns of
neural fields of coupled EI-systems in terms of their amplitude and phase
processes, we limited consideration to time intervals in which the system
remained inside some amplitude bound. There we found that quasi-cycle phases
quickly synchronized to form patterns, even with weak coupling, whereas
amplitudes formed patterns only at somewhat greater coupling strength. In the
present paper we introduced a more realistic bounding device than the sigmoid,
which we suggest can be construed as neural plasticity, and which can produce
rapid phase pattern formation, while at the same time amplitudes form a
bounded, apparently stationary, stochastic field with no apparent patterning.
This would seem to be a significant improvement in this direction over
previous approaches.

We have demonstrated, via numerical simulations, the important role of
amplitude magnitudes in the formation and sustaining of phase patterns for the
system given by (10). The only way phase patterns were stable in simulations
of (10) was if the amplitudes grew without bound. A static systemic inhibitory
mechanism ((28)) was used to attempt to bound the amplitudes while sustaining
pattern formations. Despite successfully bounding the amplitudes, such an
intervention was found to prevent and destroy phase patterns. However, a
plastic systemic inhibitory mechanism ((29)) was found to fully resolve the
issue – bounding amplitudes and allowing for sustained phase pattern
formation.

Nonetheless, there are still issues to resolve. One issue is that a mechanism
must be found in actual neural settings which fulfills the role of a plastic
systemic inhibitory mechanism as envisioned in the present paper. The
saturation-type mechanism introduced in this paper simply serves as a proof of
concept that a plausible biological function in the brain could act as such an
inhibitory mechanism, while the binary-type mechanism demonstrates that even a
crude mechanism can work. Furthermore, the conjectured connections to Vogels
_et al._ (2011), Sprekeler (2017), and Hellyer  _et al._ (2016) need to be
confirmed. Even if a plastic systemic inhibitory mechanism is discovered in
the brain that bounds amplitudes and allows phase patterns, it could possibly
operate independently of the functions described in Vogels  _et al._ (2011),
Sprekeler (2017), and Hellyer  _et al._ (2016). Another issue to resolve is to
explore the mathematics behind the stability in our system with plastic
systemic inhibition, as discussed in Section IV.2. Indeed, more work is
required to understand exactly why and how plastic systemic inhibition was the
key to bounding amplitudes while allowing phase patterns, and whether such a
mechanism is present in actual biological contexts. We hope the present paper
will stimulate such work and further interest in stochastic neural field
models.

*

##  Appendix A Change of Variables with Itô’s formula

Before carrying out the change of variables for the amplitude and phase
processes, we compute the noise terms,

$\mathbb{E}=\mathbb{Q}^{-1}\mathbb{N}=\begin{pmatrix}\frac{-\sigma_{E}}{\omega\tau_{E}}&\frac{-1+S_{EE}+\lambda_{k}\tau_{E}}{\omega
S_{IE}\tau_{E}\sigma_{I}}\\\ {0}&\frac{\sigma_{I}}{S_{IE}},\end{pmatrix}$ (36)

so that,

$\mathbb{E}d\mathbb{W}_{k}(t)=\begin{pmatrix}adW_{k}^{E}+bdW_{k}^{I}\\\
cdW_{k}^{I}\end{pmatrix},$ (37)

where,

$a=\frac{-\sigma_{E}}{\omega\tau_{E}}$ (38)
$b=\frac{-1+S_{EE}+\lambda\tau_{E}}{\omega S_{IE}\tau_{E}}\sigma_{I},$ (39)

and,

$c=\frac{\sigma_{I}}{S_{IE}}.$ (40)

In the calculations that follow, we will see the terms $(du_{k})^{2}$,
$(dv_{k})^{2}$, and $du_{k}\cdot dv_{k}$, which we will first compute here. In
computing each of $(du_{k})^{2}$ and $(dv_{k})^{2}$, there are three terms:
one term multiplied by $(dt)^{2}=0$, a second term multiplied by
$dt\cdot(adW_{k}^{E}+bdW_{k}^{I})=0$ (resp. $dt\cdot cdW_{k}^{I}=0$), and a
third term multiplied by $(adW_{k}^{E}+bdW_{k}^{I})^{2}$ (resp.
$(cdW_{k}^{I})^{2}$). These third terms simplify, respectively as,

$\begin{split}&(adW_{k}^{E}+bdW_{k}^{I})^{2}\\\
&=a^{2}(dW_{k}^{E})^{2}+b^{2}(dW_{k}^{I})^{2}+2abdW_{k}^{E}dW_{k}^{I}\\\
&=(a^{2}+b^{2})dt,\end{split}$ (41)

and,

$(cdW_{k}^{I})^{2}=c^{2}dt.$ (42)

Thus

$(du_{k})^{2}=(a^{2}+b^{2})dt$ (43)

and,

$(dv_{k})^{2}=c^{2}dt.$ (44)

When computing $du_{k}\cdot dv_{k}$, there are four terms: one multiplied by
$(dt)^{2}=0$, another multiplied by $dt\cdot adW_{k}^{E}+bdW_{k}^{I}=0$,
another multiplied by $dt\cdot cdW_{k}^{I}=0$, and the fourth term is
$(adW_{k}^{E}+bdW_{k}^{I})\cdot(cdW_{k}^{I})$, which we expand as,

$\begin{split}&(adW_{k}^{E}+bdW_{k}^{I})\cdot(cdW_{k}^{I})\\\
&=acdW_{k}^{E}\cdot dW_{k}^{I}+bcdW_{k}^{I}\cdot dW_{k}^{I}\\\
&=bcdt,\end{split}$ (45)

so that,

$du_{k}\cdot dv_{k}=bcdt.$ (46)

Now to change variables we use Itô’s formula, which says that for a smooth
function $f:\mathbb{R}^{2}\rightarrow\mathbb{R}$, we have,

$df\begin{pmatrix}u_{k}\\\ v_{k}\end{pmatrix}=(\nabla
f)^{\top}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}+\frac{1}{2}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}^{\top}Hf\begin{pmatrix}du_{k}\\\ dv_{k}\end{pmatrix},$
(47)

where $Hf$ is the Hessian matrix,

$Hf\begin{pmatrix}u_{k}\\\
v_{k}\end{pmatrix}=\begin{pmatrix}\frac{\partial^{2}f}{\partial
u_{k}^{2}}&\frac{\partial^{2}f}{\partial u_{k}\partial v_{k}}\\\
\frac{\partial^{2}f}{\partial v_{k}\partial
u_{k}}&\frac{\partial^{2}f}{\partial v_{k}^{2}}.\end{pmatrix}$ (48)

###  A.1 The Amplitude Process

For $Z_{k}\begin{pmatrix}u_{k}\\\
v_{k}\end{pmatrix}=\sqrt{u_{k}^{2}+v_{k}^{2}}$ we compute,

$\nabla Z_{k}\begin{pmatrix}u_{k}\\\
v_{k}\end{pmatrix}=\frac{1}{Z_{k}}\begin{pmatrix}u_{k}\\\ v_{k}\end{pmatrix},$
(49)

and so,

$\begin{split}&(\nabla Z_{k})^{\top}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}=\frac{1}{Z_{k}}(u_{k}du_{k}+v_{k}dv_{k})\\\
&=\frac{1}{Z_{k}}\bigg{(}(-\lambda u_{k}^{2}+\omega
u_{k}v_{k}+u_{k}\big{(}\sum_{j}\mathbb{C}_{k,j}u_{j}\big{)}-\omega
u_{k}v_{k}-\lambda
v_{k}^{2}+v_{k}\sum_{j}\mathbb{C}_{k,j}v_{j})\bigg{)}dt+\frac{u_{k}(adW_{k}^{E}+bdW_{k}^{I})+v_{k}cdW_{k}^{I}}{Z_{k}}\\\
&=\bigg{(}-\lambda
Z_{k}+\frac{1}{Z_{k}}\sum_{j}\mathbb{C}_{k,j}(u_{k}u_{j}+v_{k}v_{j})\bigg{)}dt+dR_{k},\end{split}$
(50)

where

$dR_{k}:=\frac{u_{k}(adW_{k}^{E}+bdW_{k}^{I})+v_{k}cdW_{k}^{I}}{Z_{k}}=a\cos(\theta_{k}(t))dW_{k}^{E}+(b\cos(\theta_{k}(t))+c\sin(\theta_{k}(t)))dW_{k}^{I}.$
(51)

Thus,

$(\nabla Z_{k})^{\top}\begin{pmatrix}du_{k}\\\ dv_{k}\end{pmatrix}=-\lambda
Z_{k}+\sum_{j}\mathbb{C}_{k,j}Z_{j}(\cos\theta_{k}\cos\theta_{j}\\\
+\sin\theta_{k}\sin\theta_{j})dt+dR_{k}.$ (52)

Now we compute,

$\frac{1}{2}HZ_{k}\begin{pmatrix}u_{k}\\\
v_{k}\end{pmatrix}=\frac{1}{2Z_{k}}\mathbb{I}_{2}-\frac{1}{2Z_{k}^{3}}\begin{pmatrix}u_{k}^{2}&u_{k}v_{k}\\\
u_{k}v_{k}&v_{k}^{2}\end{pmatrix},$ (53)

where $\mathbb{I}_{2}$ is the $2\times 2$ identity matrix. We have, using (43)
and (44),

$\begin{pmatrix}u_{k}\\\
v_{k}\end{pmatrix}^{\top}(\frac{1}{2Z_{k}}\mathbb{I}_{2})\begin{pmatrix}u_{k}\\\
v_{k}\end{pmatrix}=\frac{1}{2Z_{k}}\big{(}(du_{k})^{2}+(dv_{k})^{2}\big{)}=\frac{a^{2}+b^{2}+c^{2}}{2Z_{k}}dt.$
(54)

Now we compute,

$\frac{1}{2Z_{k}^{3}}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}^{\top}\begin{pmatrix}u_{k}^{2}&u_{k}v_{k}\\\
u_{k}v_{k}&v_{k}^{2}\end{pmatrix}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}=\frac{1}{2Z_{k}^{3}}(u_{k}^{2}(du_{k})^{2}+v_{k}^{2}(dv_{k})^{2}+2u_{k}v_{k}du_{k}\cdot
dv_{k}).$ (55)

Using (43), (44), and (46), we have that,

$\begin{split}&\frac{1}{2Z_{k}^{3}}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}^{\top}\begin{pmatrix}u_{k}^{2}&u_{k}v_{k}\\\
u_{k}v_{k}&v_{k}^{2}\end{pmatrix}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}=\frac{(a^{2}+b^{2})u_{k}^{2}+c^{2}v_{k}^{2}+2u_{k}v_{k}bc}{2Z_{k}^{3}}dt\\\
&=\frac{(a^{2}+b^{2})\cos(\theta_{k}(t))^{2}+c^{2}\sin(\theta_{k}(t))^{2}+bc\sin(2\theta_{k}(t))}{2Z_{k}}dt.\end{split}$
(56)

So, using the substitution
$\cos(\theta_{k}-\theta_{j})=\cos\theta_{k}\cos\theta_{j}+\sin\theta_{k}\sin\theta_{j}$
the amplitude process is given by,

$dZ_{k}=\bigg{(}\frac{a^{2}+b^{2}+c^{2}-(a^{2}+b^{2})\cos(\theta_{k}(t))^{2}-c^{2}\sin(\theta_{k}(t))^{2}-bc\sin(2\theta_{k}(t))}{2Z_{k}}-\lambda
Z_{k}+\sum_{j}\mathbb{C}_{k,j}Z_{j}\cos(\theta_{k}-\theta_{j})\bigg{)}dt+dR_{k}.$
(57)

###  A.2 The Phase Process

For $\theta_{k}=\arctan\frac{v_{k}}{u_{k}}$ we compute,

$\nabla\theta_{k}=\frac{1}{Z_{k}^{2}}\begin{pmatrix}-v_{k}\\\
u_{k}\end{pmatrix},$ (58)

and so,

$(\nabla\theta_{k})^{\top}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}=\frac{1}{Z_{k}^{2}}\bigg{(}\lambda v_{k}u_{k}-\omega
v_{k}^{2}-\omega u_{k}^{2}-\lambda
u_{k}v_{k}+\sum_{j}\mathbb{C}_{k,j}(-v_{k}u_{j}+u_{k}v_{j})\bigg{)}dt+dS_{k}(t),$
(59)

where,

$dS_{k}(t):=\frac{-v_{k}(adW_{k}^{E}+bdW_{k}^{I})+u_{k}cdW_{k}^{I}}{Z_{k}^{2}}=\frac{-a\sin(\theta_{k}(t))dW_{k}^{E}+(c\cos(\theta_{k}(t))-b\sin(\theta_{k}(t)))dW_{k}^{I}}{Z_{k}}.$
(60)

And so,

$(\nabla\theta_{k})^{\top}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}=\bigg{(}-\omega+\sum_{j}\mathbb{C}_{k,j}\frac{Z_{j}}{Z_{k}}(-\sin\theta_{k}\cos\theta_{j}+\cos\theta_{k}\sin\theta_{j})\bigg{)}dt+dS_{k}(t).$
(61)

Now we compute the Hessian,

$H\theta_{k}=\frac{1}{Z_{k}^{2}}\begin{pmatrix}\frac{2u_{k}v_{k}}{Z_{k}^{2}}&1-\frac{2u_{k}^{2}}{Z_{k}^{2}}\\\
-1+\frac{2v_{k}^{2}}{Z_{k}^{2}}&\frac{-2u_{k}v_{k}}{Z_{k}^{2}},\end{pmatrix}$
(62)

and find, using (46),

$\begin{split}&\frac{1}{2}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}^{\top}H\theta_{k}\begin{pmatrix}du_{k}\\\
dv_{k}\end{pmatrix}\\\
&=\frac{1}{2Z_{k}^{2}}\Bigg{(}((du_{k})^{2}+(dv_{k})^{2})\big{[}\frac{2u_{k}v_{k}}{Z_{k}^{2}}+\frac{-2u_{k}v_{k}}{Z_{k}^{2}}\big{]}dt+du_{k}dv_{k}\big{[}1-\frac{2u_{k}^{2}}{Z_{k}^{2}}-1+\frac{2v_{k}^{2}}{Z_{k}^{2}}\big{]}\Bigg{)}\\\
&=bc\frac{v_{k}^{2}-u_{k}^{2}}{Z_{k}^{4}}dt=bc\frac{1-2\cos(\theta_{k}(t))^{2}}{Z_{k}^{2}}dt.\end{split}$
(63)

And so, using the substitution
$\sin(\theta_{j}-\theta_{k})=-\sin\theta_{k}\cos\theta_{j}+\cos\theta_{k}\sin\theta_{j}$
the phase process is given by,

$d\theta_{k}=\bigg{(}bc\frac{1-2\cos(\theta_{k}(t))^{2}}{Z_{k}^{2}}-\omega+\sum_{j}\mathbb{C}_{k,j}\frac{Z_{j}}{Z_{k}}\sin(\theta_{j}-\theta_{k})\bigg{)}dt+dS_{k}.$
(64)

###  A.3 Differences between Itô Transformations

The amplitude and phase processes given in (57) (16) and (64) (17) are
different from those derived in Greenwood and Ward (2019). This arises from
different treatments of lower order terms and of the Brownian motions
$dW_{k}^{E}$ and $dW_{k}^{I}$. Here we retained lower order terms whereas in
Greenwood and Ward (2019) they were dropped. Initially in Greenwood and Ward
(2019) and here the noise terms are expressed in vector form as
$\mathbb{E}d\mathbb{W}_{k}$, where the matrix $\mathbb{E}$ is upper-
triangular. But in the Itô transformation in Greenwood and Ward (2019), to
make the derivation simpler, it was assumed that $\mathbb{E}=\mathbb{I}$, and
thus the stochastic terms had coefficients of 1 (i.e. $a=c=1$ from (18) and
(19), whereas here we retained the coefficient $b$ (20), and did not set
$a=c=1$. We can account for the discrepancies between the Itô derivation of
Greenwood and Ward (2019) and the present one by interpreting the noises in
Greenwood and Ward (2019) as being independent with equal coefficients after
the change of basis via the matrix $\mathbb{Q}$ (13) has taken place, instead
of being independent in the original system (i.e. (8)).

It is simple to verify that the deterministic term in (57) simplifies to the
equation of the coupled amplitude processes in Greenwood and Ward (2019) in
the case where we set $a=c=1$. The noise term $dR_{k}$ in (57) can be
identified as a Brownian motion with coefficient 1. The details are
complicated but the main idea is that the increments of $R_{k}$ are given by,

$\begin{split}&R_{k}(t+\Delta t)-R(t)\\\ &=\int_{t}^{t+\Delta
t}\frac{u_{k}(t)}{Z_{k}(t)}dW_{k}^{E}(t)+\int_{t}^{t+\Delta
t}\frac{v_{k}(t)}{Z_{k}(t)}dW_{k}^{I}(t),\end{split}$ (65)

where the stochastic integrals can be written as the limits of the sums,

$\begin{split}&\sum_{i}\frac{u_{k}(t_{i})}{Z_{k}(t_{i})}\big{(}W_{k}^{E}(t_{i+1})-W_{k}^{E}(t_{i})\big{)}\\\
&+\sum_{i}\frac{v_{k}(t_{i})}{Z_{k}(t_{i})}\big{(}W_{k}^{I}(t_{i+1})-W_{k}^{I}(t_{i})\big{)}\end{split}$
(66)

where the mesh $\\{t_{i}\\}$ increases to the mesh $\mathbb{Q}\cap[t,t+\Delta
t]$. Noting that the increments in $W_{k}^{E}$ and $W_{k}^{I}$ in (66) are
both $\mathcal{N}(0,t_{i+1}-t_{i})$ (and are independent by assumption), we
have,

$\begin{split}&\frac{u_{k}(t_{i})}{Z_{k}(t_{i})}\big{(}W_{k}^{E}(t_{i+1})-W_{k}^{E}(t_{i})\big{)}\\\
&+\frac{v_{k}(t_{i})}{Z_{k}(t_{i})}\big{(}W_{k}^{I}(t_{i+1})-W_{k}^{I}(t_{i})\big{)}\\\
&\sim\mathcal{N}(0,\frac{u_{k}^{2}+v_{k}^{2}}{Z_{k}}(t_{i+1}-t_{i}))\\\
&=_{d}\mathcal{N}(0,t_{i+1}-t_{i}).\end{split}$ (67)

Hence in the limit $\\{t_{i}\\}\rightarrow\mathbb{Q}\cap[t,t+\Delta t]$ (66)
has a $\mathcal{N}(0,\Delta t)$ distribution, and so it is easy to see that
$R_{k}$ is actually a Brownian motion. Similar calculations show that taking
$a=c=1$ and $b=0$ renders (64) equivalent to the phase process equation in
Greenwood and Ward (2019).

## Competing interests

The authors declare that they have no competing interests.

## Authors’ contributions

All authors contributed to the conceptualization and writing of the paper. The
numerical simulations were accomplished by CLM.

## acknowledgements

This research was supported by grants from the Natural Sciences and
Engineering Research Council (NSERC) of Canada to Yaniv Plan (via support for
CLM), and to LMW (A9958). We would like to thank Yaniv Plan for contributing
funding despite not being directly affiliated with this paper.

## References

  * Sacerdote and Giraudo (2013) L. Sacerdote and M. T. Giraudo, in  _Stochastic Biomathematical Models_ , edited by L. N. in Mathematics (Springer, New York, 2013) pp. 99–142.
  * Izhikevich (2003) E. M. Izhikevich, IEEE Transactions on Neural Networks 14, 1569 (2003).
  * Hodgkin and Huxley (1952) A. L. Hodgkin and A. F. Huxley, Journal of Physiology 117, 500 (1952).
  * Wilson and Cowan (1972) H. R. Wilson and J. D. Cowan, Biophysical Journal 12, 831 (1972).
  * Faugeras and Inglis (2015) O. Faugeras and J. Inglis, Journal of Mathematical Biology 71, 259 (2015).
  * Kang  _et al._ (2010) K. Kang, M. Shelley, J. A. Henrie, and R. Shapley, Journal of Computational Neuroscience 29, 495 (2010).
  * Greenwood  _et al._ (2015) P. E. Greenwood, M. D. McDonnell, and L. M. Ward, Neural Computation 27, 74 (2015).
  * Wallace  _et al._ (2011) E. Wallace, M. Benayoun, W. van Drongelen, and J. D. Cowan, PLoS ONE 6 (2011).
  * Bressloff (2010) P. C. Bressloff, Physical Review E 82, 051903 (2010).
  * Baxendale and Greenwood (2011) P. H. Baxendale and P. E. Greenwood, Journal of Mathematical Biology 63, 433 (2011).
  * Greenwood  _et al._ (2016) P. E. Greenwood, M. D. McDonnell, and L. M. Ward, Journal of Coupled Systems and Multiscale Dynamics , 12 (2016).
  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer, New York, 1984).
  * Greenwood and Ward (2019) P. E. Greenwood and L. M. Ward, Mathematical Biosciences and Engineering 16, 6769 (2019).
  * Daffertshofer and van Wijk (2011) A. Daffertshofer and B. C. M. van Wijk, Frontiers in Neuroinformatics 5, 1 (2011).
  * Vogels  _et al._ (2011) T. P. Vogels, H. Sprekeler, F. Zenke, C. Clopath, and W. Gerstner, Science 334, 1569 (2011).
  * Sprekeler (2017) H. Sprekeler, Current Opinion in Neurobiology 43, 198 (2017).
  * Hellyer  _et al._ (2016) P. J. Hellyer, B. Jachs, C. Clopath, and R. Leech, Neuroimage 124, 85 (2016).
  * Richman and Moorman (2000) J. S. Richman and J. R. Moorman, American Journal of Physiology: Heart and Circulatory Physiology 278, H2039 (2000).
  * Ramdani  _et al._ (2009) S. Ramdani, F. Bouchara, and J. Lagarde, CHAOS 19, 013123 (2009).
  * Martínez-Cagigal (2018) V. Martínez-Cagigal,  _Sample Entropy_ (Mathworks, 2018).
  * Abeysuriya  _et al._ (2018) R. G. Abeysuriya, J. Hadida, S. N. Sotiropoulos, S. Jbabdi, R. Becker, B. A. E. Hunt, M. J. Brookes, and M. W. Woolrich, PLOS Computational Biology 14 (2018).
  * Kasai  _et al._ (2010) H. Kasai, M. Fukuda, S. Watanabe, A. Hayashi-Takagi, and J. Noguchi, Trends in Neurosciences 33, 121 (2010).
  * Ipsen and Peterson (2020) J. R. Ipsen and A. D. H. Peterson, Phys. Rev. E. 101, 052412 (2020).
  * Wainrib and Touboul (2013) G. Wainrib and J. Touboul, Phys. Rev. Letters 110, 1 (2013).
  * Rajan and Abbott (2006) K. Rajan and L. F. Abbott, Phys. Rev. Letters 97, 4 (2006).