On complex dynamics in a Purkinje and a ventricular cardiac cell model

ContentslistsavailableatScienceDirect

Commun Nonlinear Sci Numer Simulat

journalhomepage:www.elsevier.com/locate/cnsns

Research paper

On complex dynamics in a Purkinje and a ventricular cardiac cell model

André H. Erhardt

^a,1,∗

, Susanne Solem

^b,1

aDepartment of Mathematics, University of Oslo, P.O. Box 1053 Blindern, Oslo 0316, Norway

bDepartment of Mathematical Sciences, Norwegian University of Science and Technology, Trondheim 7491, Norway

a r t i c l e i n f o

Article history:

Received 5 May 2020 Revised 12 August 2020 Accepted 27 August 2020 Available online 2 September 2020 2010 MSC:

37G15 37N25 35Q92 65P30 92B05 Keywords:

Nonlinear dynamics Cardiac cells

Reaction–diffusion system Andronov–Hopf and period doubling bifurcation

Deterministic chaos

a b s t r a c t

Cardiacmusclecellscan exhibitcomplexpatterns includingirregular behavioursuchas chaos or (chaotic) early afterdepolarisations (EADs), which can lead to sudden cardiac death.Suitablemathematicalmodelsand theiranalysishelptopredicttheoccurrenceof suchphenomenaandtodecodetheirmechanisms.Thefocusofthispaperistheinves- tigationofdynamicsofcardiacmusclecellsdescribedbysystemsofordinarydifferential equations.Thisisgenericallyperformedby studyingaPurkinje cellmodel andamodi- fiedventricularcellmodel.Wefindchaoticdynamicswithrespecttotheleakcurrentin thePurkinje cell model,and EADsand chaos withrespect toareduced fastpotassium currentand anenhanced calciumcurrentinthe ventricular cell model — featuresthat havebeenexperimentallyobservedandareknowntoexistinsomemodels,butarenew tothemodelsunderpresentconsideration.Wealsoinvestigatetherelatedmonodomain modelsofbothsystemstostudysynchronisationandthebehaviourofthecellsonmacro- scaleinconnectionwiththediscoveredfeatures.Themodelsshowqualitativelythesame behaviourtowhathasbeenexperimentallyobserved.However,forcertainparameterset- tingsthedynamicsoccurwithinanon-physiologicalrange.

© 2020TheAuthor(s).PublishedbyElsevierB.V.

ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Nowadays,mathematicalmodellingandnumericalsimulationsareessentialandstandardapproachestostudyandanal- yse realworld problemsand phenomenain life science. One major aim is the understanding of complexdynamics and behaviourofthesesystems.Forthispurpose,bifurcationtheoryhasproventobeaveryhelpfulandpowerfultoolinorder toinvestigatedynamicalsystemsandtheir(complex)dynamics,see[1–5]foran overview.Furthermore,numericalbifurca- tionanalysishasbecomeaprofitabletoolinthestudyof(forinstance)climate,neuronalandcardiacmodels.

Cardiacmusclecellscanexhibitcomplex patternsofoscillations likespikingandbursting,whichisrelatedtoioncur- rentinteractionsoftheconsideredcell.Asidefromnormalactionpotentialsofacardiacmusclecell,certainkindsofcardiac arrhythmiacanoccur. Thisincludesspecifictypesofabnormalheartrhythms, whichcanleadto suddencardiac death.In

∗ Corresponding author.

E-mail addresses: andreerh@math.uio.no (A.H. Erhardt), susanne.solem@ntnu.no (S. Solem).

1The authors have contributed equally to the content of this manuscript.

https://doi.org/10.1016/j.cnsns.2020.105511

1007-5704/© 2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

addition,irregular behaviour,such as(deterministic)chaos orchaotic earlyafterdepolarisations, hasbeenobservedinex- perimentalaswell asincomputational studies,see[6,7] andthereferencestherein.It isthereforehighlyinteresting and importanttounderstandthecomplexbehaviourandmechanismofsuch biologicalphenomena.Moreover, cardiacdynam- icsorheart rhythmscan be quitesensitive totheinfluence ofcertain drugs,whichhasbeen investigatedexperimentally butalso incomputational studies,seee.g. [8–10].Inlater years,the focushas beentocontinuously move towards inter- disciplinaryresearch, includingbiology,computer science,andmathematics, totackletheseissues.As aconsequence,the numberofexistingmathematicalmodelsbasedonexperimentaldataisalsocontinuouslyincreasing.

Thedevelopmentofagoodandprecisemathematicalmodelisessentialtodesignnumericalexperimentsforthestudy of cardiac dynamics,butalso forthe investigation ofthe influence ofcertain externaleffects such as drugsor oxidative stress. Tothisend, mathematicalanalysisisa keyto decodeoccurringphenomenaandtovalidate aderived modelinall details.Thenewlygainedinformationoftheconsideredmodel,canthenbeeitherusedtoimprovethemodelortoproceed withtheoriginalaspiration,e.g.theinvestigationofoptimalpropertiesofdrugs[9].

To this end, bifurcation theory has been utilised to investigate the dynamics of cardiac muscle cells in recentyears, see e.g. [11–14]. Continuing on this lineof research, this paper highlights how useful bifurcation theory can be forthe understanding of complexcardiac dynamicsand howit can be applied to find hiddenfeatures and dynamicsof cardiac singlecellmodelsdescribedbyordinarydifferentialequations(ODEs)throughnumericalinvestigations.

Intheend,itisthesynchronisationofalargegroupofcellsthatdecideswhetheracardiacarrhythmiaspreadsordies out. Forthisreason, a briefstudyof howthe micro-scale singlecellfeatures of thesemodels affectthe behaviour ofan ensembleofcellsatthemacro-scalelevel(cm)isprovided.Thisisdonebyanup-scalingoftheODEsystemtoaPDE–ODE monodomainmodel.

Alloftheabove willbedone byan in-depthmathematicalandnumericalinvestigationofthetwocardiac cellmodels introducedin[15,16],whereoneisamodelofaPurkinjecell,andtheotheramodelofahumanventricularcell.Inpartic- ular,wefindnewfeaturesoftheconsideredmodels,suchaschaosandearlyafterdepolarisations,andinvestigatehowthis affectsgroupsofcellsatthetissuelevel.

Wefindchaotic dynamicsinbothmodelsconsidered.Similarchaoticdynamicstowhatwediscovercanbeobservedin experiments,cf.[6].However,thedynamicsseemstoappearinanon-physiologicalrange.Thiseitherrequirestheimprove- mentofthemodelsortheexperimentalvalidation.

Inaddition, wediscoverEADsinthe humanventricular cellmodel.Thisbehaviour doesseemtobe within thephysi- ologicalrange[17],whichisavalidationofthemodelinthiscase.Finally,weshow thatthe(in)validityofthemodels in termsofbeingwithinthephysiologicalrangecarriesovertothesynchronisationeffectsinthecorrespondingmonodomain models.Thesefindingsclearlyhighlighttheadvantagesofbifurcationtheoryintheanalysisofcardiacmusclecelldynamics bydetectingunexpectedormaybenon-physiologicalbehaviourofthemodel.

1.1. Outlineofthepaper

InSection 2we start witha mathematicalandbiological descriptionofthe modelsandproblemsunderinvestigation.

Abriefbackgroundonthemodellingofcardiac musclecellsandonup-scalingtoamonodomainmodelatthetissuelevel (cf.[18–20])isprovided.Furthermore,we performastabilityanalysisoftheODEsystemfromNoble[15],andshowhow toextendthisanalysistothemacro-scalemonodomainmodel.Theapproachisexplainedindetailforthefourdimensional model[15].However,theansatzcanalsobeusedformorecomplexmodels,see[21,22].

Thestructures ofthesystemsin[15,16]aresimilar.Nevertheless,thebehaviouranddynamicsthattheydisplay canbe quitedisparateduetodifferentcomplexityandparametersettings.InSection3,weapplyanumericalbifurcationanalysisin ordertoderiveacompleteunderstandingofthedynamicsofthemodelfromNoble[15]withrespecttocertainparameters.

Theanalysisisthenextendedtothecorrespondingmacro-scalemonodomainmodel.BasedontheresultsinSection3,we thencontinueouranalysisbystudyingatendimensionalversionofthemodelfromBernusetal.[16]inSection4.Finally, wecloseourpaperwithadiscussioninSection5andthenaconclusion.

2. Biologicalandmathematicalbackground

The historyof mathematicalmodelling ofaction potentials(APs) ofexcitable biologicalcells likeneurons andcardiac cellsstartswiththefamousandpioneeringHodgkin–Huxley(HH)modelfrom1952[23].In[23],theauthorsestablisheda mathematicalapproachthatcanbeusedtomodelAPsofexcitablebiologicalcellsbyasystemofODEs.Thefirstmodelof acardiaccellistheNoblemodelfrom1962[15]ofagenericPurkinjecell.In1991,LuoandRudypublishedanionicmodel forcardiacactionpotentialinguineapigventricularcells[24].Inthelastdecades,therehasbeenanimmensedevelopment inthemodellingofcardiacmusclecells,seee.g.[25–28].Theseconductance–basedmodelsrepresentaminimalbiophysical interpretationofanexcitablebiologicalcellinwhichcurrentflowacrossthemembraneisduetochargingofthemembrane capacitance andmovement of ions across ion channels,cf. Fig.1. Ion channels are selectivefor particular ionic species.

Ingeneral, an APisa temporary,characteristicvariance inthemembranepotential ofanexcitablebiological cellfromits restingpotential.ThemolecularmechanismofanAPisbasedontheinteractionofvoltage-sensitiveionchannels.Thereason fortheformationandthespecial propertiesoftheAPisestablishedin thepropertiesofdifferentgroupsofionchannels intheplasmamembrane.Aninitialstimulusactivatestheionchannelsassoonasacertainthresholdpotentialisreached.

Fig. 1. (a) Scheme of a cardiac muscle cell [26] , where SR denotes the sarcoplasmic reticulum, I NaCa= Na ⁺/ Ca ²⁺exchanger current and I NaK= Na ⁺/ K ⁺pump current. (b) Physical system of a cardiac muscle cell. The dashed lines denote the ion currents, which are not included due to the lack of space.

Then, theseionchannels breakopen and/or up allowing an ioncurrent flow, which changes themembrane potential. A normalAPis always uniformandthe cardiacmuscle cell APis typicallydivided intofour phases:the restingphase, the upstrokephase,the(long)plateauphaseandtherepolarisationphase.Thismechanismisbasedonseveraldifferentcurrents.

Oneexample isthepotassium currentI_K whichis usually dividedintoa fast(I_Kr)anda slowcurrent(I_Ks), cf.schemein Fig.1(a).

Thiselectrophysiologicalbehaviourcanbedescribedbytheordinarydifferentialequation:

Cm

dV

dt =−Iion+I_stimulus,

whereV denotes the voltage(in mV) and t thetime (in ms), whileI_ion isthe sum ofall transmembrane ioniccurrents.

I_stimulusrepresentstheexternallyappliedstimulusandCmdenotesmembranecapacitance.

The model in [15] contains a singlepotassium current I_K, while the authors of [16] merged a fast (I_Kr) and a slow current(I_Ks) toderivetheir model,whichisbasedon thesystemin[29].Inthispaper,wewillslightlymodifythemodel fromBernus etal.[16],i.e.wereplaceI_Kby thecurrentsI_Kr andI_Ks fromPriebeandBeuckelmann [29].Furthermore,the differentioncurrentsmaydepend ondifferentgatingvariables,individual ionicconductancesG_i andNernstpotentialsE_i, i=Na,K,Ca,etc.,cf.Section2.1.

Wewanttohighlightthatcardiaccellmodelsusually havedifferenttimescalesandmayexhibit so-calledmixed-mode oscillations[30],cf.[31],and/orchaoticbehaviour,cf.[12,32–35],whichcanbelinkedtocertaincardiacarrhythmia.

Forinstance,ifthere are depolarisingvariations ofthemembrane voltage,we are speakingaboutafterdepolarisations (ADs).TheseADsaredividedintoearly(EADs)anddelayedafterdepolarisations (DADs).Thisdivisiondependsonthetim- ing obtaining the AP. EADsoccur eitherin theplateau orin the repolarisationphase ofthe APand are benefitedby an elongationoftheAP,whileDADsoccuraftertherepolarisationphaseiscompleted.EADsareresulting,forexample,froma reductionoftherepolarisingK⁺currentsoranenhancementinCa²⁺ currents,seee.g.[33].Triggersforthisarecongenital disordersofionchannelsortheingestionofmedicaments.Ingeneral,EADsareadditionalsmallamplitudespikes(mathe- maticallyspeakingmixed-modeoscillations),i.e.pathologicalvoltageoscillations,duringtheplateauorrepolarisationphase.

Theyarecausedby ionchannel diseases,oxidative stress ordrugs.Furthermore,thepresenceofEADsstrongly correlates withtheonsetofdangerouscardiacarrhythmias,includingtorsadesdepointes(TdP),whichisaspecifictypeofabnormal heartrhythmthat canleadtosuddencardiac death,see[17,36,37].Thus,itishighlyimportanttounderstandthecomplex behaviourofsuchbiologicalphenomena[38].

Finally,wewanttopointoutthatnotalloftheexistingcardiaccellmodelsmayincludecomplexdynamicssuchasEADs orchaos.

2.1. APurkinjecardiaccellmodel

First,wefocusonthemodelfromNoble[15],whichreadsasfollows:

dV

dt =−I_Na+I_K+I_L

Cm =:F, (1)

withthemembranecapacitanceCm=12 ^μF

cm² andtheioncurrentsI_Na (sodium),I_K(potasssium),andI_L (leak current),de- scribedbyI_Na=(^GNam³h+0.14)(^V−E_Na),

IK=

GK1n⁴+GK2exp

−V+90 50

+GK2

80 exp

_V+90

60

(

^V−EK

)

,

andI_L=G_L(^V−E_L),respectively. Theindividual ionicconductances aregivenby G_K1=G_K2=1.2 _cm^mS₂, G_Na=400 _cm^mS₂ and G_L=0.075 _cm^mS₂,whiletheNernstpotentialsaregivenbyE_K=−100mV,E_Na=40mV andE_L=−60mV.Furthermore,the differentgatingvariablesm,handnsatisfythedifferentialequation

dy

dt =ay

(

¹−y

)

−byy=ay−

(

^ay+by

)

^y=y_∞

(

^V

)

−y

τ

y

(

^V

)

^{, (2)}

whereyrepresentsthegatingvariablesm,handn,whiley_∞:=y_∞(^V)=_ay^a₊^y_by denotestheequilibriumofthegatingvariable yand

τ

^y:=

τ

^y(^V)= _ay+¹_by itsrelaxationtimeconstantwith

a_h=0.17exp

−V+90 20

, b_h= 1 1+exp

−^V+₁₀⁴²

, am= 0.1

(

^V+48

)

1−exp

−^V+₁₅⁴⁸

, bm= 0.12

(

^V+8

)

exp

_V+8

5

−1, an= 0.0001

(

^V+50

)

1−exp

−^V+50₁₀

, bn=0.002exp

−V+90 80

.

Noticethat thegatingvariables areimportantfortheactivation andinactivationofthedifferentioncurrents,whichis neededfortheioncurrentflows andtheresultingactionpotential. TheNoblemodel(1)describesthelonglastingaction andpace-maker potentials ofthe Purkinje fibresof theheart based on theHodgkin–Huxley formulation[23].While the sodiumcurrentisvery similartothe onefromHodgkinandHuxley[23],the potassiumcurrentdiffersinits formulation and thecalcium current ismissing. Nevertheless,the solutions to system(1) closely resemblesthe Purkinje fibre action andpace-makerpotentials.It isshownthatits behaviourin responsetoapplied currentsandtochangesinionic perme- abilitycorresponds fairly well withthatobserved experimentally.Furthermore,the Noblemodel(1) isoneofthe earliest mathematical modelsofcardiac APswhichis ableto produceacertain type ofcardiac arrhythmia,so-calledalternansin APduration (APD).ThesealternansinAPDisaresultofaperioddoublingbifurcationwithrespecttothecyclelength,see [39,40].Aperioddoublingbifurcationisacreationordestructionofaperiodicorbitwithdoubletheperiodoftheoriginal orbit.

In[15]theauthornumericallystudied theinfluenceoftheconductanceG_L oftheleakcurrentonthetrajectoryofthe Noble model(1).We willextendthisnumericalstudyby usingbifurcationanalysisto deriveamoredetailedinsightinto thebehaviourofthesolutionsto(1)byvaryingG_L,seeSection3.1.

2.1.1. Stabilityanalysis

The steady state or equilibriumof system (1) is determined by h≡ h_∞(V), m ≡ m_∞(V), n ≡ n_∞(V) andsolving the algebraicequation

F

(

V,h_∞

(

^V

)

,m_∞

(

^V

)

,n_∞

(

^V

))

=0. (3) This yields X_∞=(^V∞,h_∞(^V∞),m_∞(^V∞),n_∞(^V∞))^T, whereX_∞ is depending onseveral systemparameters, cf. system(1). Furthermore,theJacobianoftherighthandsideofsystem(1)evaluatedatX_∞isgivenby

J=

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎝

∂

^F

∂

^V

∂

^F

∂

^h

∂

^F

∂

^m

∂

^F

∂

ⁿ

1

τ

h

∂

^h∞

∂

^{V −}

1

τ

h

0 0

1

τ

^m

∂

^m∞

∂

^{V 0 −}

τ

^{1m 0}

1

τ

ⁿ

∂

ⁿ∞

∂

^{V 0 0 −}

1

τ

ⁿ

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎠

X∞

, (4)

whereweusedthefactthat

∂ ∂

^V

y_∞−y

τ

^y

y≡y∞

=

∂^y∞

∂^V

τ

^y−

(

^y∞−y

)

^∂τ_∂V^y

τ

y²

y≡y∞

= 1

τ

^y

∂

^y∞

∂

^{V .}

Note that the location andstability ofX_∞ is dependingon the different systemparameters, e.g. varying G_L changes the locationandthestabilityofX_∞,whilevarying C_m changesonlythestability.Todeterminethestabilityoftheequilibrium onehastocalculatethesolution(s)ofthecharacteristicpolynomial

det

(

^J−

λ

¹⁴

)

⁼

λ

^4+a1

λ

^3+a2

λ

^2+a3

λ

^{+a4=0, (5)}

i.e.theeigenvaluesofJ,where a1:= 1

τ

h + 1

τ

m+ 1

τ

n −

∂

F

∂

^V,

a2:= 1

τ

h

1

τ

m−

∂

^h∞

∂

^V

∂

^F

∂

^{h −}

∂

^F

∂

^V

+ 1

τ

n

1

τ

h−

∂

ⁿ∞

∂

^V

∂

^F

∂

^{n −}

∂

^F

∂

^V

+ 1

τ

m

1

τ

n−

∂

^m∞

∂

^V

∂

^F

∂

^m−

∂

^F

∂

^V

a3=1

τ

h

1

τ

^m

τ

^{n −}

₁

τ

^m+

1

τ

ⁿ

_F

∂

^{V +}

∂

^h∞

∂

^V

∂

^F

∂

^h

− 1

τ

ⁿ

∂

ⁿ∞

∂

^V

∂

^F

∂

^{n −}

1

τ

^m

∂

^m∞

∂

^V

∂

^F

∂

^m

− 1

τ

m

τ

n

∂

^F

∂

^{V +}

∂

ⁿ∞

∂

^V

∂

^F

∂

^{n +}

∂

^m∞

∂

^V

∂

^F

∂

^m

and

a4:=− 1

τ

h

τ

^m

τ

ⁿ

∂

^h∞

∂

^V

∂

^F

∂

^{h +}

∂

^m∞

∂

^V

∂

^F

∂

^m+

∂

ⁿ∞

∂

^V

∂

^F

∂

^{n −}

∂

^F

∂

^V

.

Inaddition,theRouth–Hurwitzcriterion impliesthatall characteristicexponents

λ

i,i=1,...,4havenegativerealpartsif andonlyiftheconditions

1=a₁>0,

2=a₁a₂−a₃>0,

3=a₃·

2>0and

4=

3−a²₁a₄>0

holdtrue,see[3,4,41].Furthermore,ifall Hurwitzminorssatisfyi >0fori=1,· · ·,n−1andⁿ=0,wherendenotes thedimensionofthesystem,weknowthatthesystemexhibitsanAndronov–Hopfbifurcation.Using4=0,weget

0=

λ

^4+a1

λ

^3+a2

λ

^2+a3

λ

^+a4=

λ

^4+a1

λ

^3+a2

λ

^2+a3

λ

^+a1a2_a^−a3

1

a3

a1

=

λ

^2+a_a³

1

λ

^2+a1

λ

^+a1a2_a^−a3

1

,

i.e.theequilibriumhasapairofpurelyimaginaryeigenvalues

λ

1,2=i

ω

0 with

ω

0=a₃

a₁ >0andtwofurthereigenvalues

λ

³,4= −a₁±

a²₁−4^a1a_a^2−a3

1

2 =−

^1±

²1−4²

1

2 .

Ingeneral,anAndronov–Hopfbifurcation correspondstothebirthofalimitcycle,whentheequilibriumchangesstability viaa pairofpurelyimaginary eigenvalues.Usually,an Andronov–Hopfbifurcation isconsidered asa triggertooscillatory behaviourindynamicalsystemsandmaycausenormalAPandcardiacarrhythmiainacardiaccellmodel.

Incasethatan=0,thenthesystemexhibitsafoldorsaddle–nodeorlimitpointbifurcation,i.e.theequilibriumhasat leastoneeigenvaluewhich isequaltozero.Alimit pointbifurcation isacollision anddisappearance oftwoequilibriain dynamicalsystems,whichisaturningpoint.

2.2.AmonodomainmodelofthePurkinjecardiaccellmodel

Besides the study of cardiac cell models an important focus is the behaviour anddynamics of severalcells, i.e. the dynamicsonamacro-scale(cm),wheremanycellsareconnectedtogetherandinteractingwitheachother.Tothisend,we considerthefollowingmonodomainmodel,i.e.extensionoftheODEmodel(1)toaPDE–ODEmodelincludinganadditional diffusionterm:

Cm

∂

^V

∂

^{t =}

λ

1+

λ χ

¹

∇

·

(

^Mi

∇

^V

)

⁻

(

^INa+IK+IL

)

∂

^y

∂

^{t = y∞}

τ ⁽

^yV

(

^V

⁾

⁻

)

^{y, y=h,m,n in}

, (6)

0=

ν

^·

(

^Mi

∇

^V

)

^on

∂

^,

where^{isa boundeddomain,M}i denotesthe intracellularconductivitytensor,

λ

^{theextra- to}intracellular conductivity ratio,

χ

^{isthemembranesurfacearea perunitvolumeand}

ν

^{istheunitnormal,cf.[42–44].System(6)isamonodomain}

model,meaningthatequalanisotropyrates,i.e.Me=

λ

^Miin ^mS_cm,areassumedinthemorecomplexbidomainmodel.Here,

λ

∈RisconstantandMe denotestheintracellular conductivitytensor.Furthermore,weuse ₁^λ_+λ^M_χⁱ =360¹ mS,whichisthe diffusionconstantoriginallyusedin[16],unlessotherwisestated.

Forabetter understandingofthebehaviourofcellsona macro-scalelevel,wefollowtheapproachfromLietal.[45], ChenandJiang[46],CaoandJiang[47]toderivealinearisedsystemofmodel(6).Firstofall,weknowthatonrectangle-like domains^:=[0,1]× · · · ×[0,n]⊂RⁿtheeigenvaluesandeigenfunctionsoftheNeumannproblem

₋

_u

k

(

^x

)

=

μ

ku_k

(

^x

)

^x∈

,

ν

^·

( ∇

^uk

(

^x

))

^{=0 x∈}

∂

⁽⁷⁾

are

μ

k⁽ⁱ⁾=

π

^k

i

2

and u⁽_kⁱ⁾=cos

π

^kx

i

k=0,1,2,..., i=1,...,n,

see [48].In general, theeigenvalue of the spectral problem (7)can be derived by multiplying the first equation ofsys- tem(7)byu_k,integratingover^{,usingGreen’sformulaandapplyingtheboundarycondition.Then,onegetsfortheNeu-} mannproblem(7)thefollowingeigenvalues:

μ

k=

∇

^uk

²L²()

^uk

²L²() .

Now,weconsiderthelinearisedsystemofthemonodomainEq.(6)aroundanequilibriumX_∞oftheODEsystem(3).It hastheform

∂ ∂

^t

⎛

⎜ ⎝

V h m n

⎞

⎟ ⎠

^=D

⎛

⎜ ⎝

V h m n

⎞

⎟ ⎠

^+J

⎛

⎜ ⎝

V h m n

⎞

⎟ ⎠

^{, (8)}

whereD denotesthe4×4diffusionmatrixwithalmosteverywherezeroentriesexceptthefirstone,whichis D11=

λ

1+

λ

M_i

χ

1 Cm,

whileJ istheJacobianasstatedin(4).DefinethelinearoperatorLas

L

⎛

⎜ ⎝

V h m n

⎞

⎟ ⎠

^:=D

⎛

⎜ ⎝

V h m n

⎞

⎟ ⎠

^+J

⎛

⎜ ⎝

V h m n

⎞

⎟ ⎠

^.

Then,considerthefollowingcharacteristicequation

L

⎛

⎝ ψ

1

..

ψ

.4

⎞

⎠

=

μ

⎛

⎝ ψ

1

..

ψ

.4

⎞

⎠

,

where(

ψ

1,,

ψ

4)^TistheeigenfunctionofLcorrespondingtotheeigenvalue

μ

^.Thus,let

⎛

⎝ ψ

¹

..

ψ

.4

⎞

⎠

=^∞

k=0

⎛

⎜ ⎝

V_k h_k m_k n_k

⎞

⎟ ⎠

^cos

π

^kx

,

whereV_k,h_k,m_kandn_karetime-dependentcoefficients.Wecanthenconcludethat ∞

k=0

Jk

⎛

⎜ ⎝

V_k h_k m_k nk

⎞

⎟ ⎠

⁼

μ

^∞

k=0

⎛

⎜ ⎝

V_k h_k m_k nk

⎞

⎟ ⎠

^,

with

Jk=

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

Jk11

∂

F

∂

^h

X∞

∂

F

∂

^m

X∞

∂

F

∂

ⁿ

X∞

1

τ

h

∂

^h∞

∂

^{V −}

1

τ

h

0 0

1

τ

^m

∂

^m∞

∂

^{V 0 −}

1

τ

^{m 0}

1

τ

ⁿ

∂

ⁿ∞

∂

^{V 0 0 −}

1

τ

ⁿ

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠

, (9)

where

Jk11=

∂

F

∂

^V

X_∞

−

π

^k

2

λ

1+

λ

^M

χ

^iC1m

k=0,1,2,3,....

Keepinmindthatsystem(7)haseigenvalues 0=

μ

^0<

μ

¹⁼

π

2

<

μ

²⁼⁴

π

2

<

μ

³⁼⁹

π

2

<· · · −→∞. Hence,thelinearisedsystem(8)hasinfinitelymanyJacobiansJk.

We continue by deriving an ODE model to represent the behaviour of the linearised system (8) for each mode k= 0,1,2,3,... in closeproximity to theequilibrium V_∞. This is done by ensuring that the resulting system hasthe same equilibriumastheNoblemodel(1)andtheJacobianJk.

Then,weareinasituationwherewecananalysethebehaviourofasinglecellonamacro-scaleincludingtheinfluence ofthe diffusionterm(dependent onk) of themonodomain model(6). Thisallows usto gain intuition onhow thecells interact.TheresultingODEsystemreadsasfollows:

d dt

⎛

⎜ ⎝

V_k hk

m_k n_k

⎞

⎟ ⎠

⁼

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎝

−I_Na+I_K+I_L Cm −

π

^k

2

λ

1+

λ

M_i

χ ⁽

V_k−V_∞

)

Cm

(

^h∞

(

^Vk

)

−h_k

)

/

τ

h

(

^Vk

) (

^m∞

(

^Vk

)

−mk

)

/

τ

m

(

^Vk

)

(

ⁿ∞

(

^Vk

)

−n_k

)

/

τ

ⁿ

(

^Vk

)

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎠

, (10)

whereV_∞ is an equilibriumforEq.(3). Wehavedesignedthe locationof theequilibriumofsystem(10) tobe thesame asfor theNoble model(1). However, considering the stability analysisfrom Section 2.1,the stability ofthe equilibrium maybedifferent.The firstentryoftheJacobian J in (4)isreplacedby Jk₁₁.Thus, alsothecoefficientsa_j, j=1,...,4of thecharacteristicpolynomial(5)arechanged.Obviously,thiswillaffectthestabilityofthe system(10)dependent onthe parameters

λ

^,Mi,

χ

^,andk.Thisindicatesthatthecellularbehaviourofmodel(1)isnot(necessarily)one-to-onetransferred tothebehaviour anddynamicsofthemonodomainEq.(6).Nevertheless,itisagoodstartingpointtostudythedynamics ofa singlecellbefore extendingthe analysistothe macro-scale. Dohowevernote that the modek=0doesgive usthe samedynamicsasthatoftheODEsystem(1).

Inthediscrete setting,ifwechoosek= _h¹2,wherehdenotesthegridsize,andk=_h¹2 isthebiggest eigenvalueforthe discretelaplacian,theterm

−

( π

^k

)

²

λ

1+

λ

M_i ²

χ

1 Cm =−

π

²

h⁴

λ

1+

λ

M_i ²

χ

1 Cm =−

π

²

h⁴ 1 ²

1 4320

mScm²

μ

^F

in(10)tendsto−∞andblows upash→0,whichshouldstabilise(10)forlargemodesk(orrefinedgridsizeh).Looking atthecoefficientsa_j, j=1,...,4ofthecharacteristicpolynomial(5),wecanseethat a1,a2 anda3 willtendto∞,while a₄ willtendto−∞ashtendstozeroandktendsto∞.Thisimpliesthat

1=a₁>0,

2=a₁a₂−a₃>0,

3=a₃·

2>0and

4=

3−a²₁a₄>0

andthus,thenumericsofthePDEwillstabilisefordecreasinggridsize,whichistobeexpected.

Thesameanalysiscanbeusedon2Ddomains=[0,]×[0,∗],,^∗ >0,byconsideringtheslightlydifferentterm

−

( π

^k

)

²

λ

1+

λ

^M

χ

^iC1m

1 ²+ 1

²_∗

.

Asimilarmodificationappliestoa3Dcubeorcuboid.Foramoregeneralgeometryonecanexpectthattheadditionalterm derivingfromthelinearisationofthemonodomainmodel(6)ismorecomplicated.However, thegeneraldiscussionabove isexpectedtoremainvalid.

2.3.Aventricularcardiaccellmodel

Asmentioned,wefirstinvestigatetheNoblemodel(1).Indeed,wewillseethatthismodelhassomelimitations.There- fore,wewillalsostudyahumanventricularcellmodel,whichismoreadvancedduetothenumberofincludedioncurrents.

Thedescriptionofthesystemin[16]forepicardialcellsissimilartosystem(1),butitcontainsmoreioncurrentsandreads asfollows:

Cm

dV

dt =−Iion+Istimulus, (11)

whereI_stimulusdenotesanexternalstimulusand

Iion=ICab+INaCa+INab+ICa+IK+INaK+INa+IK1+Ito

isdependingonthefastsodiumcurrentI_Na=G_Nam³

v

²₍V−E_Na),theslowcalciumcurrentI_Ca=³₅G_Cad_∞(^V)^f(^V−E_Ca),the transientoutward currentIto=Gtor_∞(^V)^to(^V−Eto), thedelayed rectifier KcurrentI_K andthe inward rectifier K1current IK1=GK1K1_∞(^V)(^V−EK),respectively.

In [16] the authors studied system (11) witha delayed rectifier current I_K=G_Kx²(^V−E_K), while we are considering the delayed rectifier current I_K=I_Kr+I_Ks from Priebe andBeuckelmann [29]. Here, the currentI_Kr=G_Krxrrik(^V)(^V−E_K) denotestherapidlyactivatingcurrent,whileI_Ks=G_Ksx²_s(^V−E_K)^{istheslowlyactivatingcurrent.Includingthefastandslow} potassiumcurrentwillmakesthedynamicsmorerealistic.

Furthermore,system(11)contains thebackground currentsI_Cab andI_Nab,the Na⁺/Ca²⁺exchanger currentI_NaCa and the Na⁺/K⁺ pumpcurrentI_NaK,cf.Fig.1(a).Noticethatthesystemisdependingon9gatingvariables,i.e.m,v,d,f,r,to,xr,xs

andK1,satisfyingthedifferentialEq.(2),whered,randK1areassumedtobeequaltotheirsteadystates.Wewillconsider allgatingvariablesasstatevariables,thereforethedimensionofthesystemis10.

Moreover, we useCm=1.534 _cm^μF₂,cf. [29], whilethe individual ionic conductancesare given by G_Na=16 _cm^mS₂,G_Ca= 0.064 _cm^mS₂,Gto=0.3 _cm^mS₂,G_Kr=0.015 _cm^mS₂,G_Ks=0.02 _cm^mS₂ andG_K1=2.5 _cm^mS₂.TheequilibriaandtheJacobian aresimilarly determined asforthe Noblemodel (1). Theonly difference isthat we have6ODEs more toconsider. Thisincreasesthe 4×4matrixJ toa10×10matrix.FollowingthesameapproachasinSection2.2,wecanderivefromthemonodomain modelrelatedtosystem(11)withI_stimulus=40_cm^μA₂ andadurationof2s,i.e.

Cm

∂

^V

∂

^{t =−}

λ

1+

λ

1

χ ∇

^·

(

^Mi

∇

^V

)

^−Iion+I_stimulus,

∂

^y

∂

^{t =}

y_∞

(

^V

)

^−y

τ

^y

(

^V

)

^{, y=m,}

v

,d,f,r,to,xr,xs,K1 in

⁽¹²⁾

withNeumannboundarycondition

ν

·(^Mi

∇

^V)=0on

∂

^{,thefollowingODEsystem}

d dt

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎝

V_k m_k

v

k

dk

f_k r_k tok

xrk

xsk

K1k

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎠

=

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

−Iion−I_stimulus

Cm −

π

^k

2

λ

1+

λ

Mi

χ ⁽

V_k−V_∞

)

Cm

(

^m∞

(

^Vk

)

−mk

)

/

τ

m

(

^Vk

) ( v

_∞

(

^Vk

)

−

v

k

)

/

τ

v

(

^Vk

) (

^d∞

(

^Vk

)

−d_k

)

/

τ

d

(

^Vk

) (

^f∞

(

^Vk

)

−f_k

)

/

τ

f

(

^Vk

) (

^r∞

(

^Vk

)

^−rk

)

^/

τ

^r

(

^Vk

) (

^to∞

(

^Vk

)

−tok

)

/

τ

to

(

^Vk

)

(

^xr∞

(

^Vk

)

−xrk

)

/

τ

^xr

(

^Vk

) (

^xs∞

(

^Vk

)

−xsk

)

/

τ

^xs

(

^Vk

) (

^K1∞

(

^Vk

)

−K1_k

)

/

τ

^K1

(

^Vk

)

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠

(13)

whereV_∞ theequilibriumofthevoltageVofsystem(11).Notethatthestabilityanalysisforsystem(10)andtheprevious discussionalsoholdstrueforsystem(13).

2.4. Numericalmethods

Forour simulations wewill useMATLAB R2019b andtheode solver ode15switha relative toleranceof 10⁻¹³ andan absolute tolerance of 10⁻¹⁸. For the monodomain models, the pdepe solver is used. Moreover, as initial values for sys- tem(1)wewilluseV₀=−79.04mV,h₀=0.81,m₀=0.045andn₀=0.52,whileforthesecondmodel(11)weutiliseV₀=

−93.3701mV, m₀=0.0004,

v

0=0.9990, f₀=0.8797,xr₀=0.0042,to₀=0.9999,d₀=0.0000,r₀=0.0000,K1₀=0.0419 andxs₀=0.0912,aslong we donot specify anythingelse. The desiredbifurcation diagramswill be derived utilisingthe MATLABtoolboxesMATCONTandCL_MATCONT[49–51],whicharenumericalcontinuationpackagesforinteractivebifurca- tionanalysisofdynamicalsystems.

3. DynamicsoftheNoblemodel

In[15],theauthormentionedachangeinthedynamicsofsystem(1)byvaryingtheleakconductanceG_L.Wewillshow thatchangingG_L hasinfluenceontheperiodoftheAPaswellasifthesystemconvergesintoastableequilibriumornot.

InFig.2,thetrajectory ofsystem(1)isgivenfortwodifferentvaluesoftheleakconductanceG_L,i.e.G_L=0.075 _cm^mS₂ and G_L=0.18 _cm^mS₂. InFig.2(a), system(1)revealsa periodictrajectory representingtwo APsofa cardiac singlecell,while in Fig.2(b)thetrajectory needs certainamountoftime to reachastableperiodic patternalso representingAPs.Thisshows that system (1)is sensitive with respectto G_L, which mayhave an influence on the amplitudeandthe period T ofthe trajectorygivenbyV(^t+T)−V(^t)=0.Inaddition,theinitialvaluealsohasaninfluenceontheappearingpattern(atleast locally).

Fig. 2. Different action potentials: Simulation of system (1) (default setting) with two different values of the leak conductance G L, where these values are chosen according to the observation in [15] .

Fig. 3. Bifurcation diagram in 2D: Projection onto the ( G L, V )-plane showing the equilibrium curve and the first two limit cycle branches.

3.1. Bifurcationanalysisofsystem(1)withrespecttoG_L

We now systematically investigate the behaviour described above utilising numerical bifurcation analysis as done in [31,52].Ourfirststepistodetermineabifurcationdiagramwithrespecttotheleakcurrent,i.e.wechoosetheleakconduc- tanceGLasthebifurcationparameter.NoticethatwewillonlyconsiderpositivevaluesofGL,sincetheyarethephysiologi- callyrelevantones.Mathematicallyandnumericallywecaneasilyextendthebifurcationdiagramalsotonegativevaluesof G_L.

Westartby determiningtheequilibriumcurve, i.e.we calculatefordifferentvaluesofG_L thecorrespondingequilibria, whichgivesusthedesiredequilibriumcurve.Moreover,wedeterminethestabilityofeachequilibria.Theequilibriumcurve wecaneasily calculatewithaself-writtenroutine.However, toderiveadetailedenoughbifurcationdiagram efficiently,it isadvisabletousearobustcontinuationalgorithmastheonementionedinSection2.4.

Thebifurcationdiagramofsystem(1)relatedtoG_Lexhibitsanunstable(blackdashedline)andastable(blacksolidline) equilibriumbranch,seeFig.3.Theequilibriumcurve changesstabilityviaasupercritical Andronov–Hopfbifurcation(blue dot,G_L≈0.200883_cm^mS₂)withanegativefirstLyapunovcoefficient.FromthesupercriticalAndronov–Hopfbifurcationastable limitcyclebranch(solidblueline)bifurcates,whichbecomesunstable(dashedblueline)viaaperioddoublingbifurcation (solidredsquare,G_L≈0.187785_cm^mS₂).Then,thelimitcyclebranchgainsstabilityagainviaalimitpointofcyclebifurcation (solidgreensquare,G_L≈0.193546_cm^mS₂).Furthermore,thislimit cyclebranchcontainsasecondperioddoublingbifurcation, whichisalsoconnectedtothefirstoneviaasecondlimitcyclebranch,cf.Figs.3and4.

ThebifurcationdiagraminFig.3onlyincludesthefirsttwolimitcyclebranches,asincludingfurtherdetailswouldnot beparticularlyvisible.Indeed,furtherbranchesexistaswewillseebelow.

Togetherwith Figs. 4 and5,a nice graphical explanation ofthe observations in[15] appears. We canidentify values ofG_L for whichsystem (1) oscillatesor hasa stable equilibrium. Even more,we can determine the maximal amplitude

Fig. 4. Bifurcation diagram in 2D: Zoom of Fig. 3 also including the third limit cycle branch.

Fig. 5. Bifurcation diagram in 3D: Projection onto the ( G L, n, V )-plane including several trajectories of system (1) for different values of G L.

ofan oscillationcorresponding toG_L.Inaddition, wealsoget thecorresponding periodofeachlimit cycle.We havethat theperiodTforG_L=0.075 _cm^mS₂ isapproximately564.1345ms,whileforG_L=0 _cm^mS₂ we haveT≈839.5015ms andforG_L= 0.18 _cm^mS₂ wehaveT ≈ 324.2749ms.Thus, we knowthatthe patternofthetrajectory repeats fasterifG_L increasesbefore itconvergesintoastableequilibriumforG_L toolarge,whichexplainstheobservationsfromNoble[15].Furthermore,from Fig.3weobservethatthemaximalamplitudeisdecreasingwhileG_L increases.

InFig.4athirdlimitcyclebranchisalsoincludedandwefocusonasmallerrangeofG_L values,whereinterestingdy- namicsmayappear.Thenextadditionallimitcyclebranchesbehaveverysimilarlytothethirdone,i.e.theyarebifurcating fromaperioddoublingbifurcation,theylosestabilityviaafurtherperioddoublingbifurcation,andstayunstableuntilthey convergeintoaperioddoublingbifurcationofthepreviouslimitcyclebranch.Onlythesecondlimit cyclebranchbehaves slightlydifferent— itbifurcatesfromthefirstperioddoublingbifurcation(G_L≈0.187785_cm^mS₂),becomesunstableviaalimit point ofcyclebifurcation, gains stabilityagainvia asecond limit point ofcycle bifurcation andfinally,loses stabilityvia a furtherperioddoublingbifurcation (G_L≈0.187308_cm^mS₂). Then,itstays unstableuntilit convergesintoaperioddoubling bifurcationofthefirstlimitcyclebranch,cf.Fig.4.

Startingfromthesecondlimitcyclebranch,system(1)exhibitsastableperioddoublingcascade,whichisusuallyaroute tochaos,cf.e.g.[35].However,inthebifurcationdiagraminFig.5,whereseveraltrajectories(redlines)arehighlightedfor differentvaluesofG_L,noirregularorchaoticbehaviournoradditionaloscillationscanbeseen.

The reasonis quite simple: The system exhibits a stable limit cyclebranch fromthe first limit cycle branch, which

“surrounds” the interesting dynamics of the Noble model(1),cf.Figs. 3and4.Hence,for the standard initial condition we do notget anysudden changein thedynamics.The trajectory will jumpon tothe “outer” stablelimit cyclebranch and staysthere.