Abstract
Phase oscillators are a common starting point for the reduced description of many single neuron models that exhibit a strongly attracting limit cycle. The framework for analysing such models in response to weak perturbations is now particularly well advanced, and has allowed for the development of a theory of weakly connected neural networks. However, the strongattraction assumption may well not be the natural one for many neural oscillator models. For example, the popular conductance based Morris–Lecar model is known to respond to periodic pulsatile stimulation in a chaotic fashion that cannot be adequately described with a phase reduction. In this paper, we generalise the phase description that allows one to track the evolution of distance from the cycle as well as phase on cycle. We use a classical technique from the theory of ordinary differential equations that makes use of a moving coordinate system to analyse periodic orbits. The subsequent phaseamplitude description is shown to be very well suited to understanding the response of the oscillator to external stimuli (which are not necessarily weak). We consider a number of examples of neural oscillator models, ranging from planar through to high dimensional models, to illustrate the effectiveness of this approach in providing an improvement over the standard phasereduction technique. As an explicit application of this phaseamplitude framework, we consider in some detail the response of a generic planar model where the strongattraction assumption does not hold, and examine the response of the system to periodic pulsatile forcing. In addition, we explore how the presence of dynamical shear can lead to a chaotic response.
Keywords:
Phaseamplitude; Oscillator; Chaos; Nonweak coupling1 Introduction
One only has to look at the plethora of papers and books on the topic of phase oscillators in mathematical neuroscience to see the enormous impact that this tool from dynamical systems theory has had on the way we think about describing neurons and neural networks. Much of this work has its roots in the theory of ordinary differential equations (ODEs) and has been promoted for many years in the work of Winfree [1], Guckenheimer [2], Holmes [3], Kopell [4], Ermentrout [5] and Izhikevich [6] to name but a few. For a recent survey, we refer the reader to the book by Ermentrout and Terman [7]. At heart, the classic phase reduction approach recognises that if a high dimensional nonlinear dynamical system (as a model of a neuron) exhibits a stable limit cycle attractor then trajectories near the cycle can be projected onto the cycle.
A natural phase variable is simply the time along the cycle (from some arbitrary origin)
relative to the period of oscillation. The notion of phase can even be extended off
the cycle using the concept of isochrons [1]. They provide global information about the ‘latent phase’, namely the phase that
will be asymptotically returned to for a trajectory with initial data within the basin
of attraction of an exponentially stable periodic orbit. More technically, isochrons
can be viewed as the leaves of the invariant foliation of the stable manifold of a
periodic orbit [8]. In rotating frame coordinates given by phase and the leaf of the isochron foliation,
the system has a skewproduct structure, i.e. the equation of the phase decouples.
However, it is a major challenge to find the isochron foliation, and since it relies
on the knowledge of the limit cycle it can only be found in special cases or numerically.
There are now a number of complementary techniques that tackle this latter challenge,
and in particular we refer the reader to work of Guillamon and Huguet [9] (using Lie symmetries) and Osinga and Moehlis [10] (exploiting numerical continuation). More recent work by Mauroy and Mezic [11] is especially appealing as it uses a simple forward integration algorithm, as illustrated
in Fig. 1 for a Stuart–Landau oscillator. However, it is more common to sidestep the need
for constructing global isochrons by restricting attention to a small neighbourhood
of the limit cycle, where dynamics can simply be recast in the reduced form
Fig. 1. Isochrons of a Stuart–Landau oscillator model:
The assumption that phase alone is enough to capture the essentials of neural response is one made more for mathematical convenience than being physiologically motivated. Indeed, for the popular type I Morris–Lecar (ML) firing model with standard parameters, direct numerical simulations with pulsatile forcing show responses that cannot be explained solely with a phase model [16]. The failure of a phase description is in itself no surprise and underlies why the community emphasises the use of the word weakly in the phrase “weakly connected neural networks”. Indeed, there are a number of potential pitfalls when applying phase reduction techniques to a system that is not in a weakly forced regime. The typical construction of the phase response curve uses only linear information about the isochrons and nonlinear effects will come into play the further we move away from the limit cycle. This problem can be diminished by taking higher order approximations to the isochrons and using this information in the construction of a higher order PRC [17]. Even using perfect information about isochrons, the phase reduction still assumes persistence of the limitcycle and instantaneous relaxation back to cycle. However, the presence of nearby invariant phasespace structures such as (unstable) fixed points and invariant manifolds may result in trajectories spending long periods of time away from the limit cycle. Moreover, strong forcing will necessarily take one away from the neighbourhood of a cycle where a phase description is expected to hold. Thus, developing a reduced description, which captures some notion of distance from cycle is a key component of any theory of forced limit cycle oscillators. The development of phaseamplitude models that better characterise the response of popular high dimensional single neuron models is precisely the topic of this paper. Given that it is a major challenge to construct an isochronal foliation we use nonisochronal phaseamplitude coordinates as a practical method for obtaining a more accurate description of neural systems. Recently, Medvedev [18] has used this approach to understand in more detail the synchronisation of linearly coupled stochastic limit cycle oscillators.
In Sect. 2, we consider a general coordinate transformation, which recasts the dynamics of a system in terms of phaseamplitude coordinates. This approach is directly taken from the classical theory for analysing periodic orbits of ODEs, originally considered for planar systems in [19], and for general systems in [20]. We advocate it here as one way to move beyond a purely phasecentric perspective. We illustrate the transformation by applying it to a range of popular neuron models. In Sect. 3, we consider how inputs to the neuron are transformed under these coordinate transformations and derive the evolution equations for the forced phaseamplitude system. This reduces to the standard phase description in the appropriate limit. Importantly, we show that the behaviour of the phaseamplitude system is much more able to capture that of the original single neuron model from which it is derived. Focusing on pulsatile forcing, we explore the conditions for neural oscillator models to exhibit shear induced chaos [16]. Finally in Sect. 4, we discuss the relevance of this work to developing a theory of network dynamics that can improve upon the standard weak coupling approach.
2 PhaseAmplitude Coordinates
Throughout this paper, we study the dynamics prescribed by the system
The remaining coordinate axes are conveniently grouped together as the columns of
an
Here,
Fig. 2. Demonstration of the moving orthonormal coordinate system along an orbit segment.
As t evolves from
Upon projecting the dynamics onto the moving orthonormal system, we obtain the dynamics of the transformed system:
where
and Df is the Jacobian of the vector field f evaluated along the periodic orbit u. The derivation of this system may be found in Appendix A. It is straightforward
to show that
It is pertinent to consider where this coordinate transformation breaks down, that is,
where the determinant of the Jacobian of the transformation
Fig. 3. This figure shows the determinant K of the phaseamplitude transformation for the ML model. Colours indicate the value of K. The red contour indicates where
If we now consider the driven system,
where ε is not necessarily small, we may apply the same transformation as above to obtain
the dynamics in
with
and
2.1 A 2D Conductance Based Model
The ML model was originally developed to describe the voltage dynamics of barnacle giant muscle fibre [23], and is now a popular modelling choice in computational neuroscience [7]. It is written as a pair of coupled nonlinear ODEs of the form
Here, v is the membrane voltage, whilst w is a gating variable, describing the fraction of membrane ion channels that are open at any time. The
first equation expresses Kirchoff’s current law across the cell membrane, with
As the ML model is planar, ρ is a scalar, as are the functions A and
Fig. 4. Typical trajectory of the ML model of the transformed system. Left: Time evolution of θ and ρ. Right: Trajectory plotted in the
Fig. 5.
In the next example, we show how the same ideas go across to higher dimensional models.
2.2 A 4D Conductance Based Model
The Connor–Stevens (CS) model [25] is built upon the Hodgkin–Huxley formalism and comprises a fast Na^{+} current, a delayed K^{+} current, a leak current and a transient K^{+} current, termed the Acurrent. The full CS model consists of 6 equations: the membrane potential, the original Hodgkin–Huxley gating variables, and an activating and inactivating gating variable for the Acurrent. Using the method of equivalent potentials [26], we may reduce the dimensionality of the system, to include only 4 variables. The reduced system is
where
The details of the reduced CS model are completed in Appendix B.2. The solutions
to the reduced CS model under the coordinate transformation may be seen in Fig. 6, whilst, in Fig. 7, we show how this solution looks in the original coordinates. As for the ML model,
θ evolves approximately constantly throughout, though this evolution is sped up close
to
Fig. 6. Solution of the transformed CS system. Top: Time evolution of θ. Bottom: Time evolution of ρ coordinates. Upon transforming these coordinates back to the original ones, we arrive at Fig. 7. Parameter values given in Appendix B.2. In this parameter regime, the model exhibits type I firing dynamics
Fig. 7. Transformed trajectory in
3 Pulsatile Forcing of PhaseAmplitude Oscillators
We now consider a system with timedependent forcing, given by (7) with
where δ is the Dirac δfunction. This describes Tperiodic kicks to the voltage variable. Even such a simple forcing paradigm can give rise to rich dynamics [16]. For the periodically kicked ML model, shear forces can lead to chaotic dynamics as folds and horseshoes accumulate under the forcing. This means that the response of the neuron is extremely sensitive to the initial phase when the kicks occur. In terms of neural response, this means that the neuron is unreliable [27].
The behaviour of oscillators under such periodic pulsatile forcing is the subject of a number of studies; see, e.g. [2730]. Of particular relevance here is [27], in which a qualitative reasoning of the mechanisms that bring about shear in such models is supplemented by direct numerical simulations to detect the presence of chaotic solutions. For the ML model in a parameter region close to the homoclinic regime, kicks can cause trajectories to pass near the saddlenode, and folds may occur as a result [16].
We here would like to compare full planar neural models to the simple model, studied in [27]:
This system exhibits dynamical shear, which under certain conditions, can lead to chaotic dynamics. The shear parameter σ dictates how much trajectories are ‘sped up’ or ‘slowed down’ dependent on their distance from the limit cycle, whilst λ is the rate of attraction back to the limit cycle, which is independent of θ. Supposing that the function P is smooth but nonconstant, trajectories will be taken a variable distance from the cycle upon the application of the kick. When kicks are repeated, this geometric mechanism can lead to repeated stretching and folding of phase space. It is clear that the larger σ is in (15), the more shear is present, and the more likely we are to observe the folding effect. In a similar way, smaller values of λ mean that the shear has longer to act upon trajectories and again result in a greater likelihood of chaos. Finally, to observe chaotic response, we must ensure that the shear forces have sufficient time to act, meaning that T, the time between kicks must not be too small.
This stretching and folding action can clearly lead to the formation of Smale horseshoes, which are well known to lead to a type of chaotic behaviour. However, horseshoes may coexist with sinks, meaning the resulting chaotic dynamics would be transient. Wang and Young proved that under appropriate conditions, there is a set of T of positive Lebesgue measure for which the system experiences a stronger form of sustained, chaotic behaviour, characterised by the existence of a positive Lyapunov exponent for almost all initial conditions and the existence of a ‘strange attractor’; see, e.g. [2830].
By comparing with the phaseamplitude dynamics described by Eqs. (8)–(9), we see that
the model of shear considered in (15) is a proxy for a more general system, with
To gain a deeper insight into the phenomenon of shear induced chaos, it is pertinent
to study the isochrons of the limit cycle for the linear model (15), where the isochrons
are simply lines with slope
Fig. 8. Stretchandfold action of a kick followed by relaxation in the presence of shear.
The thin black lines are the isochrons of the system, which in the case of the linear model (15), are
simply straight lines with slope
This simple model with a harmonic form for
The formation of chaos in the ML model is shown in Fig. 9. Here, we plot the response to periodic pulsatile forcing, given by (14), in the
Fig. 9. Shear induced folding in the ML model, with parameters as in Fig. 5. The red curve in all panels represents the limit cycle of the unperturbed system, whilst the green dotted line represents the stable manifold of the saddle node, indicated by the orange marker. We begin distributing points along the limit cycle and then apply an instantaneous
kick taking
Compared to the phenomenological system (15), models written in phaseamplitude coordinates
as (8)–(9) have two main differences. The intrinsic dynamics (without kicks) are nonlinear,
and the kick terms appear in both equations for
It is important to emphasise that
where
In the interval
Rescaling time as
with initial conditions
Note that this has been constructed using a matched asymptotic expansion, using (17),
and is valid in the limit
Although this explicit map is convenient for numerical simulations, we prefer to work with the full stroboscopic map (22)–(23), which is particularly useful for comparing and contrasting the behaviour of different planar single neuron models with arbitrary kick strength. As an indication of the presence of chaos in the dynamics resulting from this system, we evaluate the largest Lyapunov exponent of the map (22)–(23) by numerically evolving a tangent vector and computing its rate of growth (or contraction); see e.g. [32] for details.
In Fig. 10, we compare the functions
Fig. 10. The blue curves show the change in θ under the action of a pulsatile kick in v, whilst the red dashed curves show the change in ρ under the same kick. The top plot is for the FHN model, whilst the bottom plot is for the ML model. We evaluate the effect of the kicks at
We plot in the top row of Fig. 11 the pair
Fig. 11. Panel (a) shows successive iterates of θ in system (22)–(23) with functions
Now that we know the relative contribution of kicks in v to kicks in
After transients, we observe a
4 Discussion
In this paper, we have used the notion of a moving orthonormal coordinate system around
a limit cycle to study dynamics in a neighbourhood around it. This phaseamplitude
coordinate system can be constructed for any given ODE system supporting a limit cycle.
A clear advantage of the transformed description over the original one is that it
allows us to gain insight into the effect of time dependent perturbations, using the
notion of shear, as we have illustrated by performing case studies of popular neural
models, in two and higher dimensions. Whilst this coordinate transformation does not
result in any reduction in dimensionality in the system, as is the case with classical
phase reduction techniques, it opens up avenues for moving away from the weak coupling
limit, where
One natural extension of this work is to move beyond the theory of weakly coupled oscillators to develop a framework for describing neural systems as networks of phaseamplitude units. This has previously been considered for the case of weakly coupled weakly dissipative networks of nonlinear planar oscillators (modelled by small dissipative perturbations of a Hamiltonian oscillator) [3335]. It would be interesting to develop these ideas and obtain network descriptions of the following type:
with an appropriate identification of the interaction functions
List of Abbreviations
ML Morris–Lecar
FHN FitzHugh–Nagumo
CS Connor–Stevens
LE Lyapunov exponent
Appendix A: Derivation of the Transformed Dynamical System
Starting from
we make the transformation
We proceed by projecting (29) onto
where
Thus,
where
and
Upon projecting both sides of (29) onto
whilst the righthand side becomes
since
where
It may be easily seen that
In order to evaluate the functions
where the index i denotes the column entry of ζ and
and
where j denotes the row index, we have
By the quotient rule for vectors we find that
and that
Overall, we have that
Appendix B: Gallery of Models
B.1 Morris–Lecar
The ML equations describe the interaction of membrane voltage with just two ionic
currents:
The function
B.2 Reduced Connor–Stevens Model
For the reduced CS model, we start with the full Hodgkin–Huxley model, with m, n, h as gating variables and use the method of equivalent potentials as treated in [26], giving rise to the following form for the function g:
where
Parameter values are
B.3 FitzHugh–Nagumo Model
The FHN model is a phenomenological model of spike generation, comprising of 2 variables. The first represents the membrane potential and includes a cubic nonlinearity, whilst the second variable is a gating variable, similar to w in the ML model, which may be thought of as a recovery variable. The system is
where we use the following parameter values:
Electronic Supplementary Material
Additional file 1. A movie showing the moving orthonormal system for the Connor–Stevens model. The top
panel shows a projection of the moving orthonormal system from the full
Format: MOV Size: 648KB Download file  Watch movie
Additional file 2. A movie showing the stretchandfold action brought about by shear forces. In this
movie, we show both the shear forces and the rate of attraction back to cycle are
linear. The limit cycle is first unravelled so that it may be represented by a straight
line. We choose
Format: MOV Size: 769KB Download file  Watch movie
Additional file 3. A movie showing the accumulation of folds in the kicked ML model. The thin black line
represents the underlying periodic orbit of the system
Format: MOV Size: 2.1MB Download file  Watch movie
Competing Interests
The authors confirm that they have no competing interests of which they are aware.
Authors’ Contributions
KCAW, KKL, RT and SC contributed equally. All authors read and approved the final manuscript.
References

Winfree A: The Geometry of Biological Time. 2nd edition. Springer, Berlin; 2001.

Guckenheimer J: Isochrons and phaseless sets.
J Math Biol 1975, 1:259273. Publisher Full Text

Cohen AH, Rand RH, Holmes PJ: Systems of coupled oscillators as models of central pattern generators. In Neural Control of Rhythmic Movements in Vertebrates. Wiley, New York; 1988.

Kopell N, Ermentrout GB: Symmetry and phaselocking in chains of weakly coupled oscillators.
Commun Pure Appl Math 1986, 39:623660. Publisher Full Text

Ermentrout GB: n:m phaselocking of weakly coupled oscillators.
J Math Biol 1981, 12:327342. Publisher Full Text

Izhikevich EM: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press, Cambridge; 2007.

Ermentrout GB, Terman DH: Mathematical Foundations of Neuroscience. Springer, Berlin; 2010.

Josic K, SheaBrown ET, Moehlis J: Isochron.
Scholarpedia 2006., 1(8)
Article ID 1361

Guillamon A, Huguet G: A computational and geometric approach to phase resetting curves and surfaces.
SIAM J Appl Dyn Syst 2009, 8(3):10051042. Publisher Full Text

Osinga HM, Moehlis J: A continuation method for computing global isochrons.
SIAM J Appl Dyn Syst 2010, 9(4):12011228. Publisher Full Text

Mauroy A, Mezic I: On the use of Fourier averages to compute the global isochrons of (quasi)periodic dynamics.
Chaos 2012., 22(3)
Article ID 033112

Brown E, Moehlis J, Holmes P: On the phase reduction and response dynamics of neural oscillator populations.
Neural Comput 2004, 16:673715. PubMed Abstract  Publisher Full Text

Hoppensteadt FC, Izhikevich EM: Weakly Connected Neural Networks. Springer, Berlin; 1997.

Achuthan S, Canavier CC: Phaseresetting curves determine synchronization, phase locking, and clustering in networks of neural oscillators.
J Neurosci 2009, 29(16):52185233. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yoshimura K: Phase reduction of stochastic limitcycle oscillators. In Reviews of Nonlinear Dynamics and Complexity. Volume 3. Wiley, New York; 2010:5990.

Lin KK, Wedgwood KCA, Coombes S, Young LS: Limitations of perturbative techniques in the analysis of rhythms and oscillations.
J Math Biol 2013, 66:139161. PubMed Abstract  Publisher Full Text

Demir A, Suvak O: Quadratic approximations for the isochrons of oscillators: a general theory, advanced numerical methods and accurate phase computations.
IEEE Trans ComputAided Des Integr Circuits Syst 2010, 29:12151228.

Medvedev GS: Synchronization of coupled stochastic limit cycle oscillators.
Phys Lett A 2010, 374:17121720. Publisher Full Text

Diliberto SP: On systems of ordinary differential equations. In Contributions to the Theory of Nonlinear Oscillations. Princeton University Press, Princeton; 1950:138. [Annals of Mathematical Studies 20]

Hale JK: Ordinary Differential Equations. Wiley, New York; 1969.

Ermentrout GB, Kopell N: Oscillator death in systems of coupled neural oscillators.
SIAM J Appl Math 1990, 50:125146. Publisher Full Text

Ott W, Stenlund M: From limit cycles to strange attractors.
Commun Math Phys 2010, 296:215249. Publisher Full Text

Morris C, Lecar H: Voltage oscillations in the barnacle giant muscle fiber.
Biophys J 1981, 35:193213. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rinzel J, Ermentrout GB: Analysis of neural excitability and oscillations. In Methods in Neuronal Modeling. 1st edition. MIT Press, Cambridge; 1989:135169.

Connor JA, Stevens CF: Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone soma.
J Physiol 1971, 213:3153. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kepler TB, Abbott LF, Marder E: Reduction of conductancebased neuron models.
Biol Cybern 1992, 66:381387. PubMed Abstract  Publisher Full Text

Lin KK, Young LS: Shearinduced chaos.
Nonlinearity 2008, 21(5):899922. Publisher Full Text

Wang Q, Young LS: Strange attractors with one direction of instability.
Commun Math Phys 2001, 218:197. Publisher Full Text

Wang Q, Young LS: From invariant curves to strange attractors.
Commun Math Phys 2002, 225:275304. Publisher Full Text

Wang Q: Strange attractors in periodicallykicked limit cycles and Hopf bifurcations.

Catllá AJ, Schaeffer DG, Witelski TP, Monson EE, Lin AL: On spiking models for synaptic activity and impulsive differential equations.
SIAM Rev 2008, 50:553569. Publisher Full Text

Christiansen F, Rugh F: Computing Lyapunov spectra with continuous Gram–Schmidt orthonormalization.
Nonlinearity 1997, 10:10631072. Publisher Full Text

Ashwin P: Weak coupling of strongly nonlinear, weakly dissipative identical oscillators.

Ashwin P, Dangelmayr G: Isochronicityinduced bifurcations in systems of weakly dissipative coupled oscillators.
Dyn Stab Syst 2000, 15(3):263286. Publisher Full Text

Ashwin P, Dangelmayr G: Reduced dynamics and symmetric solutions for globally coupled weakly dissipative oscillators.
Dyn Syst 2005, 20(3):333367. Publisher Full Text

Han SK, Kurrer C, Kuramoto Y: Dephasing and bursting in coupled neural oscillators.
Phys Rev Lett 1995, 75:31903193. PubMed Abstract  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(3):11011129. Publisher Full Text