JournalofMathematicalNeuroscience(2013)3:2 DOI10.1186/2190-8567-3-2 RESEARCH OpenAccess Phase-Amplitude Descriptions of Neural Oscillator Models K.C.A.Wedgwood·K.K.Lin·R.Thul· S.Coombes Received:31July2012/Accepted:18January2013/Publishedonline:24January2013 ©2013K.C.A.Wedgwoodetal.;licenseeSpringer.ThisisanOpenAccessarticledistributedunderthe termsoftheCreativeCommonsAttributionLicense(http://creativecommons.org/licenses/by/2.0),which permitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkis properlycited. Abstract Phaseoscillatorsareacommonstartingpointforthereduceddescription ofmanysingleneuronmodelsthatexhibitastronglyattractinglimitcycle.Theframe- workforanalysingsuchmodelsinresponsetoweakperturbationsisnowparticularly welladvanced,andhasallowedforthedevelopmentofatheoryofweaklyconnected neuralnetworks.However,thestrong-attractionassumptionmaywellnotbethenat- ural one for many neural oscillator models. For example, the popular conductance basedMorris–Lecarmodelisknowntorespondtoperiodicpulsatilestimulationina chaoticfashionthatcannotbeadequatelydescribedwithaphasereduction.Inthispa- per,wegeneralisethephasedescriptionthatallowsonetotracktheevolutionofdis- tancefromthecycleaswellasphaseoncycle.Weuseaclassicaltechniquefromthe theoryofordinarydifferentialequationsthatmakesuseofamovingcoordinatesys- temtoanalyseperiodicorbits.Thesubsequentphase-amplitudedescriptionisshown tobeverywellsuitedtounderstandingtheresponseoftheoscillatortoexternalstim- uli (which are not necessarily weak). We consider a number of examples of neural oscillatormodels,rangingfromplanarthroughtohighdimensionalmodels,toillus- tratetheeffectivenessofthisapproachinprovidinganimprovementoverthestandard phase-reductiontechnique.Asanexplicitapplicationofthisphase-amplitudeframe- work, we considerin some detailthe response of a genericplanar modelwherethe strong-attractionassumptiondoesnothold,andexaminetheresponseofthesystem (cid:2) K.C.A.Wedgwood( )·R.Thul·S.Coombes SchoolofMathematicalSciences,UniversityofNottingham,Nottingham,NG72RD,UK e-mail:[email protected] R.Thul e-mail:[email protected] S.Coombes e-mail:[email protected] K.K.Lin DepartmentofAppliedMathematics,UniversityofArizona,Tucson,AZ,USA e-mail:[email protected] Page2of22 K.C.A.Wedgwoodetal. toperiodicpulsatileforcing.Inaddition,weexplorehowthepresenceofdynamical shearcanleadtoachaoticresponse. Keywords Phase-amplitude·Oscillator·Chaos·Non-weakcoupling ListofAbbreviations ML Morris–Lecar FHN FitzHugh–Nagumo CS Connor–Stevens LE Lyapunovexponent 1 Introduction Oneonlyhastolookattheplethoraofpapersandbooksonthetopicofphaseoscil- lators in mathematical neuroscience to see the enormous impact that this tool from dynamicalsystemstheoryhashadonthewaywethinkaboutdescribingneuronsand neuralnetworks.Muchof thisworkhas itsroots inthetheoryof ordinarydifferen- tialequations(ODEs)andhasbeenpromotedformanyyearsintheworkofWinfree [1],Guckenheimer[2],Holmes[3],Kopell[4],Ermentrout[5]andIzhikevich[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 non-linear dynamical system (as a model of a neuron) exhibits a stablelimitcycleattractorthentrajectoriesnearthecyclecanbeprojectedontothe cycle. Anaturalphasevariableissimplythetimealongthecycle(fromsomearbitrary origin)relativetotheperiodofoscillation.Thenotionofphasecanevenbeextended off the cycle using the concept of isochrons [1]. They provide global information aboutthe‘latentphase’,namelythephasethatwillbeasymptoticallyreturnedtofor atrajectorywithinitialdatawithinthebasinofattractionofanexponentiallystable periodicorbit.Moretechnically,isochronscanbeviewedastheleavesoftheinvariant foliationofthestablemanifoldofaperiodicorbit[8].Inrotatingframecoordinates givenbyphaseandtheleafoftheisochronfoliation,thesystemhasaskew-product structure, i.e. the equationof the phase decouples.However,it is a major challenge tofindtheisochronfoliation,andsinceitreliesontheknowledgeofthelimitcycleit canonlybefoundinspecialcasesornumerically.Therearenowanumberofcom- plementarytechniquesthattacklethislatterchallenge,andinparticularwereferthe readertoworkofGuillamonandHuguet[9](usingLiesymmetries)andOsingaand Moehlis[10](exploitingnumericalcontinuation).MorerecentworkbyMauroyand Mezic[11]isespeciallyappealingasitusesasimpleforwardintegrationalgorithm, as illustrated in Fig. 1 for a Stuart–Landau oscillator. However, it is more common to side-step the need for constructing global isochrons by restricting attention to a smallneighbourhoodofthelimitcycle,wheredynamicscansimplyberecastinthe reducedformθ˙=1,whereθ isthephasearoundacycle.Thisreductiontoaphase descriptiongivesanicesimpledynamicalsystem,albeitonethatcannotdescribeevo- lutionoftrajectoriesinphase-spacethatarefarawayfromthelimitcycle.However, JournalofMathematicalNeuroscience(2013)3:2 Page3of22 Fig.1 Isochronsofa Stuart–Landauoscillatormodel: x˙=λx/2−(λc/2+ω)y− λ(x2+y2)(x−cy)/2, y˙=(λc/2+ω)x+λy/2− λ(x2+y2)(cx+y)/2.The blackcurverepresentsthe periodicorbitofthesystem, whichissimplytheunitcircle forthismodel.Thebluecurves aretheisochronsobtainedusing the(forward)approachof MauroyandMezic[11].The greendotsareanalytically obtainedisochronalpoints[15]. Parametervaluesareλ=2.0, c=1.0,andω=1.0 the phase reduction formalism is useful in quantifying how a system (on or close toacycle)respondstoweakforcing,viatheconstructionoftheinfinitesimalphase responsecurve(iPRC).Foragivenhighdimensionalconductancebasedmodelthis canbesolvedfornumerically,thoughforsomenormalformdescriptionsclosedform solutionsarealsoknown[12].TheiPRCatapointoncycleisequaltothegradient of the (isochronal) phase at that point.This approach forms the basis for construct- ing models of weakly interacting oscillators, where the external forcing is pictured as a function of the phase of a firing neuron. This has led to a great deal of work onphase-lockingandcentralpatterngenerationinneuralcircuitry(see,forexample [13]). Note that the work in [9] goes beyond the notion of iPRC and introduces in- finitesimalphaseresponsesurfaces(allowingevaluationofphaseadvancementeven when the stimulus is off cycle), and see also the work in [14] on non-infinitesimal PRCs. Theassumptionthatphasealoneisenoughtocapturetheessentialsofneuralre- sponse is one made more for mathematical convenience than being physiologically motivated.Indeed,forthepopulartypeIMorris–Lecar(ML)firingmodelwithstan- dardparameters,directnumericalsimulationswithpulsatileforcingshowresponses thatcannotbeexplainedsolelywith aphase model[16]. The failure of a phasede- scription 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 phaseresponsecurveusesonlylinearinformationabouttheisochronsandnon-linear effectswillcomeintoplaythefurtherwemoveawayfromthelimitcycle.Thisprob- lem 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 perfectinformationaboutisochrons,thephasereductionstillassumespersistenceof thelimit-cycleandinstantaneousrelaxationbacktocycle.However,thepresenceof nearby invariant phase-space structures such as (unstable) fixed points and invari- ant 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 Page4of22 K.C.A.Wedgwoodetal. neighbourhood of a cycle where a phase description is expected to hold. Thus, de- veloping a reduced description, which captures some notion of distance from cycle isakeycomponentofanytheoryofforcedlimitcycleoscillators.Thedevelopment of phase-amplitude models that better characterise the response of popular high di- mensional single neuron models is precisely the topic of this paper. Given that it is amajorchallengetoconstructanisochronalfoliationweusenon-isochronalphase- amplitudecoordinatesasapracticalmethodforobtainingamoreaccuratedescription ofneuralsystems.Recently,Medvedev[18]hasusedthisapproachtounderstandin moredetailthesynchronisationoflinearlycoupledstochasticlimitcycleoscillators. InSect.2,weconsiderageneralcoordinatetransformation,whichrecaststhedy- namicsofasystemintermsofphase-amplitudecoordinates.Thisapproachisdirectly takenfromtheclassicaltheoryforanalysingperiodicorbitsofODEs,originallycon- sidered for planar systems in [19], and for general systems in [20]. We advocate it hereasonewaytomovebeyondapurelyphase-centricperspective.Weillustratethe 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 transfor- mations and derive the evolution equations for the forced phase-amplitude system. Thisreducestothestandardphasedescriptionintheappropriatelimit.Importantly, weshowthatthebehaviourofthephase-amplitudesystemismuchmoreabletocap- ture 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 shearinducedchaos[16].FinallyinSect.4,wediscusstherelevanceofthisworkto developingatheoryofnetworkdynamicsthatcanimproveuponthestandardweak couplingapproach. 2 Phase-AmplitudeCoordinates Throughout this paper, we study the dynamics prescribed by the system x˙ =f(x), x ∈Rn, with solutions x =x(t) that satisfy x(0)=x ∈Rn. We will assume that 0 the system admits an attracting hyperbolic periodic orbit (namely one zero Floquet exponentand the others havingnegativereal part), withperiod (cid:5),such that u(t)= u(t +(cid:5)). A phase θ ∈[0,(cid:5)) is naturally defined from θ(u(t))=tmod(cid:5). It has longbeenknowninthedynamicalsystemscommunityhowtoconstructacoordinate system based on this notion of phase as well as a distance from cycle; see [20] for adiscussion.Infact,ErmentroutandKopell[21]madegooduseofthisapproachto derive the phase-interaction function for networks of weakly connected limit-cycle oscillatorsinthelimitofinfinitelyfastattractiontocycle.However,thisassumption isparticularlyextremeandunlikelytoholdforabroadclassofsingleneuronmodels. Thus,itisinterestingtoreturntothefullphase-amplitudedescription.Inessence,the transformationtothesecoordinatesinvolvessettingupamovingorthonormalsystem aroundthelimitcycle.Oneaxisofthissystemischosentobeinthedirectionofthe tangent vector along the orbit, and the remaining are chosen to be orthogonal. We introducethenormalisedtangentvectorξ as (cid:3) (cid:3) (cid:2)(cid:3) (cid:3) ξ(θ)= du (cid:3)(cid:3)du(cid:3)(cid:3). (1) dθ dθ JournalofMathematicalNeuroscience(2013)3:2 Page5of22 Fig.2 Demonstrationofthe movingorthonormalcoordinate systemalonganorbitsegment. Astevolvesfromt0tot1,the coordinatesvarysmoothly.In thisplanarexample,ζ always pointstotheoutsideoftheorbit Theremainingcoordinateaxesareconvenientlygroupedtogetherasthecolumnsof ann×(n−1)matrixζ.Inthiscasewecanwriteanarbitrarypointx as x(θ,ρ)=u(θ)+ζ(θ)ρ. (2) Here,|ρ|representstheEuclideandistancefromthelimitcycle.Acaricature,inR2, ofthecoordinatesystemalonganorbitsegmentisshowninFig.2.Throughtheuse of the variable ρ,we canconsider pointsawayfrom theperiodicorbit. Rather than beingisochronal,linesofconstantθ aresimplystraightlinesthatemanatingfroma pointontheorbitinthedirectionofthenormal.Thetechnicaldetailsofspecifying theorthonormalcoordinatesformingζ arediscussedinAppendixA. Uponprojectingthedynamicsontothemovingorthonormalsystem,weobtainthe dynamicsofthetransformedsystem: θ˙=1+f (θ,ρ), ρ˙=A(θ)ρ+f (θ,ρ), (3) 1 2 where (cid:4) (cid:5) dζ f (θ,ρ)=−hT(θ,ρ) ρ+hT(θ,ρ) f(u+ζρ)−f(u) , (4) 1 dθ (cid:4) (cid:5) dζ f (θ,ρ)=−ζT ρf +ζT f(u+ζρ)−f(u)−Dfζρ , (5) 2 1 dθ (cid:6)(cid:3) (cid:3) (cid:7) (cid:6) (cid:7) h(θ,ρ)= (cid:3)(cid:3)(cid:3)du(cid:3)(cid:3)(cid:3)+ξT dζρ −1ξ, A(θ)=ζT −dζ +Dfζ , (6) dθ dθ dθ 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 f (θ,ρ)→0 as |ρ|→0, f (θ,0)=0 and that ∂f (θ,0)/∂ρ =0. In the 1 2 2 above, f captures the shear present in the system, that is, whether the speed of θ 1 increasesordecreasesdependentonthedistancefromcycle.Aprecisedefinitionfor shearisgivenin[22].Additionally,A(θ)describestheθ-dependentrateofattraction orrepulsionfromcycle. It is pertinent to consider where this coordinate transformation breaks down, thatis,wherethedeterminantoftheJacobianofthetransformationK =det[∂x/∂θ ∂x/∂ρ] vanishes. This never vanishes on-cycle (where ρ =0), but may do so for some|ρ|=k>0.Thissetsanupperboundonhowfarawayfromthelimitcyclewe candescribethesystemusingthesephase-amplitudecoordinates.InFig.3,weplot thecurvealongwhichthetransformationbreaksdownfortheMLmodel.Weobserve that,forsomevaluesofθ,k isrelativelysmaller.Thebreakdownoccurswherelines Page6of22 K.C.A.Wedgwoodetal. Fig.3 Thisfigureshowsthe determinantKofthe phase-amplitudetransformation fortheMLmodel.Colours indicatethevalueofK.Thered contourindicateswhereK=0, andthuswherethecoordinate transformationbreaksdown. Notehowallofthevaluesfor whichthisoccurshaveρ<0. Parametervaluesasin AppendixB.1 ofconstantθ cross,andthusthetransformationceasestobeinvertible,andtheseval- uesofθ correspondtopointsalongwhichtheorbithashighcurvature.Wenotethat thisissueislessproblematicalinhigherdimensionalmodels. Ifwenowconsiderthedrivensystem, x˙=f(x)+εg(x,t), x∈Rn, (7) whereε isnotnecessarilysmall,wemayapplythesametransformationasaboveto obtainthedynamicsin(θ,ρ)coordinates,whereθ ∈[0,(cid:5))andρ∈Rn−1,as (cid:8) (cid:9) θ˙=1+f (θ,ρ)+εhT(θ,ρ)g u(θ)+ζ(θ)ρ,t , (8) 1 (cid:8) (cid:9) ρ˙=A(θ)ρ+f (θ,ρ)+εζTB(θ,ρ)g u(θ)+ζ(θ)ρ,t , (9) 2 with dζ B(θ,ρ)=I − ρhT(θ,ρ), (10) n dθ and I is the n×n identity matrix. Here, h and B describe the effect in terms of θ n andρ thattheperturbationshave.DetailsofthederivationaregiveninAppendixA. For planar models, B =I . To demonstrate the application of the above coordinate 2 transformation,wenowconsidersomepopularsingleneuronmodels. 2.1 A2DConductanceBasedModel TheMLmodelwasoriginallydevelopedtodescribethevoltagedynamicsofbarnacle giant muscle fibre [23], and is now a popular modelling choice in computational neuroscience[7].Itiswrittenasapairofcouplednon-linearODEsoftheform Cv˙=I(t)−gl(v−vl)−gKw(v−vK)−gCam∞(v)(v−vCa), (cid:8) (cid:9) (11) w˙ =φ w∞(v)−w /τw(v). Here,v isthemembranevoltage,whilstw isagatingvariable,describingthefrac- tionofmembraneionchannelsthatareopenatanytime.Thefirstequationexpresses JournalofMathematicalNeuroscience(2013)3:2 Page7of22 Fig.4 TypicaltrajectoryoftheMLmodelofthetransformedsystem.Left:Timeevolutionofθ andρ. Right:Trajectoryplottedinthe(v,w)phaseplane.Weseethatwhenρhasalocalmaximum,theevolution ofθslowsdown,correspondingtowheretrajectoriespassneartothesaddlenode.Parametervaluesasin AppendixB.1 Kirchoff’scurrentlawacrossthecellmembrane,withI(t)representingastimulusin theformofaninjectedcurrent.ThedetailedformofthemodeliscompletedinAp- pendixB.1.TheMLmodelhasaveryrichbifurcationstructure.Roughlyspeaking, byvaryingaconstantcurrentI(t)≡I ,oneobserves,indifferentparameterregions, 0 dynamicalregimescorrespondingto sinks, limitcycles,and Hopf,saddle-nodeand homoclinic bifurcations, as well as combinations of the above. These scenarios are discussedindetailin[7]and[24]. As the ML model is planar, ρ is a scalar, as are the functions A and f . This 1,2 allows us to use the moving coordinate system to clearly visualise parts of phase spacewheretrajectoriesareattractedtowardsthelimitcycle,andpartsinwhichthey move away from it, as illustrated in Fig. 4. The functions f and A, evaluated at 1,2 ρ =−0.1 are shown in Fig. 5. The evolution of θ is mostly constant, however we clearlyobserveportionsofthetrajectorieswherethisisslowed,alongwhichρ˙≈0. In fact, this corresponds to where trajectories pass near to the saddle node, and the dynamicsstall.Thisoccursaroundθ =17,andinFig.5weseethatbothA(θ)and f (θ,ρ)areindeedcloseto0.Thereducedvelocitiesoftrajectoriesherehighlights 1 the importance of considering other phase space structures in forced systems, the details of which are missed in standard phase only models. Forcing in the presence ofsuchstructuresmaygiverisetocomplexandevenchaoticbehaviours,asweshall seeinSect.3. Inthenextexample,weshowhowthesameideasgoacrosstohigherdimensional models. 2.2 A4DConductanceBasedModel TheConnor–Stevens(CS)model[25]isbuiltupontheHodgkin–Huxleyformalism + + andcomprisesafastNa current,adelayedK current,aleakcurrentandatransient + K current, termed the A-current. The full CS model consists of 6 equations: the membranepotential,theoriginalHodgkin–Huxleygatingvariables,andanactivating and inactivating gating variable for the A-current. Using the method of equivalent Page8of22 K.C.A.Wedgwoodetal. Fig.5 f1,f2,andAfortheMLmodel,evaluatedatρ=−0.1.Weclearlyseethedifferenceintheorder ofmagnitudebetweenf1 andf2 forsmallρ.Notethat,althoughtheaverageofAoveroneperiodis negative,itispositiveforanon-trivialintervalofθ.Thiscorrespondstomovementclosetothestable manifoldofthesaddlenode.ParametervaluesasinAppendixB.1 potentials [26], we may reduce the dimensionality of the system, to include only 4 variables.Thereducedsystemis Cv˙=−F(v,u,a,b)+I, u˙=G(v,u), X˙ = X∞(v)−X, X∈{a,b}, (12) τ (v) X where F(v,u,a,b)=gKn4∞(u)(v−vK)+gNah∞(u)m3∞(v)(v−vNa) +gl(v−v)+g a3b(v−v ). (13) l l a a ThedetailsofthereducedCSmodelarecompletedinAppendixB.2.Thesolutions tothereducedCSmodelunderthecoordinatetransformationmaybeseeninFig.6, whilst,inFig.7,weshowhowthissolutionlooksintheoriginalcoordinates.Asfor theMLmodel,θ evolvesapproximatelyconstantlythroughout,thoughthisevolution is sped up close to θ =(cid:5). The trajectories of the vector ρ are more complicated, but note that there is regularity in the pattern exhibited, and that this occurs with approximately the same period as the underlying limit cycle. The damping of the JournalofMathematicalNeuroscience(2013)3:2 Page9of22 Fig.6 SolutionofthetransformedCSsystem.Top:Timeevolutionofθ.Bottom:Timeevolutionofρ coordinates.Upontransformingthesecoordinatesbacktotheoriginalones,wearriveatFig.7.Parameter valuesgiveninAppendixB.2.Inthisparameterregime,themodelexhibitstypeIfiringdynamics amplitudeofoscillationsinρoversuccessiveperiodsrepresentstheoverallattraction tothelimitcycle,whilsttheregularbehaviourofρ representsthespecificrelaxation to cycle as shown in Fig. 7. Additional file 1 shows a movie of the trajectory in (v,u,b) coordinateswiththemovingorthonormalsystemsuperimposed,aswellas thesolutioninρ forcomparison. 3 PulsatileForcingofPhase-AmplitudeOscillators Wenowconsiderasystemwithtime-dependentforcing,givenby(7)with (cid:10)(cid:8) (cid:9) g(x,t)= δ(t−nT),0,...,0 T, (14) n∈Z whereδ istheDiracδ-function.ThisdescribesT-periodickickstothevoltagevari- able. Even such a simple forcing paradigm can give rise to rich dynamics [16]. For theperiodicallykickedMLmodel,shearforcescanleadtochaoticdynamicsasfolds 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 neuralresponse,thismeansthattheneuronisunreliable[27]. Thebehaviourofoscillatorsundersuchperiodicpulsatileforcingisthesubjectof anumberofstudies;see,e.g.[27–30].Ofparticularrelevancehereis[27],inwhich a qualitative reasoning of the mechanisms that bring about shear in such models is Page10of22 K.C.A.Wedgwoodetal. Fig.7 Transformedtrajectoryin(v,u,b)spaceofthephase-amplitudedescriptionofthereducedCS model.Thedottedblacklineisthephaseamplitudesolutiontransformedintheoriginalcoordinates,whilst thecolouredorbitistheunderlyingperiodicorbit,wherethecolourcorrespondstothephasealongthe orbit.Thesolutionofthephase-amplitudedescriptionofthemodelin(θ,ρ)coordinatesisshowninFig.6 supplementedbydirectnumericalsimulationstodetectthepresenceofchaoticsolu- tions.FortheMLmodelinaparameterregionclosetothehomoclinicregime,kicks can cause trajectories to pass near the saddle-node, and folds may occur as a result [16]. We here would like to compare full planar neural models to the simple model, studiedin[27]: (cid:10) θ˙=1+σρ, ρ˙=−λρ+εP(θ) δ(t−nT). (15) n∈Z 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 thatthefunction P is smoothbutnon-constant,trajectorieswillbetakena variable distance from the cycle upon the application of the kick. When kicks are repeated, thisgeometricmechanismcanleadtorepeatedstretchingandfoldingofphasespace. Itisclearthatthelarger σ isin(15),themoreshearispresent,andthemorelikely wearetoobservethefoldingeffect.Inasimilarway,smallervaluesofλmeanthat theshearhaslongertoactupontrajectoriesandagainresultinagreaterlikelihoodof chaos.Finally,toobservechaoticresponse,wemustensurethattheshearforceshave sufficienttimetoact,meaningthatT,thetimebetweenkicksmustnotbetoosmall.