Synchronization dynamics of two coupled neural oscillators receiving shared and unshared noisy stimuli

The response of neurons to external stimuli greatly depends on the intrinsic dynamics of the network. Here, the intrinsic dynamics are modeled as coupling and the external input is modeled as shared and unshared noise. We assume the neurons are repetitively firing action potentials (i.e., neural oscillators), are weakly and identically coupled, and the external noise is weak. Shared noise can induce bistability between the synchronous and anti-phase states even though the anti-phase state is the only stable state in the absence of noise. We study the Fokker-Planck equation of the system and perform an asymptotic reduction ρ ₀. The ρ ₀ solution is more computationally efficient than both the Monte Carlo simulations and the 2D Fokker-Planck solver, and agrees remarkably well with the full system with weak noise and weak coupling. With moderate noise and coupling, ρ ₀ is still qualitatively correct despite the small noise and coupling assumption in the asymptotic reduction. Our phase model accurately predicts the behavior of a realistic synaptically coupled Morris-Lecar system.

1. Z(t) is also called the infinitesimal phase-response curve (PRC) and the adjoint of the limit cycle because it is a periodic solution to \(\frac{dZ}{dt}=-D_XF(X_0(t))^T Z(t)\). A negative PRC value is a delay in phase upon receiving a perturbation at time t, while a positive PRC value is a phase advance.

2. Necessary and sufficient conditions for Schur stability (i.e., all roots have modulus less than 1) are quite complicated, but they are stated in Hinrichsen and Pritchard (2005).

We thank Brent Doiron for useful discussions. CL is supported by an NSF Postdoctoral Fellowship # DMS0703502. GBE is supported by an NSF grant # DMS0513500.

Authors and Affiliations

1. Department of Mathematics, University of Pittsburgh, Pittsburgh, PA, 15260, USA

   Cheng Ly

2. University Professor of Computational Biology, Professor, Department of Mathematics, University of Pittsburgh, Pittsburgh, PA, 15260, USA

   G. Bard Ermentrout

Corresponding author

Correspondence to Cheng Ly.

Action Editor: David Golomb

Appendix A: Solving for a _n is a well-posed problem for a particular PRC and coupling function

We seek bounded solutions to the system (18):

for n ≥ 1, given a ₀. Recall that a _n  = a _− n for all n.

1.1 The case with no coupling (K = 0)

The recursive relation with coupling is non-autonomous, but is autonomous without coupling. Setting K = 0, the system (18) is:

The general solution is:

where a ₀ and a ₁ must be specified by initial conditions. a ₀ is already determined (\(\frac{1}{4\pi^2}\)) so we are only left with one degree of freedom a ₁. The a _n terms must go to 0 because we seek bounded solutions (this is also necessary for ρ ₀ to be a convergent Fourier series), thus k ₂ must vanish (because \(\frac{1}{c}+\sqrt{\frac{1}{c^2}-1}>1\) for 0 < c < 1), which means a ₁ is can no longer be freely chosen. The constant k ₁ can be solved for since we know a ₀, therefore, the unique solution is:

which is in agreement with our numerical implementation (plots not shown) where we truncate the system at a large enough N and solve a matrix-vector equation (see Section 3.2). Notice the rate of decay of a _n to 0 only depends on the amount of shared noise c in this case; a _n goes to 0 faster as c gets smaller. For the degenerate case of only shared inputs c = 1, a _n  = a ₀ for all n (i.e., ρ(ϕ) = δ(ϕ)).

1.2 The case with coupling (K > 0)

This system at hand is (for n ≥ 1):

For n > > 1, the recursive relation nearly reduces to the case with no coupling (K = 0) above. We now state a useful theorem due to Hinrichsen and Pritchard (2005). Let

Theorem 1

(Hinrichsen & Pritchard) Let \(p(\lambda,\tilde{a})\) be a polynomial of degree \(\tilde{n}\) with distinct roots r _j each of multiplicity m _j and q(λ,a) be a polynomial of degree \(n\geq\tilde{n}\) . Then for any ε > 0 so that the closed discs \(\overline{ B_{\varepsilon}(r_j) }\) are mutually disjoint, there exists a δ(ε) > 0 such that for all \({\bf a}\in D_{\delta(\varepsilon)}(\tilde{a})\) , there are m _j roots of q(λ,a) in B _ε (r _j ) and the other \(n-\tilde{n}\) roots lie outside of \(B_{\varepsilon^{-1}}(0)\).

The general solution to Eq. (30) is

where a ₀, a ₁, a ₂ and a ₃ must be specified initial conditions. For large enough n,

where λ _j are the roots to the polynomial:

As n→ ∞ (ε→0), by Theorem 1, three of the roots are arbitrarily close to: 0, \(\frac{1}{c}\pm\sqrt{\frac{1}{c^2}-1}\) (i.e., the roots of cλ ³ − 2λ ² + cλ), while the 4^th root goes to ∞ (in modulus). Since the solution must be bounded, two of the four k _j (say k ₃ and k ₄ without loss of generality) must vanish (corresponding to the |λ ₄|→ ∞ and \(\lambda_3\approx\frac{1}{c}+\sqrt{\frac{1}{c^2}-1}\)). Note that λ ₁ is arbitrarily close to 0 but never equal to 0 as long as K > 0 and for n large enough, \(\lambda_2\approx\frac{1}{c}-\sqrt{\frac{1}{c^2}-1}\implies|\lambda_2|<1\). Equation (30) for n = 1, has three unknowns (a ₁, a ₂, and a ₃), keeping in mind that a _− 1 = a ₁ and that a ₀ is given. But a ₂ and a ₃ (without loss of generality) cannot be chosen arbitrarily because k ₃ = 0 and k ₄ = 0. So the n = 1 equation determines a ₁. Therefore, there is a unique solution.

Disregarding the transient behavior for intermediate n, we know for large enough n, the solution is asymptotically

and a _n →0 as n→ ∞.

Appendix B: Sufficient conditions for well-posedness: solving for a _n with general coupling H and PRC Δ

In general, the infinite system Eq. (26) for a _n with the requirement that a _n be bounded may not be well posed. However, for reasonable coupling functions H and PRCs Δ, the arguments in Appendix A give insight to our problem. Recall that our system is (for n ≥ 1):

1.1 Coupling and PRC with finite Fourier series expansion

Assume the Fourier series expansion of H and Δ have at most 2M + 1 terms. Then the system Eq. (26) is a 2M+1 term inhomogeneous recurrence relation. The first M − 1 equations (n = 1,2,..,M − 1) have to be altered appropriately because a _n  = a _− n . The general solution is of the form:

where the initial conditions a ₀, a ₁, ..., a _2M − 1 are specified.

For n > M, the right hand side of the system is 0:

Let \({\tilde d}_m=d_{-m}d_m\). When n > > 1, the coupling term are arbitrarily small, and our system is asymptotically close to the recurrence relation (by Theorem 1):

Thus,

where λ _j are the roots of the polynomial:

For bounded solutions, a sufficient condition is that the roots of Eq. (33) have modulus^{Footnote 2} less than 1: |λ _j | < 1. To have a unique bounded solution (well posed), M of the roots have to have modulus greater than 1 and consequently, the other M roots have modulus less than 1. The M roots that blow up have to have coefficients (k _j ) vanish. We are left with M initial conditions a ₀, a ₁,..., a _M − 1 (note that a _M can be uniquely determined in n = 1 equation once the initial conditions are specified). These are chosen to give bounded solutions. By Theorem 1, a _n →0 as n→ ∞.

1.2 Coupling and PRC with infinite Fourier series expansion

The Fourier series of H and Δ must converge, thus the terms must decay fast enough to 0:

For n large enough, the system is asymptotically close to Eq. (32):

The large n behavior is a recurrence relation with infinitely many terms. Since many of the coefficients of the infinite recursive relation are very small, the systems behaves like that of the previous section. The most general conditions for well-posedness is a difficult question that is beyond the scope of this paper.

Appendix C: Morris-Lecar equations

The two Morris-Lecar cells are:

with \(\left<\xi_1(t),\xi_1(t')\right>=\delta(t-t')\), \(\left<\xi_2(t),\xi_2(t')\right>=\delta(t-t')\), and \(\left<\xi_1(t),\xi_2(t')\right>=c\delta(t-t')\), c ∈ [0,1]. The auxiliary functions are:

The parameter values used in the figures are: \(C=20\frac{\mu\text{F}}{\text{cm}^2}\), \(I_{app}=48\,\frac{\mu\text{A}}{\text{cm}^2}\), g _l  = 2 \(\frac{\text{mS}}{\text{cm}^2}\), ε _l  = − 60 mV, g _k  = 8 \(\frac{\text{mS}}{\text{cm}^2}\), ε _k  = − 84 mV, g _ca  = 4 \(\frac{\text{mS}}{\text{cm}^2}\), ε _ca  = 120 mV, g = 0.01 \(\frac{\text{mS}}{\text{cm}^2}\), ε _s  = − 84 mV, \({\tilde \sigma}=1.632\). For the auxiliary functions: v _a  = − 1.2 mV, v _b  = 18 mV, v _c  = 12 mV, v _d  = 17.4 mV, \(\varphi=\frac{1}{15}\,\)ms^− 1, v _t  = 48 mV, v _s  = 1 mV, α = 1, β = 2, \(\varphi_s=\frac{1}{5}\,\)ms^− 1.

Reprints and permissions