Abstract
We consider a coupled, heterogeneous population of relaxation oscillators used to model rhythmic oscillations in the preBötzinger complex. By choosing specific values of the parameter used to describe the heterogeneity, sampled from the probability distribution of the values of that parameter, we show how the effects of heterogeneity can be studied in a computationally efficient manner. When more than one parameter is heterogeneous, full or sparse tensor product grids are used to select appropriate parameter values. The method allows us to effectively reduce the dimensionality of the model, and it provides a means for systematically investigating the effects of heterogeneity in coupled systems, linking ideas from uncertainty quantification to those for the study of network dynamics.
Keywords:
heterogeneity; neural oscillators; preBötzinger complex; model reduction; bifurcation; computation1 Introduction
Networks of coupled oscillators have been studied for a number of years [17]. One motivation for these studies is that many neurons, when isolated (and possibly injected with a constant current), either periodically fire action potentials [8,9] or periodically move between quiescence and repetitive firing (the alternation being referred to as bursting [10,11]). In either case, the isolated neuron can be thought of as an oscillator. Neurons are typically coupled with many others via either gap junctions [12] or chemical synapses [1315]; hence, a group of neurons can be thought of as a network of coupled oscillators.
As an idealisation, one might consider identical oscillators; in which case, the symmetry of the network will often determine its possible dynamics [16,17]. However, natural systems are never ideal, and thus, it is more realistic to consider heterogeneous networks. Also, there is evidence in a number of contexts that heterogeneity within a population of neurons can be beneficial. Examples include calcium wave propagation [18], the synchronisation of coupled excitable units to an external drive [19,20], and the example we study here: respiratory rhythm generation [13,21].
One simple way to incorporate heterogeneity in a network of coupled oscillators is to select one parameter which affects the individual dynamics of each oscillator and assign a different value to this parameter for each oscillator [3,15,22,23]. Doing this raises natural questions such as from which distribution should these parameter values be chosen, and what effect does this heterogeneity have on the dynamics of the network?
Furthermore, if we want to answer these questions in the most computationally efficient way, we need a procedure for selecting a (somehow) optimal representative set of parameter values from this distribution. In this paper, we will address some of these issues.
In particular, we will show how  given the distribution(s) of the parameter(s) describing the heterogeneity  the representative set of parameter values can be chosen so as to accurately incorporate the effects of the heterogeneity without having to fully simulate the entire large network of oscillators.
We investigate one particular network of coupled relaxation oscillators, derived from a model of the preBötzinger complex [13,14,24], and show how the heterogeneity in one parameter affects its dynamics. We also show how heterogeneity in more than one parameter can be incorporated using either full or sparse tensor product grids in parameter space.
Our approach thus creates a bridge between computational techniques developed in the field of uncertainty quantification [25,26] involving collocation and sparse grids on the one hand, and network dynamics on the other. It also helps us build accurate, reduced computational models of large coupled neuron populations.
One restriction of our method is that it applies only to states where all oscillators are synchronised (in the sense of having the same period) or at a fixed point. Synchronisation of this form typically occurs when the strength of coupling between oscillators is strong enough to overcome the tendency of nonidentical oscillators to desynchronise due to their disparate frequencies [2,3,27] and is often the behaviour of interest [6,13,14,23].
We present the model in Section 2 and show how to efficiently include parameter heterogeneity in Section 3. In Section 4, we explore how varying heterogeneity modifies bifurcations and varies the period of the collective oscillation. Sections 5 and 6 show how to deal with two and more, respectively, heterogeneous parameters. We conclude in Section 7.
2 The model
Our illustrative example is a network of model neurons thought to describe at some level the dynamics of the preBötzinger complex, governed by the following equations:
for
as considered in the work of Rubin and Terman [14]. Here,
The functions
Note that the synaptic coupling is excitatory. These parameters are the same as that
used in the work of Rubin and Terman [14] except that they [14] used
If the values of the applied current
Fig. 1. Solutions of Equations 1 and 2. These are the solutions when the
3 Managing heterogeneity
3.1 The continuum limit
The key observation behind our approach can be seen in Figure 2, where we plot the
Fig. 2. Solutions of Equations 1 and 2 at one instant in time.
where
The results for
3.2 Stochastic Galerkin
One approach to studying Equations 9 and 11, motivated by techniques developed in
the context of uncertainty quantification [25,26], is to expand the functions
where
3.3 Stochastic collocation
An alternative, motivated by the stochastic collocation method [25], is to simply discretise in the μ direction, obtaining N different
It is important to realize that the number (N) of neurons simulated in this approach may well be much smaller than the number of neurons in the ‘true’ system, considered to be in the thousands. Notice also that these neurons are ‘mathematically’ coupled to one another via the discretisation of the integral (Equation 11), which is an approximation of the continuum limit.
Using the values of
3.3.1 Period calculation
Firstly, we consider the period of the collective oscillations seen in Figure 1. The analogue of finite differences, or the method of lines, is to uniformly discretise
the interval
which, after defining
Fig. 3. Error in the calculated period of the synchronised oscillators. Error in the calculated
period of the synchronised oscillators as a function of the number of neurons simulated
(N) for the midpoint rule (red stars) and the Gaussian quadrature (blue circles). Also shown (dashed) is a line corresponding to error approximately
However, by choosing nonuniformly spaced values of
then the GaussLegendre quadrature rule [31] is
Convergence of the error in the period with N is shown in Figure 3 (blue circles), where we see the very rapid convergence expected from a spectral
method. For
3.3.2 Hopf bifurcations
By decreasing or increasing
Fig. 4. Error in calculation error of the value of
The expected behaviour (very rapid convergence for Gaussian quadrature and the error
scaling as
Fig. 5. Error in calculation error of the value of
3.4 Summary
In this section, we have shown that a judicious choice of the values of the heterogeneous
parameter, combined with a scheme for the Gaussian quadrature, allows us to calculate
quantities of interest (such as the period of oscillation and the parameter value
at which a Hopf bifurcation occurs) much more parsimoniously than a naive implementation
of uniformly spaced
Having demonstrated this, we now fix
In the next section, we investigate the effects of varying
4 The effects of heterogeneity
4.1 A single neuron
In order to investigate the effects of heterogeneity, we first examine a single uncoupled
neuron (i.e.
Fig. 6. The bifurcation behaviour, V as functions of
The behaviour in the left panel of Figure 6 can also be understood by looking at the
Fig. 7. The phase plane for a single uncoupled neuron. hnullcline (dashed, on which
4.2 A coupled population of neurons
Now, consider a coupled heterogeneous population with
Fig. 8. The bifurcation behaviour of a heterogeneous population. Top: mean voltage at a fixed point (solid stable, dashed unstable), mean ± one standard deviation (dotted), and the maximum and minimum of the mean of V over one period of oscillation (circles), as a function of
To better understand these results, we can follow the Hopf bifurcations as two parameters are varied. Figure 9 (top) shows the two curves of Hopf bifurcations in the
Fig. 9. Hopf bifurcation curves and period of the stable periodic orbit for three different
values of
Now, consider varying
Fig. 10. Hopf bifurcation curves and period of the stable periodic orbit for three different
values of
5 Two heterogeneous parameters
Now, consider the case where both
The values of
We initially choose
i.e.λ is normally distributed. Then, as mentioned, the continuum variables V and h are written in the form
Keeping the GaussLegendre rule in the μ direction, this gives
The simplest approach to this integral is the Monte Carlo method [30], where we simply randomly choose M values of λ from the unit normal distribution and calculate an approximation to the integral as the following:
Here, the weights in the λ direction are all equal to
i.e.Q is the cumulative density function for λ, and then for a general function f, the integral
can be written as
Thus, we define
for
where
(The first few probabilists’  as opposed to physicists’  Hermite polynomials are
Fig. 11. Examples of the heterogeneity grid values of
An example of the
The result of using these three different methods to allocate the
Fig. 12. Error in the calculated period using three different methods (see text) for
We see that as expected, the GaussHermite quadrature performs the best, with the
error saturating between
6 Sparse grids
The process described above can obviously be generalised to more than two randomly,
but independently, distributed parameters. The distribution of each parameter determines
the type of quadrature which should be used in that direction, and the parameter values
and weights are formed from tensor products of the underlying onedimensional rules.
However, the curse of dimensionality will restrict how many random parameters can
be accurately sampled. If we use N points in each of D random dimensions, the number of neurons we need to simulate is
One way around this problem is to use sparse grids [33,34], as introduced by Smolyak [35]. The basic idea is to use sparse tensor products, chosen in such a way as to have similar accuracy to the corresponding full tensor product, but with fewer grid points, and thus (in our case) fewer neurons to simulate. A general theory exists [33,34], but to illustrate the idea, suppose we have two uncorrelated random parameters, each is distributed uniformly between −1 and 1. A full tensor product for the GaussLegendre quadrature using 11 points in each direction is shown in Figure 13.
Fig. 13. Full tensor product using 11 points in each direction (121 points in total). The points
are the roots of
To form a twodimensional sparse grid using the GaussLegendre quadrature, we first write the onedimensional integration rule for integrating a function f as
where
i.e.
where
The grid for this rule is shown in Figure 14 (top), along with grids corresponding to several of its components.^{a} Figure 14 (bottom) shows the grid for rule
Fig. 14. Grids for rules
Rules such as these can be constructed for an arbitrary number of spatial dimensions,
using a variety of quadrature rules (and possibly different rules in different dimensions).
Their advantage becomes apparent as the dimension of the space to be integrated over
(or in our case, the number of heterogeneous parameters) is increased. To illustrate
this, we consider as an example the model, Equations 1 and 2 with
Fig. 15. Error in calculation of period. This happens when four distinct parameters are simultaneously heterogeneous (independently of one another) for both full and sparse grids. See text for details. N is the number of neurons simulated.
We see that for fixed N, the sparse grid calculation is approximately two orders or magnitude more accurate than the full grid  implying, in turn, that the way we select the reduced number of neurons we retain to simulate the full system is critical. This relative advantage is expected to increase as the number of distributed parameters increases. As an example of the growth in the number of grid points, a level 6 calculation in 10 dimensions uses fewer than one million points, and the resulting system can be easily simulated on a desktop PC. (Note that the grid points and weights are calculated before the numerical integration starts, so the computational cost in producing data like that shown in Figure 15 is almost entirely due to numerical integration of the ODEs, which is proportional to the number of grid points, i.e. neurons, used.)
7 Discussion
In this paper, we have presented and demonstrated the use of a computationally efficient method for systematically investigating the effects of heterogeneity in the parameters of a coupled network of neural oscillators. The method constitutes a model reduction approach: By only considering oscillators with parameter values given by roots of families of orthogonal polynomials (Legendre, Hermite, …), we can use the Gaussian quadrature to accurately evaluate the term coupling the oscillators, which can be thought of as the discretisation of an integral over the heterogeneous dimension(s).
Effectively, we are simulating the behaviour of an infinite number of oscillators by only simulating a small number of judiciously selected ones, modifying appropriately the way they are coupled. When the oscillators are synchronised, or at a fixed point, the function to be integrated is a smooth function of the heterogeneous parameter(s), and thus, convergence is very rapid. The technique is general (although subject to the restriction immediately above) and can be used when there is more than one heterogeneous parameter, via full or sparse tensor products in parameter space. For a given level of accuracy, we are simulating far fewer neurons than might naively be expected. The emphasis here has been on computational efficiency rather than a detailed investigation of parameter dependence.
The model we considered involved coupling only through the mean of a function, s, of the variable
The primary study with which we should compare our results is that of Rubin and Terman
[14]. They considered essentially the same model as Equations 1 and 2 but with heterogeneity
only in the
Our approach can be thought of as a particular parametrisation of this snake, which takes into account the probability density of the heterogeneity parameter(s); we also showed a systematic way of extending this onedimensional snake to two and higher dimensions. Another paper which uses some of the same ideas as presented here is that of Laing and Kevrekidis [3]. There, the authors considered a finite network of coupled oscillators and used a polynomial chaos expansion of the same form as Equation 12. However, instead of integrating the equations for the polynomial chaos coefficients directly, they used projective integration [39] to do so, in an ‘equationfree’ approach [40] in which the equations satisfied by the polynomial chaos coefficients are never actually derived. They also chose the heterogeneous parameter values randomly from a prescribed distribution and averaged over realisations of this process in order to obtain ‘typical’ results. Similar ideas had been explored earlier by Moon et al.[27], who considered a heterogeneous network of phase oscillators.
Assisi et al.[22] considered a heterogeneous network of coupled neural oscillators, deriving equations of similar functional form to Equations 9 and 11. Their approach was to expand the variables in a way similar to Equation 12 but using a small number of arbitrarily chosen ‘modes’ rather than orthogonal polynomials. Their choice of modes, along with the fact that their neural model consisted of ODEs with polynomial right hand sides, allowed them to analytically derive the ODEs satisfied by the coefficients of the modes. This approach allowed them to qualitatively reproduce some of the behaviour of the network such as the formation of two clusters of oscillators. However, in the general case modes should be chosen as orthogonal polynomials, the specific forms of which are determined by the distribution of the heterogeneous parameter(s) [25,26].
The network we considered was alltoall coupled, and the techniques presented should be applicable to other similar systems. The only requirement is that the relationship between the heterogeneity parameter(s) and the state of the system (possibly after transients) be smooth (or possibly piecewise smooth). An interesting extension is the case when the network under consideration is not alltoall. Then, the effects of degree distribution may affect the dynamics of individual oscillators [38,41,42], and if we have a way of parameterising this type of heterogeneity, it might be possible to apply the ideas presented here to such networks. Degree distribution is a discrete variable, and corresponding families of orthogonal polynomials exist for a variety of discrete random variables [25,26].
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
CRL, YZ and BS performed calculations relating to the results presented. CRL and IGK wrote the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The works of CRL and BS were supported by the Marsden Fund Council from government funding, administered by the Royal Society of New Zealand. The works of IGK and YZ were supported by the AFOSR and the US DOE (DESC0005176 and DESC00029097).
End notes

^{a}These sparse grids were computed using software from http://people.sc.fsu.edu/~jburkardt/ webcite.
References

Baesens C, Guckenheimer J, Kim S, MacKay RS: Three coupled oscillators: modelocking, global bifurcations and toroidal chaos.
Physica D, Nonlinear Phenom 1991, 49(3):387475. Publisher Full Text

Matthews PC, Mirollo RE, Strogatz SH: Dynamics of a large system of coupled nonlinear oscillators.
Physica D, Nonlinear Phenom 1991, 52(23):293331. Publisher Full Text

Laing CR, Kevrekidis IG: Periodicallyforced finite networks of heterogeneous globallycoupled oscillators: a lowdimensional approach.
Physica D, Nonlinear Phenom 2008, 237(2):207215. Publisher Full Text

Martens EA, Laing CR, Strogatz SH: Solvable model of spiral wave chimeras.

Ermentrout B, Pascal M, Gutkin B: The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators.
Neural Comput 2001, 13(6):12851310. PubMed Abstract  Publisher Full Text

Pikovsky A, Rosenblum M, Kurths J: Synchronization. Cambridge University Press, Cambridge; 2001.

Acebrón JA, Bonilla LL, Pérez Vicente CJ, Ritort F, Spigler R: The Kuramoto model: a simple paradigm for synchronization phenomena.
Rev Mod Phys 2005, 77:137185. Publisher Full Text

Hassard B: Bifurcation of periodic solutions of the HodgkinHuxley model for the squid giant axon.
J Theor Biol 1978, 71(3):401420. PubMed Abstract  Publisher Full Text

Ermentrout GB, Terman D: Mathematical Foundations of Neuroscience. Springer, Heidelberg; 2010.

Coombes S, Bressloff PC: Bursting: The Genesis of Rhythm in the Nervous System. World Scientific, Singapore; 2005.

Izhikevich EM: Neural excitability, spiking and bursting.
Int J Bifurc Chaos 2000, 10(6):11711266. Publisher Full Text

Coombes S: Neuronal networks with gap junctions: a study of piecewise linear planar neuron models.
SIAM J Appl Dyn Syst 2008, 7:1101. Publisher Full Text

Butera RJ, Rinzel J, Smith JC: Models of respiratory rhythm generation in the preBötzinger complex. II. Populations of coupled pacemaker neurons.
J Neurophysiol 1999, 82:398. PubMed Abstract  Publisher Full Text

Rubin J, Terman D: Synchronized activity and loss of synchrony among heterogeneous conditional oscillators.
SIAM J Appl Dyn Syst 2002, 1:146174. Publisher Full Text

Golomb D, Rinzel J: Dynamics of globally coupled inhibitory neurons with heterogeneity.
Phys Rev E 1993, 48(6):48104814. Publisher Full Text

Golubitsky M, Stewart I, Buono PL, Collins JJ: Symmetry in locomotor central pattern generators and animal gaits.
Nature 1999, 401(6754):693695. PubMed Abstract  Publisher Full Text

Ashwin P, Swift J: The dynamics of n weakly coupled identical oscillators.
J Nonlinear Sci 1992, 2:69108. Publisher Full Text

Gosak M: Cellular diversity promotes intercellular Ca^{2+} wave propagation.
Biophys Chem 2009, 139:5356. PubMed Abstract  Publisher Full Text

Pérez T, Mirasso CR, Toral R, Gunton JD: The constructive role of diversity in the global response of coupled neuron systems.
Philos Trans R Soc A, Math Phys Eng Sci 2010, 368(1933):56195632. Publisher Full Text

Tessone CJ, Mirasso CR, Toral R, Gunton JD: Diversityinduced resonance.

Purvis LK, Smith JC, Koizumi H, Butera RJ: Intrinsic bursters increase the robustness of rhythm generation in an excitatory network.
J Neurophysiol 2007, 97(2):15151526. PubMed Abstract  Publisher Full Text

Assisi CG, Jirsa VK, Kelso JAS: Synchrony and clustering in heterogeneous networks with global coupling and parameter dispersion.

White JA, Chow CC, Ritt J, SotoTreviño C, Kopell N: Synchronization and oscillatory dynamics in heterogeneous, mutually inhibited neurons.
J Comput Neurosci 1998, 5:516. PubMed Abstract  Publisher Full Text

Butera RJ, Rinzel J, Smith JC: Models of respiratory rhythm generation in the preBötzinger complex. I. Bursting pacemaker neurons.
J Neurophysiol 1999, 82:382. PubMed Abstract  Publisher Full Text

Xiu D: Fast numerical methods for stochastic computations: a review.

Xiu D, Karniadakis GE: Modeling uncertainty in flow simulations via generalized polynomial chaos.
J Comput Phys 2003, 187:137167. Publisher Full Text

Moon SJ, Ghanem R, Kevrekidis IG: Coarse graining the dynamics of coupled oscillators.

Rubin JE: Bursting induced by excitatory synaptic coupling in nonidentical conditional relaxation oscillators or squarewave bursters.

Dunmyre JR, Rubin JE: Optimal intrinsic dynamics for bursting in a threecell network.
SIAM J Appl Dyn Syst 2010, 9:154187. Publisher Full Text

Quarteroni A, Sacco R, Saleri F: Numerical Mathematics. Springer, Heidelberg; 2007.

Trefethen LN: Is Gauss quadrature better than ClenshawCurtis?
SIAM Rev 2008, 50:67. Publisher Full Text

Moehlis J: Canards in a surface oxidation reaction.
J Nonlinear Sci 2002, 12(4):319345. Publisher Full Text

Gerstner T, Griebel M: Numerical integration using sparse grids.
Numer Algorithms 1998, 18(3):209232. Publisher Full Text

Barthelmann V, Novak E, Ritter K: High dimensional polynomial interpolation on sparse grids.
Adv Comput Math 2000, 12:273288. Publisher Full Text

Smolyak SA: Quadrature and interpolation formulas for tensor products of certain classes of functions.

Somers D, Kopell N: Waves and synchrony in networks of oscillators of relaxation and nonrelaxation type.
Physica D, Nonlinear Phenom 1995, 89(12):169183. Publisher Full Text

Abrams DM, Strogatz SH: Chimera states in a ring of nonlocally coupled oscillators.
Int J Bifurc Chaos 2006, 16:2137. Publisher Full Text

Ko TW, Ermentrout GB: Partially locked states in coupled oscillators due to inhomogeneous coupling.

Kevrekidis IG, Gear CW, Hyman JM, Kevrekidis PG, Runborg O, Theodoropoulos C: Equationfree, coarsegrained multiscale computation: enabling macroscopic simulators to perform systemlevel analysis.

Xiu D, Kevrekidis IG, Ghanem R: An equationfree, multiscale approach to uncertainty quantification.
Comput Sci Eng 2005, 7(3):1623. Publisher Full Text

Rajendran K, Kevrekidis IG: Coarse graining the dynamics of heterogeneous oscillators in networks with spectral gaps.

Tsoumanis AC, Siettos CI, Bafas GV, Kevrekidis IG: Computations in social networks: from agentbased modeling to coarsegrained stability and bifurcation analysis.
Int J Bifurc Chaos 2010, 20(11):36733688. Publisher Full Text