Assessment of Autonomic Control and Respiratory Sinus Arrhythmia Using Point Process Models ...

13 Pages·2009·0.66 MB·English

Assessment of Autonomic Control and Respiratory Sinus
Arrhythmia Using Point Process Models of Human Heart
Beat Dynamics
The MIT Faculty has made this article openly available. Please share
how this access benefits you. Your story matters.
Citation Zhe Chen, E.N. Brown, and R. Barbieri. “Assessment of
Autonomic Control and Respiratory Sinus Arrhythmia Using Point
Process Models of Human Heart Beat Dynamics.” Biomedical
Engineering, IEEE Transactions on 56.7 (2009): 1791-1802.
Print. © 2009 Institute of Electrical and Electronics Engineers
As Published http://dx.doi.org/10.1109/tbme.2009.2016349
Publisher Institute of Electrical and Electronics Engineers
Version Final published version
Accessed Thu Apr 11 13:56:06 EDT 2019
Citable Link http://hdl.handle.net/1721.1/51871
Terms of Use Article is made available in accordance with the publisher's policy
and may be subject to US copyright law. Please refer to the
publisher's site for terms of use.
Detailed Terms
IEEETRANSACTIONSONELECTROMAGNETICCOMPATIBILITY,VOL.56,NO.7,JULY2009 1791
Assessment of Autonomic Control and Respiratory
Sinus Arrhythmia Using Point Process Models
of Human Heart Beat Dynamics
ZheChen∗,Member,IEEE,EmeryN.Brown,Fellow,IEEE,andRiccardoBarbieri,SeniorMember,IEEE
Abstract—Tracking the autonomic control and respiratory si- lated by the autonomic nervous system [1]. It has long been
nus arrhythmia (RSA) from electrocardiogram and respiratory understoodthatthehealthyheartisinfluencedbymultipleneu-
measurementsisanimportantproblemincardiovascularcontrol.
ralandhormonalinputsthatresultinvariationsofdurationinthe
Weproposeapointprocessadaptivefilteralgorithmbasedonan
interbeatintervals(R-Rintervals).Thesynergicinteractionbe-
inverse Gaussian model to track heart beat intervals that incor-
poratesrespiratorymeasurementsasacovariateandprovidesan tweenthetwobranchesoftheautonomicnervoussystemtothe
analyticformforcomputingadynamicestimateofRSAgain.We hearthasamajorinfluenceinregulatingthecardiacdynamics
useKolmogorov–Smirnovtestsandautocorrelationfunctionanal- and physiological mechanism of HRV. In particular, parasym-
yses to assess model goodness-of-fit. We illustrate the properties
pathetic influences decrease the firing rate of pacemaker cells
ofthenewdynamicestimateofRSAintheanalysisofsimulated
in the heart’s sinus-atrial (SA) node, whereas sympathetic in-
heart beat data and actual heart beat data recorded from sub-
jectsinafour-stateposturalstudyofheartbeatdynamics:control, fluences have the opposite effect [9]. In cardiovascular physi-
sympathetic blockade, parasympathetic blockade, and combined ology,itisknownthatlungvolumetendstobecorrelatedwith
sympatheticandparasympatheticblockade.Inadditiontogiving variations in the timing of heart beat or HRV. Typically, HR
an accurate description of the heart beat data, our adaptive fil-
slowsdownduringexpirationandspeedsupduringinspiration.
ter algorithm confirms established findings pointing at a vagally
Thisphenomenonisknownastherespiratorysinusarrhythmia
mediated RSA and provides a new dynamic RSA estimate that
canbeusedtotrackcardiovascularcontrolbetweenandwithina (RSA)[28].RSAisprimarilymediatedbymodulationofvagal
broadrangeofpostural,pharmacological,andageconditions.Our outflow to the SA node. Quantification of RSA provides im-
paradigmsuggestsapossibleframeworkfordesigningadevicefor portantinformationaboutsomeofthemechanismsinvolvedin
ambulatory monitoring and assessment of autonomic control in
cardiorespiratorycoupling[20],[28].Inclinicalpractice,RSAis
bothlaboratoryresearchandclinicalpractice.
oftentreatedasanindirectandnoninvasivemeasureofparasym-
Index Terms—Adaptive filters, autoregressive (AR) processes, patheticcardiaccontrol[20],[24],eveninthepresenceofpaced
heart rate variability (HRV), point processes, respiratory sinus
breathing[29],andmayalsobeconsideredasareliableindica-
arrhythmia(RSA).
torofcardiacdysfunction[26].Aquitecomprehensivereview
I. INTRODUCTION onRSAmaybefoundin[10].
Acentralgoalinbiomedicalengineeringappliedtocardiovas-
HEARTrate(HR)andHRvariability(HRV)areimportant
cularcontrolistodevelopquantitativemeasuresandinformative
quantitative markers of cardiovascular control, as regu-
indexesthatcanbeextractedfromphysiologicalmeasurements.
Specifically,amajorchallengeincardiovascularengineeringis
Manuscript received September 4, 2008; revised December 12, 2008 and
February 5, 2009. First published March 4, 2009; current version published todevelopstatisticalmodelsandapplysignalprocessingtools
June12,2009.ThispaperwaspresentedinpartattheProceedingoftheIn- toinvestigatevariouscardiovascular–respiratoryfunctions[8],
ternationalConferenceonAcoustics,Speech,andSignalProcessing(ICASSP)
such as HRV, RSA, and baroreflex. In the literature, numer-
2008,LasVegas,NV.ThisworkwassupportedbytheNationalInstitutesof
Health(NIH)underGrantR01-HL084502,GrantR01-DA015644,GrantDP1- ous methods have been proposed for quantitative HRV analy-
OD003646, Grant GM-26691, and Grant AG-9550, and by the Division of sis[1],[25],includingpointprocessanalysis[2],[3],frequency-
ResearchunderGrantRR-79.Asteriskindicatescorrespondingauthor.
∗Z. Chen is with the Neuroscience Statistics Research Laboratory, Mas- domainanalysis[7],andnonlineardynamicsanalysis[19].RSA
sachusetts General Hospital, Harvard Medical School, Boston, MA 02114 was defined in [22] as “the difference between the maximum
USA, and also with the Department of Brain and Cognitive Sciences, Mas- HRrateaftertheonsetofinspiratoryflowandtheimmediately
sachusetts Institute of Technology, Cambridge, MA 02139 USA (e-mail:
minimumHR”;whereasRSAwascalculatedin[35]usingthe
[email protected]).
E.N.BrowniswiththeNeuroscienceStatisticsResearchLaboratory,Mas- formula:100×(meanlongestR-R-meanshortestR-R)/mean
sachusetts General Hospital, Harvard Medical School, Boston, MA 02114 R-Rinterval.Saulandcolleagues[27]proposedatransferfunc-
USA, and also with the Harvard–MIT Division of Health Science and
tionanalysisapproach forevaluating theRSA,whichrequires
Technology, and the Department of Brain and Cognitive Sciences, Mas-
sachusetts Institute of Technology, Cambridge, MA 02139 USA (e-mail: directmodelingoftheSAnode.Abivariateautoregressive(AR)
[email protected]). model was proposed in [4] and [5] to evaluate a time-varying
R.BarbieriiswiththeNeuroscienceStatisticsResearchLaboratory,Mas-
index of RSA (within a temporal window). However, none of
sachusettsGeneralHospital,HarvardMedicalSchool,Boston,MA02114USA,
andalsowithMassachusettsInstituteofTechnology,Cambridge,MA02139 these RSA indexes provide atruly instantaneous evaluation of
USA(e-mail:[email protected]). thecardiorespiratorydynamics.
Colorversionsofoneormoreofthefiguresinthispaperareavailableonline
SeveralissuesinRSAassessmentfromR-waveeventsareyet
athttp://ieeexplore.org.
DigitalObjectIdentifier10.1109/TBME.2009.2016349 tobeaddressed[10].First,estimatesofRSAhavebeenderived
0018-9294/$25.00©2009IEEE
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
1792 IEEETRANSACTIONSONELECTROMAGNETICCOMPATIBILITY,VOL.56,NO.7,JULY2009
TABLEI
COMPARISONOF4TWO-PARAMETERPROBABILITYMODELSFORHEARTBEATR-RINTERVALANDHR
fromeitherHRorheartperioddata.Theformerismorecom- gain across different posture and pharmacological conditions.
monly computed in clinical practice, whereas the latter would Finally,discussionsandconclusionaregiveninSectionIV.
bepreferredonbiometricgrounds,especiallywhentheinterest
is in indexing parasympathetic control because of the relative II. METHODS
linearitybetweenvagalfrequencyandheartperiod.Second,the
Inthissection,wepresenttheheartbeatintervalandtheHR
R-Rintervalseriesareunevenlyspacedintime.Directapplica-
probabilitymodels,theinstantaneousestimatesofHRandHRV
tiontothesedatatospectralanalysisisnotappropriate,andis
fromtheheartbeatinterval(R-R)modelparameters,theexten-
usuallysolvedbyuseofinterpolatingfilters.Inaddition,longer
sion to a bivariate model with RP as covariate to derive RSA
heartperiodsmaysignificantlydecreasetheNyquistfrequency
measures,thepointprocessadaptivefilteringalgorithmforin-
under fast respiratory oscillation, giving rise to possible alias-
stantaneousassessmentoftheHRandRSAindexes,andfinally
ingeffects.Third,standardtime-seriesanalysisusuallyassumes
goodness-of-fit tests to evaluate how well these estimates de-
thatthedatashowatleastweaksensestationarity,thusrequir-
scribethestochasticstructureofthewaveeventsextractedfrom
ing particular care in choosing appropriate data segments for
an ECG. A detailed description of the conceptual framework
analysis,orrequiringremovalofnonstationarytrends.
ofthehistory-dependentpointprocessmodelandtheadaptive
To address these issues, we investigate different probability
paradigmcanbefoundinourpreviouspublications[2],[3].
models for human heart beat intervals with an adaptive point
process filtering paradigm [3] and illustrate the analysis with
A. PointProcessProbabilisticModelsofHeart
bothsyntheticdata,andECGandlungvolumerecordingsfrom
BeatInterval
apreviousstudy[31]–[33]underanautonomicblockadeassess-
ment protocol. Furthermore, we extend the inverse Gaussian The R-wave events mark the electrical impulses from the
probability model to take into account the influence of respi- heartsconductionsystemthatrepresentventricularcontractions.
ration (RP) on HRV, based on which, we derive an analytic Hence,theyareasequenceofdiscreteoccurrencesincontinuous
formforcomputingthefrequencyresponseofRSA.Modeling time, and as such, form a point process. Suppose that we are
accuracy is evaluated via goodness-of-fit tests, including the givenasetofdiscreteR-waveevents{u }J detectedfromthe
j j=1
Kolmogorov–Smirnov(KS)test.Toourknowledge,thisisthe ECG,letRRj =uj −uj−1 >0denotethejthR-Rinterval,or
first effort in the literature that investigates the dynamic RSA equivalently, the waiting time until the next R-wave event. By
effect within a point process framework. In line with the most treating the R-waves as discrete events, we propose different
conventional guidelines [1], [10], our paradigm resolves the parametric point process probability models (Table I) in the
oldconflictbetweenheartperiodandHR,obviatestheneedof continuous-timedomain.
interpolationwiththepotentialtosolvepossiblealiasingprob- As an example, assuming history dependence, the waiting
lems, allows for instantaneous measures at virtually any time timet−u untilthenextR-waveeventmaybemodeledasan
j
resolution offering (in contrast to the interpolated R-R inter- inverseGaussianmodelasfollows[18]:
valvalues)amorerigorousfrequencyanalysis,andovercomes (cid:1) (cid:2) (cid:1) (cid:2)
θ 1/2 θ(t−u −µ )2
nonstationarityissuesassociatedwithwindow-basedestimation p(t)= exp − j t , (t>u )
2πt3 2µ2(t−u ) j
models. t j
(1)
Thepaperisorganizedasfollows.SectionIIfirstintroduces
where u denotes the previous R-wave event occurred before
the point process framework, then presents several probability j
time t, θ >0 denotes the shape parameter, and µ denotes the
modelsoftheheartbeatdynamics,andfinallyproposestheex- t
instantaneous R-R mean value. Because the parasympathetic
tendedbivariateprobabilitymodelbyinclusionofRP(e.g.,lung
andsympatheticinputstotheSAnodecanoccuronamillisec-
volume)measurementsforthepurposeofquantifyinginstanta-
ondtimescale,buttheireffectscanlastforseveralseconds,the
neousRSA.InSectionIII,bothsyntheticandrealexperimental
intervalsmustbemodeledasdependentontherecenthistoryof
recordings areusedtoillustrateandvalidate theinstantaneous
theSAnodeinputs
RSAgainascomputedbythenovelpointprocessalgorithm.In
addition,statisticaltestsareconductedontheautonomicblock- (cid:3)p
ade protocol to evaluate intersubject statistical trends of RSA µt ≡µRR(t)=a0 + aiRRt−i. (2)
i=1
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
CHENetal.:ASSESSMENTOFAUTONOMICCONTROLANDRESPIRATORYSINUSARRHYTHMIA 1793
Namely,themeanvalueismodeledbyaunivariatep-orderAR Notethat,ingeneral,themeasurementsofRR(beatpercycle)
process,whichisassumed(approximately)tobeinfluencedby and RP have different sampling frequencies. In practice, our
thepastpR-Rintervalvalues.Inordertotrackthenonstationary framework allows for two ways to tackle this issue. The first
behaviorofheartbeatdynamics,allparametersinvolvedin(1) approach is to resample the RP measurements to synchronize
and (2) are adaptive; hence, the instantaneous mean µ (t) with the heart beat and obtain the RP values at the beat time;
RR
is time-varying, which is determined by the time-varying AR thesecondoneistotreatthemasseparatemeasurements(with
coefficients {a (t)}p . On the other hand, the instantaneous different sampling rates) and conduct frequency analysis with
i i=0
varianceoftheinverseGaussianmodelcanbederivedas[18] extra caution. The second approach could be useful to avoid
aliasing effects in the case of low Nyquist frequency in the
µ3 (t)
σ2 (t)= RR . (3) presenceoflongheartbeats.
RR θ(t)
Intermsoffrequencyanalysis,theRSAeffectisreflectedby
In a similar fashion, we can derive the mean and variance of the fact that the R-R interval has a spectral component mod-
R-Rintervalforotherprobabilitymodels,suchastheGaussian, ulated by the respiratory variable [27]. With a linear system
lognormal,andgammamodels(seeTableIforabriefsummary). assumption, RSA can be estimated with transfer function or
frequency response analysis using standard signal processing
tools[5],[27].
B. InstantaneousIndexesofHRandHRV
Given the parametric AR model (7), we can evaluate the
HR is defined as the reciprocal of the R-R interval. For t
frequencyresponsefortheR-Rintervalitself
measured in seconds, a new variable r =c(t−u )−1 (where (cid:4)
j (cid:4)
cof=-va6r0iasb/mleisnf)ocramnublae,dtehfienHeRdipnrobbeaabtsilpiteyrpm(irn)u=te.pB(yc(tthe−chuaj)n−g1e)- H1(ω)= 1−(cid:7)pi=11ai(k)z−i(cid:4)(cid:4)z=ej2πfs (8)
isgivenby
aswellasthefrequencyresponsefortheRSA
(cid:4) (cid:4) (cid:4)
(cid:4)dt(cid:4) (cid:7) (cid:4)
p(r)=(cid:4) (cid:4)p(t) (4) q b (k)z−j(cid:4)
dr H12(ω)= (cid:7)j=1 j z=(cid:4)(cid:4)ej2πfs (9)
and the mean and the standard deviation (std) of HR r can 1− p a (k)z−i(cid:4)
i=1 i
bederived(seeTableI).Essentially,theinstantaneousindexes z=ej2πfs
of HR and HRV are characterized by the mean µ and std where fs is the beat rate of the RR and RP are sampled at the
HR
σ , respectively. In the case of inverse Gaussian model, we same frequency as the beat.1 With the estimated time-varying
HR ARcoefficients {a (k)}and{b (k)}attimet=k∆,wemay
have[2],[3] i j
evaluate the dynamic power spectrum (parametric autospec-
µHR =µ˜−1 +θ˜−1 (5) trum)orthegain(amplitude)inthefrequencydomain[7]
(cid:5) (cid:6) (cid:4) (cid:4)
σHR = 2µ˜µ˜θ˜+2θ˜ 1/2 (6) PRR(ω,t)=(cid:4)σRR(t)(cid:4)(cid:4)H1(cid:4)(ω,t)(cid:4)(cid:4) (10)
(cid:4) (cid:4)
whereµ˜ =c−1µ andθ˜=c−1θ. RSAgain(ω,t)=(cid:4)H12(ω,t)(cid:4). (11)
RR
Sincetwomajorrhythmsincardiovascularvariabilityanalysis
C. BivariateProbabilisticModelforRSAAssessment aretheoneoccurringatthefrequencyoftheMayerwaves(LF,
0.04–0.15 Hz) and the one triggered by RP (HF, 0.15–0.5 Hz,
In the probability models considered thus far, we have only
±0.04Hzaroundtherespiratoryrate)[1],wecancomputethe
usedtheR-Rintervaltimeseriestoestimatetheinstantaneous
power or the gain across these frequencies over time for both
meanµ .PhysiologysuggeststhatHRisinfluencedbyother
RR
(10)and(11).
physiological covariates, such as changes in lung volume due
Hence,from(10)and(11),wecanestimatetherelevant(in-
torespiratoryactivity[4],[5].Thisfactfurthermotivatesusto
stantaneousandmean)statistics,suchastheLFpower,theHF
incorporate RP as a covariate into the model. Specifically, for
power,andtheRSAgaininHF.Wemayalsocomputethedy-
theinverseGaussianmodel,wemayreplacetheinstantaneous
namicLF/HFpowerratio.Asmall(orlarge)LF/HFratioindi-
meanin(1)by
catesrelativelypredominantvagal(orsympathetic)control[1].
(cid:3)p (cid:3)q
µt ≡µRR(t)=a0 + aiRRt−i + bjRPt−j (7) D. AdaptivePointProcessFiltering
i=1 j=1
Inpractice,wecanbinacontinuous-timepointprocesswitha
whereRPt−j denotesthepreviousjthRPmeasurementbefore certainbinsize∆.Thebinsizehastobesmallenoughtonotonly
timet.Now,theinstantaneousmeanµ isdescribedbyabi-
RR containoneeventatmostinsideeachbin,butalsotocharacterize
variateAR-typemodel.Itshouldalsobenotedthattheextended
bivariate model shall not be limited to the inverse Gaussian
distribution;itisalsostraightforwardtoextendittotheGaus- 1Here,weassumethatRRandRPmeasurementshavethesamesampling
ratein(7).Inpractice,theRPtimeseriestypicallyhasahighersamplingrate,
sianorlognormaldistribution,dependingonwhicheverismore
butitcanalwaysberesampledandinterpolatedtoobtaindatapointsatthetime
desirable. ofheartbeats.
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
1794 IEEETRANSACTIONSONELECTROMAGNETICCOMPATIBILITY,VOL.56,NO.7,JULY2009
the dynamics at the timescale of interest. It is known from cdffromtheuniformdensitydefinedas(j−0.5)/J.Ideally,if
pointprocesstheory[2],[3],[14]thattheconditionalintensity themodeliscorrect,thepointsshouldlieonthe45◦line,andthe
function(CIF)λ(t)isrelatedtotheintereventprobabilityp(t) 95% confidence interval lines are y =x±(1.36/(J −1)1/2).
withaone-to-onerelationship TheKSdistance,definedasthemaximumdistancebetweenthe
KSplotandthe45◦line,isusedtomeasurelack-of-fitbetween
p(t)
λ(t)= (cid:8) . (12) themodelandthedata.
1− t p(τ)dτ
uj ofItnheadgdist:ioAnC, Fw(ema)ls=o(c1o/mJp−utemt)h(cid:7)e aJu−tmocgorgrelatio.nIffuthnectgions
TheestimatedCIFcanbeusedtoevaluatethegoodness-of-fitof j j=1 j j+m j
are independent, they are also uncorrelated; hence, ACF(m)
theprobabilitymodelfortheheartbeatdynamics.Thequantity
λ(t)∆approximatelyyieldstheprobabilityofobservingabeat shallbesmall(around0andwithinthe95%confidenceinterval
1.96/(J −1)1/2)forallvaluesofm.
duringthe[t,t+∆)interval.
Letξdenoteavectorthatcontainsallofunknownparameters
inanyparametricprobabilitymodel(inthecaseofinverseGaus- III. RESULTS
sianmodel,ξ =[{ai}pi=0,{bj}qj=1,θ]T),wecanrecursivelyes- In this section, we conduct a probabilistic analysis of heart
timatethemviaadaptivepointprocessfiltering[3] beatdatawiththestochasticpointprocessparadigm.Themajor
advantage of casting the heart beat interval within the point
ξ =ξ (13)
k|k−1 k−1|k−1 processframeworkistoallowforthepossibilitytomodeland
Pk|k−1 =Pk−1|k−1 +W (14) evaluate the instantaneous statisticsof HR, HRV, and RSA, at
arbitraryfine-timeresolution.
ξk|k =ξk|k−1 +Pk|k−1(∇logλk)[nk −λk∆] (15) Westartbyvalidatingthemodelwithsyntheticdata,andil-
(cid:5) (cid:6)
Pk|k = Pk−|k1−1 +∇λk∇λTk λ∆ −∇2logλk[nk −λk∆] −1 laulgstorraittihnmg taochoimghplaigrihstotnhewfiathst-etrsaticmkiantegsafbriolimtyaofwthinedpoowin-btapsreod-
k
cess filter. After describing experimental protocol and model
(16)
initialization,weconfirmthebetterperformanceoftheinverse
where P and W denote the parameter and noise covariance GaussianmodelwithRPascovariateasthemodelwithoverall
matrices, respectively; ∆=0.005 s denotes the time bin size; best fit. We then focus on instantaneous RSA estimation from
and∇λk =∂λk/∂ξk and∇2λk =∂2λk/∂ξk∂ξTk denotesthe the real data, and illustrate several examples where the point
first-andsecond-orderpartialderivativesoftheCIFwithrespect process filter estimates novel dynamic signatures of RSA. Fi-
toξattimet=k∆,respectively.Theindicatorvariablenk =1 nally, we further validate our indexes by presenting a concise
ifaheartbeatoccursintime((k−1)∆,k∆]and0otherwise. group study reporting significant differences of RSA with age
Remark: To avoid numerical problems that might occur in andchangesinpostureordrug.
the matrix inverse and to increase the numerical stability, one
canreplacetheoriginalposteriorcovarianceupdate(16)witha A. SimulatedRSAStudyWithSyntheticData
Fisher’sscoringstep[30]
(cid:5) (cid:6) Toillustrateandvalidateourproposedalgorithminestimating
Pk|k = Pk−|k1−1 +∇λk∇λTk λ∆k −1. (17) tHheerRe,StAhe,wReRfiarnsdtcRoPndtuimcteasseirmieuslaatrieognesntuedraytewditbhysaynctlhoesteidc-dloaotap.
bivariateAR(8)model[5].AllRPmeasurementsareassumed
In addition, if the prediction covariance is badly conditioned,
tobesynchronizedatthetimeofheartbeats.Thesimulatedtime
thepredictioncovarianceisretainedastheposteriorcovariance.
serieslastabout23min.
The instantaneous R-R interval is comodulated by its own
E. Goodness-of-FitTests
historyaswellastheRPseries
Model goodness-of-fit is assessed based upon the time-
rescaling theorem [14]. Given a point process specified by (cid:3)8 (cid:3)8
J discrete even(cid:8)ts: 0<u1 <···<uJ <T define the random RRk =a0 + aiRRk−i + bjRPk−j +ek (18)
variablesz = uj λ(τ)dτ forj =1,2,...,J −1.Then,the i=1 j=1
j uj−1
random variables z s are independent, unit-mean exponen- where e denotes a Gaussian random noise component with
j k
tiallydistributed.Byintroducingthevariableoftransformation zero mean and std 20. The synthetic RR and RP time series
v =1−exp(−z ),v sareindependent,uniformlydistributed weregeneratedinawaythattheRRhasoneLFat0.1Hzand
j j j
withintheregion[0,1].Letg =Φ−1(v )(whereΦ(·)denotes one HF at 0.2 Hz (they are both fixed along time). The RSA
j j
thecumulativedensityfunction(cdf)ofthestandardGaussian gain achieves the maximum gain 1 around 0.2 Hz for the first
distribution), then g s will be independent standard Gaussian 10 min, and then increases by a factor of 2 in the remaining
j
random variables. The KS test is used to compare the cdf of time.Spectralanalysisalsoconfirmedthatthesetwotimeseries
v againstthatoftherandomvariablesuniformlydistributedin havethemaximumcoherenceat0.2Hz.
j
[0,1].TheKSstatisticisthemaximumdeviationoftheempirical We used the inverse Gaussian distribution as the probabil-
cdffromtheuniformcdf.Tocomputeit,thev saresortedfrom ity model for the heart beat interval, with its instantaneous
j
thesmallesttothelargestvalue,andplottedagainstvaluesofthe meanmodeledby(7).TheinitialbivariateARcoefficientsare
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
CHENetal.:ASSESSMENTOFAUTONOMICCONTROLANDRESPIRATORYSINUSARRHYTHMIA 1795
Fig.2. Comparisonbetweentheadaptivepointprocessfilter(darkblacktrace)
andRLSfilter(graytrace)intheestimatedRSAcurve(averagedovermultiple
independentMonteCarloruns).Toppanel:thecompletetrace.Bottompanel:
thezoom-intracefrom450to750s.TheRSAunitismillisecondsperliter.
TABLEII
F(iings.ta1n.tanEesotuimsaHteRdVt)i,maencdoRurSsAesgoafiµnRstRat,isσtRicRs,foµrHoRne(iMnsotanntetaCneaorulosrHeRal)i,zσatHioRn MOANVTEERACGAIRNLGOORVSEAR4G0AIINNDSETPAETNISDTEINCST(RMEEAALNIZ±ATSIOTNDS,)TCIMOEMAPANRDISTORNIAL
ofthesimulatedheartbeatdata.Inthefirstpanel,thesuperimposedcurveisthe COMPUTEDFROMADAPTIVEPOINTPROCESSANDRLSFILTERS
originalR-Rintervaltimeseries. INSIMULATIONEXPERIMENTS
estimated from the first 100 s of the R-R and RP measure-
ments.ThenoisecovariancematrixW in(14)wassetasW =
diag{3×10−3,1×10−7,...,1×10−6,...,1×10−2}.Fig.1
shows the estimated time courses of several statistical indexes
to a lack of dynamics across time and a slow tracking perfor-
fromonesimulatedheartbeatRRandRPtimeseries.
mancetoreachthesteadystateoftheRSAgain.Ofnote,ifwe
Asexemplarycomparison,weimplementedarecursiveleast-
passasmoothingwindowonourestimates,wecouldgetsimi-
squares(RLS)filter[21]usingamovingwindow-basedmethod
lartrendsasintheRLScase,wheresuchwindowisimplicitly
(basedonRRandRPtimeseries)toestimatethetime-varying
includedintheestimationprocess.Inotherwords,RLSmaybe
bivariateARcoefficients.TheforgettingfactoroftheRLSfilter
considered a window-smoothed, unevenly sampled version of
waschosentobe0.98inordertowellbalancethebias-variance
thepointprocessfilter.
tradeoff [12]. Specifically, we conducted 40 Monte Carlo ex-
Insummary,inlightofmanyexperimentalobservations(in-
perimentsforbothadaptivepointprocessandRLSfilters,and
cludingadditionalsimulationsnotreportedhere),itappearsthat
compared their averaged RSA estimates (see the top panel of
notonlythepointprocessfilterproducesabettercharacteriza-
Fig. 2). Table II summarizes the comparison of Monte Carlo
tion in the instantaneous estimates of HR and HRV indexes
meanandstdstatisticsoftheseindexes inthefirstandsecond
(see[3]),butitsbivariateextensioncansimultaneouslyprovide
halfofthedata.Themeanvaluesofthepointprocessalgorithm
anaccurateestimateoftheRSAgain(seeMonteCarlocompar-
seem to be more accurate, with an averaged variance that is
ison). More importantly, we have illustrated how the adaptive
higherbyafactorof3.ThelowervariancevaluesfortheRLS
pointprocessfilteriscapableoftrackingsharpdramaticchanges
aremainlyattributabletothewindowingeffect.Acarefulexam-
in RSA gain. As a further validation for the model fit, the KS
ination of the estimates reveals that the adaptive point process
plot and the autocorrelation plot were generated and shown in
filter has a faster tracking performance (to approach the true
Fig. 3. As seen from the figure, the KS fit is within the 95%
RSAgainvalue,i.e.,from1to2ms/Linoursimulationexam-
confidenceinterval,indicatingthatthemodelusedhereisquite
ple) than the RLS filter (see the bottom panel of Fig. 2). Our
satisfactory;intheautocorrelationplot,nearlyallofpointsare
algorithmreachesthe95%lowerthresholdtowardtheexpected
alsolocatedinsidethe95%confidencebounds.
valueat615s(15saftertheabruptchangeingain),whereasthe
RLS estimate reaches the same threshold only at 720 s (120 s
B. ExperimentalProtocol:AutonomicBlockade
afterthechange).Thisisnotverysurprising,sincetheRLSfil-
terusesanintrinsicexponentialmovingwindow-basedmethod Parasympathetic tone usually dominates in healthy, resting
andupdatesvaluesonlyatthebeats,inamuchgreatertimescale individuals. When a resting subject is given atropine (ATR,
thanthepointprocessfilter[6].Thesecharacteristicsalsoyield a muscarinic receptor antagonist that blocks parasympathetic
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
1796 IEEETRANSACTIONSONELECTROMAGNETICCOMPATIBILITY,VOL.56,NO.7,JULY2009
10–16 belong to the ATR group and subjects 20–27 belong to
the PROP group. In addition, subjects 11, 12, 14, 16, 21, 22,
24,26,and27belongtotheoldagegroupandsubjects10,13,
15,20,23,and25belongtotheyoungagegroup.Fig.5shows
a snapshot of R-R intervals and respiratory recordings from
a representative subject with varying conditions (in this case,
PROPisadministeredfirst).Asexpected,notethatthegeneral
RRdecrease,bothincontrolandPROP,isaccompaniedbyan
increaseinventilationwhenthesubjectisstanding(observedin
onlysomeofthesubjects).AsATRisadministered,theintrinsic
Fig.3. KSplotandtheautocorrelationplotforthesimulatedheartbeatdata. HR shows the expected acceleration with markedly reduced
dynamics, which are not affected by change in posture. Of
note,insomesubjects,weobservedmildincreasesinHRwith
changeofpostureduringDB,thismaybeduetowaningofthe
firstadministereddrugeffects:inthesecases,evenaverymild
autonomic modulation may overcome any mechanical effect
ontheSAnode[11].Goingbacktoourrepresentativesubject,
notethehighnonstationaritiesinthecontrolPROPepoch,with
importantbradycardicevents,supposedlyduetoexclusivevagal
modulation. This example well summarizes some of the main
issues involved in HRV studies that motivated our research:
1)presenceoffastdynamics;2)presenceofhighnonstationari-
tieseveninsupposedlysteady-stateconditions;and3)dynamic
Fig.4. Diagramoftheautonomicblockadeprotocol. changes in ventilation to be accounted for accurate RSA
assessment.
effects), HR usually increases substantially. In contrast, if the
subject is given propranolol (PROP, a receptor antagonist that
C. ModelandParameterInitialization
blockssympatheticeffects),HRusuallydecreasesonlyslightly.
If both sympathetic and parasympathetic effects are blocked, The order of the AR model was determined based on the
HRiscalledtheinstrinsicHR[9]. Akaike information criterion (AIC) (by prefitting a subset of
The experimental data considered in our study were previ- the data) as well as the KS distance in the post hoc analysis.
ously presented in [31]–[33]. A schematic diagram is shown In all univariate AR cases, the order p=8 was chosen from
in Fig. 4 to describe different stages of the protocol. In each {2,4,6,8,10}. In bivariate AR analyses, the order p=q =8
epoch, 5 min segments of continuous ECG and lung volume was used. The initialAR coefficients are estimated by solving
were recorded. In the drug administered state, either ATR the Yule–Walker equation using about 60–80 s of the initial
(0.04 mg/kg intravenous (iv) over 5 min, parasympathetic recordings[4].
blockade) or PROP (0.2 mg/kg iv over 5 min, sympathetic
blockade)wasdeliveredtothesubject.Inthedoubleblockade D. ModelComparison
(DB), the inputs from both sympathetic and parasympathetic
Inordertochooseareasonablesubsetofparametricmodels,
branches of the autonomic nervous system were suppressed.
weperformedapreliminaryhistogramanalysisandprobability
SubjectswererandomizedtofirstreceiveeitherATRorPROP,
fitfortherecordedR-Rintervaltimeseries(seeFig.6resultfor
and then the alternate drug; after 10 min, all measurements
therepresentativesubjectinFig.5).We,consequently,selected
were repeated. All segments were recorded in a steady-state
fourcharacteristicprobabilitystructures:gamma,Gaussian,log-
condition (or as stationary as possible). R-wave spike trains
normal,andinverseGaussian.
weredetectedfrom360HzECGrecordings.Thelungvolume
Thecomparativeresultsofthefourprobabilitymodelscon-
dataweredigitallyrecordedat360Hz,measuringthecalibrated
sidered for the heart beat data are presented in Table III. The
outputscorrespondingtoribcageandabdominalcompartment
inverseGaussianmodelachievestheoverallbestfitintermsof
volumechangesassociatedwithRP.AfterverifyingtheNyquist
smaller KS distance, especially during the control and PROP
frequency condition for the R-R intervals, the long volume
epochs, in both supine and upright positions. The lognormal
recordingswereresampledatthebeattimes.
modelachievesbetterperformanceduringtheDBepochs.The
A total of 17 healthy volunteers participated in the study.
gammamodelhastheworstperformanceamongthefourprob-
Fifteensubjects (sixyoung subjects and nine oldsubjects) are
ability models tested here. All models perform rather unsatis-
included in the present study. Here, for convenience, we have
factorilyduringtheATRepochs.
renamedthesubjectswithnumbers.2Morespecifically,subjects
When the bivariate (inverse Gaussian) model was further
applied to the heart beat time series together with the lung
2Comparingwiththeoriginalstudy[31]:10=NO,11=AS,12=CA,13=
volume recordings, it was found that the inclusion of the RP
DB,14=GB,15=GJ,16=RB,20=TM,21=EO,22=GH,23=HU,
24=JW,25=KN,26=ME,and27=RN. covariateimprovestheKSfit(secondlastcolumn,TableIII)in
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
CHENetal.:ASSESSMENTOFAUTONOMICCONTROLANDRESPIRATORYSINUSARRHYTHMIA 1797
Fig.5. SnapshotofR-Rintervals(inmilliseconds)andlungvolumeRPmeasures(RP,calibratedandzeromean)undersixdifferentconditions(subject20).Top
sixpanels:supineposture.Bottomsixpanels:uprightposture.NotethatallRR(aswellasallRP)plotsarevisualizedwithinthesamescales
Fig.6. HistogramanalysisandprobabilityfitforthecontrolandDBconditionsinsupineposition(subject20;seeFig.2).Intheprobabilityfitplots,ifthedata
fitthetestedprobabilitydistribution,thedatapointswillmatchthestraightdashedline.
themajorityofthesubjectsunderallthreepharmacologicalcon- (especially in the control cases), and adding covariates makes
ditions(ATR,PROP,andDB).Fig.7illustratesacomparative theonlineadaptationmorechallenging.
exampleusing(2)and(7)inthe“upright+PROP”condition. Asseeninourexperiments,thebivariatemodelhasthebest
In a few subjects that were examined, it was found that the performance in control and PROP conditions, followed by the
inclusion of RP measurements does not improve the KS fit. DBcondition,andhastheworstfitintheATRcondition.Specif-
ThismightbeduetothefactthattheR-Rintervalshistoryitself ically,forthegoodness-of-fitinthecontrolcondition,about85%
might be already sufficient to characterize the instantaneous ofthetime-rescaledpointswereinsidethe95%confidenceinter-
mean statisticfor the inverse Gaussian probability distribution valoftheKSplot,whichimpliesthatroughly85%ofthebeats
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
1798 IEEETRANSACTIONSONELECTROMAGNETICCOMPATIBILITY,VOL.56,NO.7,JULY2009
TABLEIII
PERFORMANCECOMPARISONBETWEENFOURUNIVARIATEPROBABILITYMODELS
groupandayoungsubjectfromthePROPgroup,respectively.
Sincetheplotsarescaledtoevidencethedynamics,meanvalues
arereportedinthefiguretoallowforamoredirectcomparison
among epochs. The first observation is that whenever ATR is
administered,togetherwiththesharpmeanRSAdecreasedueto
theabsenceofparasympatheticmodulation,theRSAdynamics
almostdisappeartothepointthatquantizationeffectsduetothe
bin resolution become evident. On the other hand, interesting
trendscanalsobeobservedinthesetwosubjectsinthepresence
ofparasympatheticmodulation.Inparticular,weobservethesat-
Fig.7. ComparisonofinverseGaussianmodelswiththemeanastheunivari- urationeffectsinCUfortheATRsubjectandinPROP-upright
ate(toprow)andbivariate(bottomrow)ARmodelsforsubject20(upright,
forthePROPsubject.NotealsothefasteroscillationsintheCS
PROP).Leftpanel:estimatedtime-varyingprobabilitydensityfunctionofthe
instantaneous µRR(t). Middle panel: KS plot. Right panel: autocorrelation epochfortheyoungsubject(Fig.10).Again,thesetrendscould
function.(Dashedlinesindicatethe95%confidencebounds.) notbeobservedwithstationaryorwindow-basedanalysis.
arecorrectlypredictedbythemodel(lastcolumnofTableIII).
F. StatisticalAnalysisinGroupStudy
It is also important to pinpoint that if the same criterion were
to be applied to any other window-based model, there would Upon computing the RSA gain statistics (Table IV), we
notbeasinglecasewhereitwillbesatisfied.Morespecifically, investigate if there is any statistical difference between pos-
theKScurveswouldfallfarawayfromthe45◦ diagonal,with ture/pharmacological/age conditions,despitethefactthatonly
consequentsignificantlygreaterKSdistancevalues(see[2]). alimitednumberofsubjectsareavailableinthepresentstudy.
Specifically, we applied a nonparametric Mann–Whitney test
E. InstantaneousRSAAnalysis (alsoknownastheranksumtest)tocomparetwoindependent
samples.Inthecaseofmultiplecomparison,wealsoadjustedthe
Ourpointprocessbivariateframeworkallowedforestimation
p-value according to the Bonferroni correction. We found that
of instantaneous measures of RSA. For the protocol analysis,
within the young age group, the supine versus upright posture
an updating delta interval of 5 ms was chosen. Fig. 8 shows
effect (p<0.01) as well as the control versus DB drug effect
dataandresultsfromsubject25forthreeepochs:controlsupine
(p<0.01)bothshowstatisticalsignificance.Thesetwoeffects
(CS),controlupright(CU),andPROPsupine(PS).Inaddition
failtoreachsignificancewithintheoldagegroup,possiblydue
to the expected mean RSA decrease with change in posture,
tohighintersubjectvariability.Inthecomparisonbetweentwo
there is an expected mean RSA increase with PROP adminis-
agegroupsintheCScondition,wealsofoundastatisticallysig-
tration (when compared to the respective control epoch). The
nificant difference (p<0.01). Our observations are consistent
new estimates show novel interesting dynamics. In particular,
withpreviousfindingspointingataweakerautonomiccontrol
notethesharpincreasesinRSAaround160,190,and230sin
withincreasingage[31],[32].
CS,themarkedoscillatorytrendinCU,andthesharpdecreases
(around 105, 160, 185, and 230 s) from a saturated maximum
IV. DISCUSSIONANDCONCLUSION
RSA level at 440 ms/L in PS. All these dynamic signatures
could not have possibly been inferred by looking at the RR We have presented an extended point process paradigm for
andRPtimeseries,orbyapplyingstationaryorwindow-based humanheartbeatintervalswithRPasacovariatetoassessauto-
analyses.Theseobservationspointatourestimatesaseffective nomiccontrolasquantifiedbyRSA.Ourmethodisvalidatedby
instantaneous measures for a unique novel characterization of MonteCarlosimulationsusingsyntheticdata,andthenewfast-
vagalmodulationatsmalltimescalelevels. trackinginstantaneousRSAassessmentisillustratedasapplied
Figs. 9 and 10 plot the instantaneous RSA mean gain (HF: totheexperimentalautonomicblockaderecordings.Theseex-
0.15–0.5Hz)forallrecordingsfromanoldsubjectfromtheATR amplesrevealnovelinterestingdynamictrendsreflectiveofthe
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.
CHENetal.:ASSESSMENTOFAUTONOMICCONTROLANDRESPIRATORYSINUSARRHYTHMIA 1799
Fig.8. R-Rinterval,lungvolumeRPmeasure(adjustedtozeromean),andestimateddynamicRSAmeangain(0.15–0.5Hz)inthreeconsecutiveepochs
(subject25).FromCStoCU,RRdecreasesandRPincreasessignficantly,andRSAdecreases.FromCUtoPROPupright,RRincreasesandRPdecreases,and
RSAincreasessignificantly.
Fig.9. EstimatedinstantaneousRSAgain(unit:millisecondsperliter)intheHF(0.15–0.5Hz)range(subject14fromtheold/ATRgroup).Thenumberineach
subplotindicatesthemeanvalueoftheRSAgainaveragedovertheentirerecording(whichiscomputedusingtheunnormalizedRPmeasure),allRSAunitsare
inmillisecondsperliter.Inthiscase,thefollowingRSAmeangainrelationshipholds:supine>upright(exceptDB),control>ATR,andcontrol>DB.
nonstationary nature of cardiovascular control, whereas sim- A. ChoiceofProbabilityModel
ple summary statistics confirm established findings related to
Inmodelingtheheartbeatintervalduringthecontrolepochs,
RSA measures across different posture/pharmacological/age
the inverse Gaussian model achieves the best performance,
conditions. Several points described next are worth further
which is in agreement with our earlier claims [2], [3]. The
discussion.
Gaussianmodelachievesasimilarperformancesince,whenthe
Authorized licensed use limited to: MIT Libraries. Downloaded on November 23, 2009 at 17:04 from IEEE Xplore. Restrictions apply.