{'fraction_non_alphanumeric': 0.03881946290880085, 'fraction_numerical': 0.01153948417973943, 'mean_word_length': 4.624495289367429, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Kuramoto coupling of quasi-cycle oscillators

A family of stochastic processes has quasi-cycle oscillations if the oscillations are sustained by noise. For such a family we define a Kuramoto-type coupling of both phase and amplitude processes. We find that synchronization, as measured by the phase-locking index, increases with coupling strength, and appears, for larger network sizes, to have a critical value, at which the network moves relatively abruptly from incoherence to complete synchonization as in Kuramoto couplings of fixed amplitude oscillators. We compare several aspects of the dynamics of unsynchronized and highly synchronized networks. Our motivation comes from synchronization in neural networks.


INTRODUCTION
Strong bursty oscillations observed in, for example, EEG, or in local field potentials recorded from intracranial electrodes, derive from the activity of a number of subsystems of neurons that form an interactive network. In order to understand the dynamics of such a network, we would like to have a neural system model that produces similar synchronized oscillations.
In [1] we considered a stochastic neural system comprised of a group of pairs of excitatory and inhibitory neurons (E-I pairs). We described it in terms of a pair of interacting stochastic differential equations and found, using a new result about quasi-cycles, that such a model produces bursts of narrow-band oscillations, for example bursts of gamma-frequency oscillations.
The present paper is the first to study coupling of quasicycles. Our aim is to understand, for quasi-cycles, how the degree of synchronization depends on the strength of coupling of the subsystems.
Previous work has shown that groups of neurons with group potential oscillating independently can become synchronized when driven by a common input [2,3]. The paper [4] shows that both noisy limit cycles and quasi-cycles can arise from the Wilson-Cowan model, which is similar to the model we use to generate quasi-cycles. In [5] a Kuramoto-type coupling is applied to deterministic Wilson-Cowan model limit cycles. In [6] Kuramoto-type coupling In contrast, here we omit the sigmoidal function in the initial model so that its centered dynamics for the model without noise would be simply oscillations damped to a fixed point at (0,0). We use a recent stochastic dynamics result of Baxendale and Greenwood [7] to move to an equivalent model where amplitude and phase of quasi-cycles can be conveniently viewed as the primary variables satisfying a system of stochastic differential equations (SDEs). Then we introduce Kuramoto-type coupling to the amplitude and phase equations separately. The resulting model tells us how the degree of coupling affects the degree of synchronization among quasi-cycles.
Our results apply in the following general context: consider a system of stochastic oscillators, where each subsystem is described by a pair of linear, or locally linear, SDEs. Suppose the deterministic version of this system, with zero noise, has oscillations damped to a fixed point. Then, each stochastic subsystem has sustained oscillations. We coupled these stochastic subsystems and studied their synchronization. We first write a linear model similar to that of [8], which in neuroscience is often called a Wilson-Cowan-type model. The Excitatory-Inhibitory (E-I) model for each subset of neurons (E-I pair or stochastic oscillator) appears as Here W E , W I are independent, standard Brownian motions. The parameters S EE , S II , S IE , S EI ≥ 0, are constants that represent the mean efficacies of the excitatory or inhibitory synaptic connections to post-synaptic neurons within each separate population, as indicated by the notation, with S IE representing input to inhibitory from excitatory neurons. These parameters, along with the time constants, τ E , τ I , and amplitudes of the Brownian motions, σ E , σ I , determine the oscillatory behaviour of the system and in particular its resonant frequency of oscillation. We limit our discussion here to the parameter ranges where the oscillation is narrow-band and thus has a distinct phase even though it arises from a stochastic process [1]. Note that (1) can be interpreted as applying to a single pair of neurons, or even to the subthreshold dynamics of a single stochastic neuron [9,10]. Here, however, we consider that as a single oscillatory system characterized by a particular resonant frequency, regardless of how many neurons are involved in creating that system. Note also that we can consider (1) to be a linear approximation, in the neighborhood of a fixed point, to a more elaborate dynamical system. System (1) has already been centered at the fixed point, (0,0), and can be written as what recently has been termed a linear noise model (or linear noise approximation if the original model is not linear), and When we take σ E = 0 and σ I = 0, the deterministic system obtained from Eqn. (2) has oscillations that damp to the stable point (0,0) at a rate λ. In fact, Cowan's early work [11] used the fact that asymmetric coupling of neurons would typically generate oscillations because of the presence of complex eigenvalues near the steady state. Under certain conditions the stochastic system in (2), with elements of N nonzero, has sustained oscillations, called quasi-cycles, of a narrowband nature [7,12]. This is the case when the eigenvalues of the matrix −A are complex, −λ ± iω d , with 0 < λ ω d . Here, as in [8], λ is the damping rate of the oscillation, and ω d is its natural frequency, which is positive when the eigenvalues are complex. We give the name quasi-cycle oscillator to a linear noise model (2) such that λ and ω d satisfy 0 < λ << ω d , and hence has oscillations sustained by noise.

Background II: An Asymptotically Equivalent Amplitude and Phase Model
When a system (2) is a quasi-cycle oscillator, i.e., the damping rate is substantially less than the frequency of the oscillations, and the eigenvalues of −A are complex with negative real part, the damped oscillations of the deterministic version of the system will be sustained by 'small' noise. In a neural system, this noise is always present because of fluctuating inputs occurring at synapses and because of the stochastic effects of ion channels [13]. For reasonable parameter values in (1), we have 0 < λ ω d , i.e. the range of parameter values is such that sustained oscillations are produced. This choice is justified in, e.g., [1] and [8]. Under these conditions, and according to a limit theorem of Baxendale and Greenwood [7], the process (2) is approximated as λ/ω d approaches 0, by where R is the rotation Q is a matrix that transforms A in Eqn.
(2) into a canonical form S(t) is a standard two-dimensional Ornstein-Uhlenbeck process with independent components, and is a scalar.

INTRODUCTION
The approximation we used to obtain (7) from the linear model (2) is very different from the familiar method involving expansion and dropping of higher order terms. Thus we here sketch a proof of the result leading to (7) starting from the linear noise approximation (2), keeping in mind that 0 < λ ω d will be used in an essential way [7]. We make three changes of variables in (2). First we transform the matrix −A to normal form in order to see more clearly the separate effects of the relatively slow damping characterized by λ, and the relatively fast rotation characterized by ω. Using a matrix Q as in (9) (see Eqn. (4.4) in Then where C = Q −1 N. The noise in (11) has covariance matrix Next we write where R s is the rotation (8). Using the SDE (11) and Itô's formula we obtain Finally, in order to compare the process Z(t) with a standard two-dimensional Ornstein-Uhlenbeck process, we rescale time and space by writing where To identify the SDE that U(t) satisfies it is convenient to use the integrated form whereW(t) = √ λ W(t/λ) is a standard Brownian motion. Therefore where D = (1/σ)C, so that tr(DD ) = 2. Theorem 1 of [7] says that the distribution of {U λ (t) : 0 ≤ t ≤ T }, where U λ denotes the process satisfying (18) with a fixed λ, ω, U λ (0) = x, fixed, and 0 < T fixed, converges, as λ/ω goes to 0, to the distribution of the standard 2-dimensional O-U process {S(t) : 0 ≤ t ≤ T } with S(0) = x and independent components. The proof uses the Martingale problem method (see [7]).
The approximation (7) is then obtained by reversing the three changes of variables, starting from the SDE for the standard 2-dimensional O-U process In order to obtain an error bound on the approximation one would need a uniformity result with T → ∞ as λ/ω → 0, in Theorem 1 of [7]. In [1] we applied this result to the system (2) and showed that the approximation (7) is quite good for the parameters we considered. The approximation predicts several critical quantities: the constant ratio of envelopes and phase difference of the E and I processes, their average instantaneous frequencies, and the relationship between the variance of instantaneous frequency and the amplitude for each of the two processes. Simulations of the full model defined by (1) and (2), using the same parameters used later in this paper, showed that these were all very close to the respective values predicted from the approximation. Impressively, the inverse relationship between instantaneous frequency variance and the envelope amplitude was just as predicted from the approximation (22), with large phase perturbations occurring more frequently when envelope amplitude was small. Moreover, as we show in the Results section, the approximation produces sample paths that qualitatively agree with those produced by the full model, as well as displaying the relationship between dφ and Z(t) predicted by (22). In terms of polar coordinates the approximation (7) can be written as where |S(λt)| = S 1 (λt) 2 + S 2 (λt) 2 is the amplitude of S(λt), and For a standard two-dimensional Ornstein-Uhlenbeck process S(t), the stochastic process literature provides an explicit relationship between the process Z(t) = |S(t)|, the radial process or modulus process associated with S(t), and φ(t), the phase perturbation process (see [14], [15]). The pair of processes, Z(t) and φ(t), satisfies a system of stochastic differential equations where W (t) is a standard Brownian motion on the real line and b(t) is a standard Brownian motion on the unit circle, independent of W (t). We will see in the Results section that the approximation, Equation (7), generates the same types of sample paths as the original model. Thus, from this point onward, we will take the approximation (20) as our model.

EXTENSION TO A POPULATION
Suppose we have a family of N E-I subpopulation models indexed by i = 1, 2, ...N . We will have a version of (1), and hence of (20), for each member of the family, each with its own parameter values. The σ E , σ I parameters for E and I subpopulations are all assumed to be equal for now, as are the time constants τ E , τ I . This can of course be changed to include distributions of parameters for the different subpopulations. These have been chosen to be biologically realistic insofar as possible. From (20), (21) and (22) we have, and for the amplitude of where ||Q|| is the 2-norm of Q. An advantage of V * (t) over V(t) is that when we transform to polar coordinates as in (20) we have the explicit expressions (23) and (25) for the phase and amplitude processes, which are not available for V(t), and which are a natural way to describe quasi-cycle oscillators. This allows us to couple our family of quasi-cycle oscillators using Kuramoto-type coupling.

KURAMOTO COUPLING
The history of the Kuramoto model is beautifully set out in [16]. See also [17] for a review. The model describes a dynamic system of phase changes that produces phase synchronization among a population of oscillators, with phase functions θ j , j = 1, 2, ...N , depending on a distribution g(ω) of natural frequencies, ω j , and functions Γ j,k of phase differences: The interaction functions Γ j,k could involve Fourier harmonics and an unspecified connection topology. In fact, Kuramoto studied the mean-field version of the model with sinusoidal coupling, For a collection of points e iθj moving about on the unit circle in the complex plane, the quantity is the centroid of the phases, θ j , where the radius, ρ(t), called the phase locking index, measures phase coherence, and ψ(t) is the average phase. If the phases, θ j , j = 1..., N, move in a clump, the index ρ(t) is near 1, whereas if the phases are spread out around the circle, the index ρ(t) is near zero. One can write (28) as for any θ j . Equating the imaginary parts and substituting into (27) one obtains In this form, we see that each phase θ j is pulled toward the mean phase, ψ, with a strength proportional to the phase locking index ρ. If the population is becoming more coherent as time advances, ρ(t) increases and so does Kρ(t).
Using the feedback loop, Kuramoto derived a critical value is a saturating level of coherence [16].

Kuramoto Coupling of a Family of Quasi-cycle Oscillators
Results for stochastic versions of the Kuramoto model have been obtained, e.g., by [18,19]. Our question here is whether these results for the Kuramoto model (27), in particular monotonicity of coherence with coupling strength and rapid change of coherence near a critical value of coupling strength, will continue to hold for a coupled version of the stochastic systems (23) and (25) for i = 1, ..., N . We modified the system (23) and (25) by inserting a coupling term like that in (26) into their drift coefficients, assuming Theorem 1 of [7] holds uniformly over i = 1, ..., N . The factor Z j /Z i , the ratio of amplitudes, appears in the coupling of phases as in [5,6]. It reconciles the fact that subpopulations i and j have different amplitude processes, and their ratio is a separate factor from the coupling strength determining the impact of the phase subpopulation j on the phase of subpopulation i. This point is pursued in the Discussion. The result is a system of coupled phases, where C is a matrix defining the coupling between subpop- In addition to impacting each others' phases, we would expect that the amplitudes of the various E-I subpopulations would impact each other as well. This is important because the amplitudes affect the phase slip processes (24). As in an OU process, where the process is pulled toward 0 in proportion to its distance from 0, we can conceptualize the driving impact of each subpopulation's amplitude on that of the others as a 'pull' on the others' amplitudes proportional to the difference between them. This symmetrical coupling has the advantage that it doesn't allow the amplitudes to become too large, which is necessary in order that the phase slip processes (24) be observable. We recognize that other forms of amplitude coupling are possible and that some of these might be more realistic in neural systems. We come back to this point in the Discussion. We couple the amplitudes (25) as follows: The overall degree of coupling in the entire set of subpopulations can be characterized in general by the 2-norm of the matrix C, ||C||, which when normalized with regard to the population size, N , can be seen to play the same role as the coupling coefficient, K, in the Kuramoto model.

RESULTS
In order to study the consequences of this Kuramoto-type coupling of quasi-cycle oscillators, we solved (29) and (30) numerically using the Euler-Maruyama method [20] for a range of the number of coupled E-I systems, N , and the overall coupling strength between the systems, ||C||, similar to Kuramoto's K. We assumed that the subpopulations all oscillated in a narrow band of natural frequencies distributed normally with mean 437.72 rad/s and standard deviation 1 rad/s, similar to the standard system studied for the Kuramoto model. We also began each realization with the phases of the various systems distributed uniformly between −π and π. Finally we studied all-to-all coupling matrices C i,j , all equal except that C i,i = 0, ∀i. We investigated the coherence measure, ρ, as a function of ||C||, and also studied the dynamics of the coherence by looking at the membership of the synchronous and asynchronous E-I groups via raster plots and plots of the mean frequency, amplitude, and phase in the most synchronous group compared with those of the entire set of processes for representative values of ||C||.
In our simulations we chose a natural, or resonant, frequency, ω d , that for different E-I processes was chosen from a (clipped, i.e. compact support) Gaussian distribution. The damping rate, λ and the noise amplitude, σ, also were required to vary in a consistent way because both quantities depend, in different ways, upon the same parameters from (1), namely the synaptic efficacies S i,j , i, j = E, I, the time constants τ E , τ I , and σ E , σ I (see (5), (6), (10)). For purposes of the present study, λ was computed for each separate process from the following equations given ω d for that process (many significant digits retained because of the sensitivity of the calculations), For the parameter values we used in our simulations of (30), ||Q|| ≈ 703.5. We rescaled the amplitudes into a range that allows for an observable phase slip process by using ||Q||/703.5 ≈ 1 instead of ||Q|| in those simulations. As mentioned earlier, if the amplitudes become too large, the phase slip processes (24) become negligible and we have what amounts to a standard Kuramoto model.

Uncoupled oscillations
Before we describe our results for the coupled systems defined by (29) and (30), we need to establish that the phase   (7).The natural frequency, ω, of the E-I process was set to 69.66 Hz (437.72 rad/s). Other parameters as in Table 1. and amplitude processes, (23) and (25) produce appropriate stochastic oscillatory paths. To this end we first simulated (7) to demonstrate the ability of the approximation V * (t) to generate oscillatory sample paths qualitatively similar to those generated by the full model (1). One such sample path is shown in Figure 1 and it is qualitatively the same as the sample paths produced by the full model, including the gamma bursts (episodes of high amplitude oscillations, cf. Figure 6 in [1]) produced by that model. We then simulated (29) and (30) in the absence of coupling (C = O (the zero matrix)), which are in fact the same system as (23) and (25), and inspected the paths of the individual uncoupled phase and amplitude processes. Figure 2 displays a typical phase path, and Figure 3 a typical amplitude path for a single E-I process approximated by (29) and (30). The phase and amplitude processes of Figures 2 and 3 are not those of the model (6) of [1], which is (1) of the present paper, but are samples from the almost equivalent model V * (t) of (20). The oscillatory nature (rotation in polar coordinates) of this process is revealed by the orderly phase progression. The amplitude process can be regarded as the envelope of the oscillations.  Note in Figure 2, however, that the orderly phase progression is interrupted on occasion by large increments of the phase perturbation process. These increments are related to the amplitude process as predicted by Equation (22). An example of this is shown in Figure 4 where phase is plotted with amplitude in the same realization. Large phase perturbation increments occur when amplitude is low and small increments occur when amplitude is large.  4: Amplitude, Z(t) and phase, φ(t), of a single process plotted together. The data are from the same realization as that in Figure 2 and Figure 3. Notice that large phase slips occur when the amplitude is low, around 1.03 s and 1.25 s in this figure.

Phase coherence and coupling
For the purposes of this paper we define synchronization in terms of phase coherence, which is the Kuramoto Phase Locking Index (PLI) in (28), The θ j are the phases of the N separate stochastic processes V * j (t). Figure 5 displays the average phase coherence of 10 independent realizations of (29) and (30) as a function of coupling strength, ||C||, among E-I populations with N =2, 10, 66, and 100, with (clipped between -3 and +3) normallydistributed natural frequencies. As in the basic and stochastic Kuramoto models, phase coherence is nearly constant for coupling strengths below a critical value, and increases rapidly (for N ≥ 10) with increasing coupling strength until near an asymptote at which the phases are all locked together and ρ stays very near 1, with only occasional slight departures caused by occasional large noise samples in one or more of the coupled processes. The standard deviation of ρ over the 10 realizations in each graph is roughly constant for small coupling values except near the critical coupling value, where phase coherence is rising rapidly under the mutual influence of the different E-I processes. The standard deviation of ρ, of course, decreases dramatically after PLI reaches its asymptote. The effect of increasing the number of coupled systems is to increase the critical coupling value, in terms of ||C||, and to decrease the minimum phase coherence of the system. Clearly stronger coupling is required to synchronize a larger number of independent systems.
FIG. 5: Phase coherence, ρ (phase-locking index, PLI), as a function of coupling strength, ||C||, with normally distributed natural frequencies. In all graphs, the connected points represent means over 10 independent realizations ofρ, and the solid lines at the bottom of the graphs are the standard deviations ofρ over the 10 realizations. For each realization,ρ entered into the interation average was the average value of ρ over 5000 iterations (0.25 s) of the model, beginning after iteration 5000, at which point the model was at or close to a steady state. Natural frequencies, ω d , of E-I processes were sampled from a (clipped) Gaussian distribution with a mean of 437.72 rad/s and a standard deviation of 1 rad/s, and this was done anew for each iteration and each coupling value. Starting phases were randomly distributed between −π and π, and starting amplitudes were uniformly distributed between 0 and 1.

Dynamics of phase coherence
How do the group dynamics of the unsynchronized and synchronized collections of the stochastic E-I pairs evolve, and how do they compare? To gain insight, we determined, using simulations similar to those in Figure 5, for two coupling values, ||C|| = 0 and ||C|| = 4950, and with N = 100, which particular pairs of E-I processes, characterized by their natural frequencies, participated in the group of processes that was most highly synchronized. To do this we sorted the phases of the processes into 20 equal phase bins between −π and π and chose the one with the largest membership at each iteration as the 'synchronous group' at that moment. Figures 6 and 7  the synchronous group at each time point of the simulation compared with those of the total set of 100 processes.
The raster plots of synchronous group membership clearly show that when coupling is 0 (Fig. 6), membership is highly variable and the synchronous group itself is sporadically populated, with most of the processes involved at some point regardless of natural frequency. There is no strong tendency for either extreme or central frequencies to dominate. For stronger coupling (Fig. 7), the result is similar, except that more processes are in the synchronous group more of the time, and the processes with higher natural frequencies dominate the synchronous group. Apparently the independent noises for phase and amplitude for each process preclude most processes from settling into the synchronous group for a very long time, although a few do seem to be in it most of the time, including especially those with higher natural frequencies which would tend to lead the others in phase.
In light of the heterogeneity of synchronous group mem- Coupled processes -substantial synchronization. Relatively strong coupling between processes, single realization with N = 100, ||Ci,j|| = 4950. In the raster plot in the upper left each plotted point represents membership of the indicated process in the 'synchronous group' (see text) at the indicated iteration. Red points represent processes whose natural frequency fell within the extreme two highest and two lowest of 20 equal-sized bins over the range of the normal distribution of natural frequencies; blue points represent processes whose natural frequencies fell within the two most populous bins; green points represent the remainder of processes, with natural frequencies between the extremes and the most populous bins. Upper right plot: PLI over the 100 processes at each time point. Lower left plot: Red line indicates mean phase (theta) of all 100 processes; blue line indicates mean phase of processes in the synchronous group. Lower middle plot: Red line indicates mean natural frequency of all 100 processes; blue line indicates mean natural frequency of processes in synchronous group. Note that instantaneous frequencies are not always equal to natural frequencies because of phase slips caused by driving noise. Lower right plot: Red line indicates mean amplitude of all 100 processes; blue line indicates mean amplitude of the processes in the synchronous group. bership across model time, the plots of mean frequency and mean amplitude of the processes in the synchronous group with stronger coupling are interesting. Figure 7 reveals that, on average, when coupling is relatively strong the mean natural frequency of the synchronous group is higher than that of the entire population. Figure 7 also shows that the mean amplitude of the processes in the synchronous group is higher than the mean amplitude of the entire population. It is reasonable that processes with higher amplitudes and frequencies would tend to dominate the interactions. Figure 8 displays additional detail concerning the evolution of synchronization in the case of the realization shown in Figure 7. The number of members in each phase bin is plotted in Figure 8 for a series of 0.00005 second time slices. The 'synchronous group' is represented by the highest peak of the distribution at each time point. Clearly, early in the realization, the phases of the various processes are incoherent, spread across the 20 bins roughly equally as they were at the onset of the realization. Gradually, as model  Figure 7) of the model for relatively strong coupling between processes, showing the progression of the neural oscillators from an unsynchronized to a synchronized state. Each plot represents a 5 ms period of the realization. Plots from upper left to lower right progress from early to later in the realization; specific time points are indicated by numbers that match the red numbers in the upper right plot in Figure  7, which also plots the associated coherence values.
time progresses, the interactions between the processes promoted by the moderately strong coupling creates first several clumps of synchronized processes, and then a single dominant group. Notice that, consistent with the slight drop in PLI later in the realization, at around 0.4 seconds (labeled '6' in the graph of PLI in Fig. 7), the dominant group becomes somewhat less coherent again, presumably under the influence of the noise that is driving the quasi-cycles.

DISCUSSION
We have reported results from a combination of a linear stochastic model of populations of neurons with the Kuramoto approach to synchronization of quasi-cycle oscillators. This is new because for the first time we demonstrate for quasi-cycles the existence, and study the properties, of a critical coupling leading to their synchronization. Such oscillators are likely to be representative of a broad class of biological oscillators, such as neurons and neural networks, and including gene networks. The E-I processes in our model could represent a range of possible neurological configurations, from pairs of reciprocally-interacting neurons, through populations of sparsely-firing neurons that generate an oscillating field potential, to synchronously firing populations of neurons.
We have demonstrated that for this novel quasi-cycle stochastic Kuramoto model for finite N , as for other Kuramoto models, phase coherence generally increases with coupling strength. The incoherent state at low coupling changes first to a partially coherent state and then to a nearly phase-locked state as coupling increases beyond a critical value. We have also shown some of the details of the transition in the form of raster plots and other properties of the group of most nearly phase-locked oscillators. The exploration of additional coupling matrices is left to future work. If quasi-cycles are indeed the source of neural oscillations, as argued in [1], then the results in the present paper become relevant to synchronization of neural systems in real brains.
In the present work we introduced, in addition to phase coupling, coupling of the amplitudes of the quasi-cycle oscillators. This was necessary because the amplitudes of the oscillators affect their phase slip processes, and thus also affect the synchronization of the phases. If phase slips are large and frequent (i.e., when amplitudes are very small), synchronization becomes more problematic, whereas if phase slips are small and infrequent (i.e., when amplitudes are very large) the rotation of each quasi-cycle is more regular and the situation more closely resembles interactions among a set of pure phase oscillators. One limitation of our present work is that the form of the amplitude coupling we used was symmetrical, i.e., the lower amplitude oscillators pulled the higher ones down, and vice versa, so that the amplitudes tended to stabilize in a small range that allowed for observable phase slips. We also scaled the amplitudes to maintain this range. It is possible that other forms of amplitude coupling would operate similarly. For example, we tested a model with the ratio of amplitudes, Z j (λ j t)/Z i (λ i t) inserted in place of the difference of amplitudes in (30). Results were highly similar, although the amplitudes tended to stabilize at somewhat higher levels. On the other hand, inserting only the amplitude of the 'other' coupled oscillator, Z j (λ j t), in place of its difference from Z i (λ i t) results in an explosion of amplitudes and thus an exit from the quasi-cycle regime. Similarly, a coupling factor of cos(θ j − θ i ) or Z j (λ j t)cos(θ j − θ i ) also has this effect. Thus, these couplings will not work for quasi-cycle oscillators.
All of this raises the question of what a biologicallyrealistic coupling of quasi-cycle oscillators would be, for example in the context of a neural network. Typically, neurons are coupled via synapses. Our original model (1) expresses the coupling between E and I neurons realistically as a bidirectional, short-range synaptic coupling. But when we move into the realm of the quasi-cycle oscillator we lose the ability to represent directly the long-range couplings between the neurons themselves, and can only represent couplings between the quasi-cycle oscillators' phases and amplitudes. It remains to be worked out how the couplings we have studied are related to realistic couplings between neurons or other biological entities. Long range couplings between neurons are more common between excitatory neu-rons, but can also occur between I and E neurons. Thus, the coupling we studied in the present work might best capture a network of quasi-cycle oscillators in which bidirectional E and I connections occur at two scales, say within micro-columns (the individual E-I quasi-cycle oscillators) and also between columns (the coupled set of quasi-cycle oscillators). Such a network would allow for both amplitude and phase 'pulls' in both directions, as in (29) and (30). Our conclusion would then be that the larger-scale set of oscillators synchronizes according to coupling strength in a Kuramoto-like manner.

Relation to other models
Although we have not provided mathematical results for the critical coupling of quasi-cycles, there do exist such results for related models. For example, the paper of Ott and Antonsen [21] considers what happens when a system of coupled oscillators is also coupled to an external oscillator. It is found, using an Ansatz regarding the form of the coefficients in the Fourier expansion of the Fokker-Planck (Kolmogorov) equation for the stochastic Kuramoto system, that in the large N limit there is a critical value of the coupling constant at which there is a phase transition from incoherence to partial or complete synchronization.
Sonnenschein and Shimansky-Geier [18] derived a similar result for a stochastic model with Gaussian noise, assuming that the phases of the oscillators are Gaussian in distribution. In their model, the variance of the phases of the oscillators, σ 2 (t) in their notation, goes to infinity as t → ∞ for C < 2D, and has a finite limit for C > 2D, where C is the coupling strength and D is the variance of the noise. It is quite possible that a similar balance between coupling strength and noise variance holds in our model, (29) and (30), although our noise variance in (29) involves Z i . Moreover, 2D in their model is a transition point between two regimes of phase variance, whereas in our model there is a critical coupling strength that separates incoherence from coherence. Nonetheless, if we consider a normalized coupling matrix 2-norm, ||C||/(N 2 − N ), as the coupling strength, and the N = 100 case in which ρ is near 1 at ||C|| ≈ 10 4 , then the critical coupling value is around 10 4 /(10 4 − 100) ≈ 1 in our model. Because Z i (λ i t) converges to values around 1 when coupling is large (or even moderate, see Figure 7), and the variance of the noise is 1 in our simulations, at convergence an analogue of the D of [18], (db(λ i t)/Z i (λ i t)) 2 , is approximately equal to 1/1 = 1. Their theory would predict a value of 2 × D = 2 × 1 = 2, for our critical coupling value, which is of the same order of magnitude as what we observed, although this could certainly change as N becomes even larger.
It should be noted, however, that our approximation informs us of an additional constraint on the relation between critical coupling value and the noise that is driving the quasi-cycles. Namely, σ/ √ λ must not go to either 0 or infinity. Because λ must be small relative to ω, σ too must be relatively small in order for the quasi-cycles to be of at least moderate size. Thus, because the noise is constrained, if there is a relationship between the critical coupling value and the noise, that relationship is also constrained.
The recent paper of Ton et al. [6] seems to be one of the closest to ours in some respects. The deterministic dynamics, however, is taken from a Freeman neural mass model expressed by a second order ODE rather than from a linear stochastic model that generates quasi-cycles; E-I pairs are involved but in a different way. The emphasis is on the effect of delays, which are transformed into phase shifts. Deterministic limit cycles around an unstable focus are the source of unstable oscillations, rather than noisy quasi-cycles near a stable focus as in the present paper. A commonality is that the limit cycles in [6] have a variety of radii and each coupling parameter contains a factor of the ratio of these radii, just as in [5], with justification involving, again, averaging over one period of the oscillation (see [6], Methods). We follow this procedure in introducing the ratio of amplitudes into (29), where the distances of the quasi-cycle oscillations from the fixed point vary from moment to monent for each oscillator as well as across oscillators at any given moment. The introduction of such a factor into the phase coupling is not critical in our model, however, for we found similar results when it was omitted and the coupling was exactly as in the Kuramoto model.
The object of study in [6] is the dynamics of the probability density of the vector of phases, stochasticity arising only from the distribution of the natural frequencies. They find, for fixed N , that there are stationary distributions of the vector of phases clustered around a particular vector. We speculate that results of this type might also hold for our very different model, (29) and (30), which is based on quasi-cycles.
There is a significant mathematical literature, not focused on neuron populations, about Kuramoto-type synchronization in the presence of noise. Two such recent papers are by Giacomin et al. [22] and by deVille [19]. Giacomin et al. studied long-term dynamics of infinitely-many identical, noisy, phase oscillators using the Kolmogorov equation of the stochastic system. When coupling strength is below a critical value there is an incoherent state, e.g., phases are uniformly distributed. When coupling strength is above this critical value the global attractor includes a partially synchronized state. One would conjecture that similar results would hold for infinitely-many interacting quasicycles. Even the careful statement of an analogue to the Giacomin et al. results, however, is beyond the scope of this paper.
Finally, deVille [19] studied the long-term behaviour of a Kuramoto system of oscillators perturbed by noise using Friedlin-Wentzell theory. That work concentrates on, but is not limited to, nearest-neighbour interactions. This viewpoint is very different from ours, since we begin with a