Comparing Two Parallel Genetic Algorithms For The Inverse Problem Associated To The Cardiac

AssociatedToTheCardiacBidomainEquations

CarolinaR.Xavier

ViniciusdaF.VieiraRodrigoWeberdosSantos

DavesM.S.Martins

DepartmentofComputerScience,UniversidadeFederaldeJuizdeFora

JuizdeFora,MinasGerais,Brazilrodrigo.weber@ufjf.edu.br

Abstract

Inthisworkweuseacomputationalhumanleftventricu-larwedgetosimulateregionsofabnormalintra-andextra-cellularconductivitiesthatmimicsomeknownpathologicalconditions.Weimplementandcomparetwogeneticalgo-rithmsthataimonestimatingthedistributionofintra-andextra-cellularconductivities,bycomparingcardiacsimu-lationstosomegiventransmuralelectrograms.Themeth-odsweredevelopedfordistributedsystemsandtheresultswereobtainedinaclustercomposedof8computersinter-connectedbyafastnetworkswitchingdevice.Theresultssuggestthattheproposedmethodsareabletocorrectlyestimatebothintra-andextra-cellularconductivitydistri-butionsfromtransmuralelectrogramswithanaccuracyof30%.Inaddition,theparallelimplementationbasedonthesteady-statemethodachievedbetterspeedupresultsthanthetraditionalgenerationalgeneticalgorithm.

1.Introduction

Computationalmodelingisausefultoolfortheinvesti-gationandcomprehensionofthecomplexbiophysicalpro-cessesthatunderliecardiacelectrophysiology[7,3].Mod-erncomputationalmodelsofthecardiactissuetakeintoac-countdetailedpropertiesofcardiacfiberandsheettrans-muraldistributions,extra-andintra-cellularconductivities,andthedynamicsofseveralioniccurrents.Thesecomputa-tionalmodelsareabletosimulatecardiacelectricalactivityandthecorrespondingsurfaceelectrocardiograms(ECGs)similartothoseobtainedinanimalexperiments[10].

Studiesusingone-dimensionalandtwo-dimensionalcar-diacmodelshavesuccessfullysimulatedECGsthatresem-bleseveralcardiacpathologiesbyalteringspecificmodelparameters,suchasconductivityvaluesandioniccurrent

densities.Unfortunately,themulti-physicsandmulti-scalecomplexityassociatedtothephenomenonofcardiacelec-tricalwavepropagationtranslatesmathematicallytolargenon-linearsystemsofpartialdifferentialequations.Theverylongexecutiontimesandcomputermemorydemandsofsuchmodelshavehistoricallylimitedtheuseofcardiacsimulationtothesolutionofforwardproblems,wheredif-ferentscenariosarestudiedbysolvingthePDEsystemswithdifferentparameters,initialandboundaryconditions,andindiversedomainsanddimensions.Inthiswork,weinvestigateparallelalgorithmstosolvetheinverseproblemassociatedtocardiacelectrophysiology.Inparticular,weimplementandcomparetwogeneticalgorithmsthataimonestimatingthedistributionofintra-andextra-cellularcon-ductivities,bycomparingcardiacsimulationstosomegiventransmuralelectrograms.Therearestrongevidencesthatthesephysiologicalvariableschangeundermanypatholog-icalconditions,suchasdilatedcardiomyopathy(DCM),is-chemiccardiomyopathy(ICM)andmyocarditis(inflamma-torycardiomyopathy)[5],acuteischemia[4]andchronicinfarct[1].

Thegeneticalgorithmstestparameter-candidatesbycomparingcardiacsimulationstosomegivenECGs.Here,weassumethatalltheotherparametersofthemathemat-icalmodelarefixedandknown.Inaddition,theaccu-racyoftheestimationistestedwithsimulatedECGs,forwhichthecorrectconductivityvaluesareknown.Twodif-ferentparallelgeneticalgorithmswereimplementedandcompared:aclassicalgenerationalgeneticalgorithm;andonebasedonsteady-statemethod.Preliminaryresultswereobtainedforasmalltwo-dimensionalcardiacmodel,inaclustercomposedof8computersinterconnectedbyafastnetworkswitchingdevice.Theresultssuggestthatthepro-posedmethodisabletocorrectlyestimateconductivityval-uesfromgivensurfaceECGs.Inaddition,theparallelimplementationbasedonthesteady-statemethodachieved

betterspeedupresultsthanthetraditionalgenerationalge-neticalgorithm.

2.Cardiacsimulation

Cardiacmodelingisapowerfulltooltostudyandcom-prehendthecomplexbiophysicalprocessesofcardiacelec-trophysiology[7].Fromsimulations,moredetailsaboutthecardiacelectricconductionmaybeobtained,andsomecardiacpathologiesmaybeinvestigated.ThesetofBido-mainequationsiscurrentlyoneofthemostcomplexmathe-maticalmodelthatsimulatecardiacelectrophysiology.TheBidomainmodelisbasedinthedivisionofthecardiactis-sueintotwodomains:intracellularandextracellular.Thesedomainsarecontinuous,andeachpointinthecardiacmus-cleisconsideredtobeinbothdomains,whicharelinkedbyasemipermeablemembrane,thecardiaccellmembrane.Both,theextracellularandintracellularmatricesarecon-sideredpurelyresistive.Thisassumptionisreasonabletotheextracellularmatrix,inwhichthecurrentmayflowinacontinuousspaceamongthecells,whileintheintracellu-larmatrixthemodelisjustifiedbytheexistenceofspecialproteinarrangements,thesocalledgapjunctions,whichconnecttheinteriorofneighboringcells.Thefunctionalcouplingofthesetwodomainsisaccomplishedusingnon-linearsetsofequationswhichdescribethetransmembraneioniccurrentsthataregeneratedacrossthesarcolemmaofthehumanventricularmyocyte.Thishasbeendescribedindetailbyten-Tusscheretal.[6].Inaccordancewithpub-lishedcellularelectrophysiologicaldata,threedistinctven-tricularmyocytephenotypeormathematicalmodelswereutilizedinthesecalculations:epicardial,Mandendocardialcells.Thenumericalsolutionofthislargenon-linearpartialdifferentialsystemyieldsspatialdistributionsandtemporalcharacteristicsoftheextracellularpotential(φe),intracellu-larpotential(φi)andtransmembranepotential(Vm).Santosetal.[8]havedevelopedthenumericalmethodsforeffi-cientlysolvingthisbidomainmathematicalmodel.

Wehavedevelopedamathematicalmodelforthepur-poseofsimulatingelectrophysiologicalphenomenaarisingfromasegmentorwedgeofhumanleftventriclewhichisassumedtobeinaperfusionmediumorbath(apassiveandisotropicconductor).Thesizeofthisventriculartissuewaschosentobeapprox.1.5cm(endo-toepicardium)by1.5cm(apextobase).AsshowninFigure1,theendocardial,epicardial,apexandbasesurfacesinterfacewithahomo-geneousextracellularmedium(bath),yieldinganoveralltissue-bathdimensionof3.0x3.0cm.

Allbidomainparameterswerebasedonthosereportedinpreviouswork.[7]3Dorthotropicconductivitytensorsthatvaryinspacewereusedtoreplicatethelaminarfiberstruc-tureoftheheart,whereσeandσistandfortheextracellu-larandintracellularconductivitytensors,respectively.The

cardiactissueconductivityvaluesfromtheliteraturehavebeenuniformlyrescaledtomatchthereportedapex-to-base(70cm/s)andtransmuralconductionvelocities(45cm/s)inthemammalianventricles.Thebathconductivitywassetto20mS/cm.Thecapacitanceperunitareaandthesurfacearea-to-volumeratioaresetto2µF/cm2and2000/cm,re-spectively.Thespatialandtemporaldiscretizationstepsofthenumericalmodelaresetto150µmand50µs,respec-tively.Allsimulationswerecarriedoutforaminimumof350msafterasinglecurrentstimuluswasintroducedataselectedendocardialsite.

Theextracellularpotentialsdve,i.e.thesignalswhichwouldbesensedbytransmuralleadswerecalculatedbytak-ingthedifferenceofthesimulatedextracellularpotentials(φe)attheendocardialandepicardialboundarypoints.ThepositionsofsuchvirtualelectrodesareillustratedinFigure1.

Wemodifiedtheextracellularandintracellularconduc-tivityvaluesofaspecificandsmallrectangularregionof0.3x0.7cmneartheendocardial,asillustratedinFigure1.Weconsidertheextracellular(intracellular)tensor,σi(σe),tobeisotropicallymodifiedbyascalarα(β).Thereforetheextracellularandintracellulartensorsinthethespec-ifiedregionareασiandβσe,respectively.Thiswaytheanysotropyorratherorthotropyofthetissueiskeptunmod-ified.Thesimulationwecarriedoutreproducesthecaseofachronicinfarctedregionbytaking(α,β)=(0.1,3.0),i.e.,reproducingthereductionofintracellularspaceandconduc-tivityduetomyocyteloss,andcorrespondingincreaseofextracellularspaceandconductivity,aspreviouslyreportedin[4].

3.Inverseproblemandgeneticalgorithms

Thefocusofthisworkisthesolutionofaninverseprob-lemorparameterestimationforthebidomainequations.Particularly,weconsidersomegivenseriesofextracellularpotentials,ECGs,and,fromthis,ourgoalistoestimatethetwoparameters(α,β)whichcharacterizeaspecificregionofthecardiactissue.WetakeasECGmeasurementstheresultsofachronicinfarctsimulationconsideringαwithavalueof0.1andβwithavalueof3.

Theinverseproblemcanbeformulatedas:

󰀁(minα,β)

F(α,β)(1)

󰀂np󰀂nt

F(α,β)=

i=1

j=1

(vde(α,β)(i,j)−vdeo(i,j))2

npnt

(2)

whereiisfrom1tonp,npisthenumberoftransmuralderivations.InourcaseasillustratedinFigure1np=4.vdeo(i,j)istheithobservedtransmuralelectrogramat

pointjdt,wheredtisthesamplingrateortimediscretiza-tion,jisfrom1tont,ntisthetotalnumberofdiscretiza-tions.vde(α,β)arethesimulatedtransmuralelectrogramsforagiven(α,β).

Wesolvethisoptimizationproblemusingthesearchmethodcalledgeneticalgorithm.ThisclassofalgorithmsisbasedonthenaturalselectiontheoryofCharlesDarwin[9].Thecardiacsimulatorisusediterativelyandseveraltimesbyaparallelgeneticalgorithmtosearchthebestpairofparameters(α,β)whichreproducesthegivenECGtimeseries.

Thegeneticalgorithmsearchesfortheparameters(α,β)thatminimizethefitnessfunctionFgivenbyEquation2.Eachcandidate(α,β)istreatedasachromosomeandisrepresentedbyawordof14bits.Thefirstsevenbitsrep-resentαandthelastsevenbitsrepresenttheparameterβ.Ourgeneticalgorithmbeginswithapopulationof21chro-mosomesrandomlychosen.Foreachcadidate(α,β),thegeneticalgorithmcallsthecardiacsimulatorwhichgener-atesthesetoftransmuralelectrogramsvde(α,β).Thege-neticalgorithmtakesthesimulatedresultsandcomputesthefitnessofthecandidate(α,β)bycastingEquation2.Thegenerationsofnewcandidatesisbasedontheprocessesofselection,crossoverandmutation,asinspiredbythenaturalselectiontheoryofDarwin[9].

Figure1.Measuringpointsdistributioninthecardiactissue

3.1.Master-SlaveParallelArchitecture

Geneticalgorithmsaregenerallyconsideredasnaturalparallelalgorithms,sincetheyhaveahighlyparallelizablecomputationalstructure.Analysingthegeneticalgorithm

structure,wecanconclude:

•eachpopulationindividualmaybeevaluatedindepen-dentofanyotherfactor.

•eachoperatorandgeneticoperationisindependent,be-causetheymaybeappliedinanyorder,sequentiallyornot,toanyindividual.

Ourgeneticalgorithmsweretargettodistributedenvi-ronmentsandtheirparallelimplementationsfollowthetra-ditionalmaster-slaveparalleldecomposition.Acentralma-chine,themaster,isresponsiblefortheexecutionofthegeneticalgorithm.Thismachineexecutesthetasksofse-lectionofnewcandidates,crossoverandmutation.Eachslavemachineisresponsibleforexecutinganinstanceofthecardiacsimulatorforagivencadidate(α,β),aswellasforcalculatingthefitnessfunctionofthatcandidate.Oneimportantfeatureofthiskindofparallelimplementationisthatitbehaveslikethesequentialversion.TheMaster-SlavearchitectureisshowninFigure2.

Figure2.TheMaster-SlaveArchitectureInageneticalgorithm,apopulationevolvesaccordingtoreproductionrules.Inthisworkwedescribetworeproduc-tionapproachesandcomparetheirimplementations.

3.2.ClassicalGenerationalApproach

Theclassicalgeneticalgorithm,knownasgenerationalbeginswithasetofpossiblesolutions-thepopulation.Thefitnessofeachindividualiscalculated,andbasedonthesevalues,acompletelynewpopulationisgeneratedviacrossover,mutationandselectionalgorithms.Weuseaninitialpopulationcomposedby21individualsandineachiterationanewpopulationisgeneratedandevaluated.Ourgenerationalalgorithmuseselitismtostorethebestindi-vidual,whichisnotreevaluatedbythefitnessfunction.Thecrossoverrateusedinthisgeneticalgorithmwasapproxi-mately60%andthemutationrateusedwas15%.

Aflowchartforaclassicalgenerationalgeneticalgo-rithmisshowninFigure3.

Figure3.Classicalgenerationalgeneticalgo-rithmflowchart

3.3.Steady-StateApproach

Inthetraditionalsteady-stateapproach,insteadofre-placingthewholepopulationonlyonesingleindividualisreplacedineachiteration.Theworstindividualofthepop-ulationissubstitutedbyanewindividual,generatedfrommutationand/orcrossoveroftheotherpopulationindividu-als.Asbefore,thepopulationiscomposedof21individu-als.

Inourwork,wedecidedtouseamodifiedsteady-statealgorithm,thatsubstitutethesevenworstindividualsofthepopulationandnotjusttheworstasinthetraditionalap-proach.Thisdecisionmakesabetteruseoftheavailablecomputationalresources,sincethecomputationaltestswereperformedinaclusterofeightmachines,decomposedasonemasterandsevenslaves.

Hence,inthisimplementation,the7worstindividualsarediscardedandsubstitutedbyindividualsgeneratedbyacrossoveronthe7betterindividuals.Thus,thecrossoverrateusedinourgeneticalgorithmisapproximately33%.Asbeforethemutationratewas15%.

Aflowchartforasteady-stategeneticalgorithmisshowninFigure4.

Figure4.Steady-stategeneticalgorithmflowchart

4.Implementation

ThegeneticalgorithmsweredevelopedusingtheCpro-gramminglanguage.TheMessagePassingInterfacelibrary[2]wasadoptedforthecommunicationbetweenmasterandslavemachines.Thesimulatorwasoriginallydevelopedasastand-aloneapplication.Tointegrateitwiththegeneticalgorithm,thesimulatorwasreestructuredasaCfunction,wheretheinputisthegeneticalgorithmchromosomeandtheoutputisthecalculatedfitness.

ThedevelopmentandexecutionofthetwoalgorithmsdescribedinthisworkwasdoneintheLaboratoryofCom-putationalPhysiology(FISIOCOMP),locatedinUniversi-dadeFederaldeJuizdeFora,inasmallclusterof8AMD64bitsmachines,with1Gbofmemory.Currently,theclus-terrunsaConectiva10OperatingSystem.Thecommuni-cationbetweenthemachinesisdonewitha1Gbps3Comswitch.Oneofthemachines,thefrontend,workedasthemastermachineandtheother7machinesworkedasslavemachines.

5.Results

Theclassicalgenerationalalgorithmwasexecutedfor25generationsandevaluated21individualsineachiteration.Hence,thisalgorithmevaluates525individualsduringitsexecution(25individuals×21iterations).Thisalgorithm

wasexecuted5timesandeachexecutiontookabout12hours,usingthe8availablemachines.Theaveragerela-tiveerror,i.e.thedifferencebetweentheestimated(α,β)andthecorrectparameters,wasaround30%.

Figure5showsthedifferencebetweentheobjectiveECGandtheECGestimatedbythegeneticalgorithm.

Figure5.OriginalandGAfittedECGsfortheclassicalalgorithm

Themodifiedsteady-statealgorithmwasalsoexecutedfor25iterations.Thefirstiterationevaluatedallthe21indi-vidualsofthepopulation.Thefollowingiterationsevaluatejustthe7newgeneratedindividuals.Hence,thisalgorithmevaluates189individualsduringitsexecution(21individu-als+7individuals×21iterations).

Thisalgorithmwasexecuted5timesandeachexecu-tiontookabout4.5hours,usingthe8availablemachines.Theaveragerelativeerroroftheestimatedparameterswasaround23%.

Afterimplementingandexecutingthegenerationalandthesteady-statealgorithmsandcomparingtheexecutiontimesandtheresultsfoundbythetwoalgorithmswecon-cludethatthesteady-stateapproachwasmoreefficientthanthegenerationalapproach.

Table1presentsacomparisonbetweentheexecutiontime(inhours),theaveragerelativeerroroftheestimatedparametersandnumberofevaluatedindividuals.GenerationalSteady-StateExecutionTime124.5EvaluatedIndividuals525189AverageError30%23%Table1showsthatthesteady-statealgorithmfoundbet-terpairofparameters(α,β)withshorterexecutiontimesandwithlessindividualevaluationsthanthegenerationalalgorithm.

6.Conclusions

Thisworkpresentedsomepreliminaryresultsrelatedtotheapplicationofcardiacmodelingontheestimationofcar-diactissuepropertiesfrominformationgivenbytransmuralelectrograms.Wecomparedtwoparallelimplementationsofgeneticalgorithm,hereusedasoptimizationmethods:theclassicalgenerationalmethodandthesteady-statebasedone.Motivatedbythestrongevidencesthatcardiacbulkconductivityvaluesanddistributionchangeundermanypathologicalconditions,inthisworkwehavemodeledthethecardiacconductivitychangesreportedforchronicin-farct.Theparallelgeneticalgorithmcoupledtothecar-diacbidomainmodelwasabletoestimate,fromgivenob-servedelectrograms,cardiactissueconductivityvalueswithanaccuracyof30%.Ourresultsshowthatthesteady-statealgorithmestimatedbetterpairsofparameters(α,β)withshorterexecutiontimesandwithlessindividualevaluationsthanthegenerationalalgorithm.Theexecutiontimeoftheparallelgeneticalgorithmtookaroundfourhoursina8-nodelinuxcluster.Furtherresearchisnecessaryinordertobettercharacterizethefeasibilityoftheinverseprocedureheredescribed.

7.Acknowledgements

WethankthefinancialsupportprovidedbyCNPq,project506795/2004-7andbyFAPEMIG,theMinasGeraisstatefoundationofBrazil,undertheprojectTEC-1358/05.

References

[1]M.Fallert,M.Mirotznik,S.Downing,E.Savage,K.Fos-ter,M.Josephson,andD.Bogen.Myocardialelectricalimpedancemappingofischemicsheepheartsandhealinganeurysms.Circulation,87:199–207,1993.

[2]M.P.I.Forum.Messagepassinginterface.availableat

http://www.mpi-forum.org/.2002.

[3]K.GimaandY.Rudy.Ioniccurrentbasisofelectrocardio-graphicwaveforms:Amodelstudy.CirculationResearch,90:889–896,2002.

[4]B.Hopenfeld,J.Stinstra,andR.Macleod.Mechanism

forstdepressionassociatedwithcontiguoussubendocar-dialischemia.JournalofCardiovascularElectrophysiology,15:200–1206,2004.

[5]S.Kostin,M.Rieger,S.Dammer,S.Hein,M.Richter,

W.Klovekorn,E.Bauer,andJ.Schaper.Gapjunctionremodelingandalteredconnexin43expressioninthefail-inghumanheart.MolecularandCellularBiochemistry,242:135144,2003.

[6]K.H.W.J.tenTusscher,D.Noble,andA.V.Noble,P.

J.andPanfilov.Amodelforhumanventriculartissue.J.Physiol,286:1573–1589,2004.

[7]R.WeberdosSantos,F.OtavianoCampos,L.Neu-mannCiuffo,A.Nygren,andW.Giles.Atx-iieffectsontheapparentlocationofmcellsinacomputationalhumanleftventricularwedge.JournalofCardiovascularElectro-physiology,17:S86–S95,2006.

[8]R.WeberdosSantos,G.Plank,S.Bauer,andE.Vigmond.

Parallelmultigridpreconditionerforthecardiacbidomainmodel.IEEETransactionsonBiomedicalEngineering,51:1960–1968,2004.

[9]D.Whitley.TheGENITORalgorithmandselectionpres-sure:Whyrank-basedallocationofreproductivetrialsisbest.InJ.D.Schaffer,editor,ProceedingsoftheThirdIn-ternationalConferenceonGeneticAlgorithms,SanMateo,CA,1989.MorganKaufman.

[10]G.Yan,W.Shimizu,andC.Antzelevitch.Characteristics

anddistributionofmcellsinarteriallyperfusedcanineleftventricularwedgepreparations.CirculationResearch,pages1921–1927,1998.