Stochastic non-linear oscillator models of EEG: the Alzheimer's

14 Pages·2015·0.88 MB·English

by  Hashem Ashrafiuon

View metadata, citation and similar papers at core.ac.uk brought to you by CORE
provided by Frontiers - Publisher Connector
ORIGINALRESEARCH
published:24April2015
doi:10.3389/fncom.2015.00048
Stochastic non-linear oscillator
models of EEG: the Alzheimer’s
disease case
ParhamGhorbanian1,SubramanianRamakrishnan2andHashemAshrafiuon1*
1DepartmentofMechanicalEngineering,CenterforNonlinearDynamicsandControl,VillanovaUniversity,Villanova,PA,
USA,2DepartmentofMechanicalandIndustrialEngineering,UniversityofMinnesotaDuluth,Duluth,MN,USA
In this article, the Electroencephalography (EEG) signal of the human brain is modeled
as the output of stochastic non-linear coupled oscillator networks. It is shown that
EEG signals recorded under different brain states in healthy as well as Alzheimer’s
disease (AD) patients may be understood as distinct, statistically significant realizations
of the model. EEG signals recorded during resting eyes-open (EO) and eyes-closed
(EC) resting conditions in a pilot study with AD patients and age-matched healthy
control subjects (CTL) are employed. An optimization scheme is then utilized to match
the output of the stochastic Duffing—van der Pol double oscillator network with EEG
signalsrecordedduringeachconditionforADandCTLsubjectsbyselectingthemodel
physical parameters and noise intensity. The selected signal characteristics are power
spectraldensitiesinmajorbrainfrequencybandsShannonandsampleentropies.These
Editedby: measures allow matching of linear time varying frequency content as well as non-linear
TobiasAlecioMattei,
signalinformationcontentandcomplexity.Themainfindingoftheworkisthatstatistically
Brain&SpineCenter-InvisionHealth-
KenmoreMercyHospital,USA significant unique models represent the EC and EO conditions for both CTL and AD
Reviewedby: subjects.However,itisalsoshownthattheinclusionofsampleentropyintheoptimization
RobinA.A.Ince,
process,tomatchthecomplexityoftheEEGsignal,enhancesthestochasticnon-linear
UniversityofManchester,UK
FanLiao, oscillatormodelperformance.
WashingtonUniversityinStlouis,USA
Keywords:EEG,Alzheimer’sdisease,stochasticdifferentialequations,duffing—vanderPol,entropy
*Correspondence:
HashemAshrafiuon,
DepartmentofMechanical
1. Introduction
Engineering,CenterforNonlinear
DynamicsandControl,Villanova
University,800E.LancasterAve., Quantitativeanalysisofhumanbrainelectroencephalography(EEG)recordingsaimedatenhanc-
Villanova,PA19085,USA ingourunderstandingofbraininjuriesanddisordersiscurrentlyanimportantresearcharea.In
hashem.ashrafi[email protected] additiontobeingusefulindiagnosis,suchanalysiscanprovideinsightsintotheunderlyingneuro-
physiologyoftheinjuryordisorder,therebyleadingtobettertreatmentandpreventivestrategies.
Received:28January2015 Alzheimer’s disease (AD) is the most common form of dementia and is the subject of intense
Accepted:07April2015 research. While no known cure exists, certain medications have shown promise in delaying the
Published:24April2015
symptoms(Dauwelsetal.,2010)promptingresearcherstoseekearlydiagnosisandintervention
Citation: strategies.Inthiscontext,analysisoftheEEGisapotentialnon-invasivetoolthatmayaidearly
GhorbanianP,RamakrishnanSand diagnosis of AD. However, the use of EEG signal analysis in order to improve the diagnosis of
AshrafiuonH(2015)Stochastic
ADisacomplexproblemwhere,despitesignificantadvances,anumberoffundamentalquestions
non-linearoscillatormodelsofEEG:
remainopen(Elgendietal.,2011).
theAlzheimer’sdiseasecase.
Front.Comput.Neurosci.9:48. ConsideringnowthecharacteristicsoftheEEG,sincethenon-stationarynatureofthesignal
doi:10.3389/fncom.2015.00048 is generally well-recognized (see, for instance Akin, 2002), decomposition using a fast Fourier
FrontiersinComputationalNeuroscience|www.frontiersin.org 1 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
transform (FFT) with sliding windows and the wavelet trans- and Moorman, 2000) as an improvement over approximate
formshavebeenthemostpopulartechniquesemployedtocap- entropy.
turethespectralpropertiesofEEG(DarvishiandAl-Ani,2007; In recent work, the authors proposed a phenomenological
Dauwels et al., 2010). However, linear transformation methods model of the EEG signal based on the dynamics of a stochas-
fail to address the non-linear characteristics of the EEG signal tic, coupled, Duffing- van der Pol oscillator network (Ghorba-
(Stam, 2005). Therefore, non-linear dynamic approaches have nianetal.,2015).Anoptimizationschemewasadoptedtomatch
beenattemptedaswell,mostlyinvolvingcomputationallycom- modeloutputwithactualEEGdataobtainedfromhealthysub-
plex time series analysis (Jeong, 2004). Several other aspects jectsinthetwodistinctrestingeyes-open(EO)andeyes-closed
of non-linear modeling and analysis in this context have also (EC) conditions and it was shown that the actual EEG signals
been studied in the literature (see, for instance, Stam, 2005 for in both cases were distinct realizations of the model with qual-
a review). These include frameworks based on a neural mass itativelydifferentnon-lineardynamiccharacteristics.Moreover,
model (Valdes et al., 1999; Huang et al., 2011), coupled oscilla- themodeloutputandtheactualEEGdatawereshowntobein
tors(Baieretal.,2005;Leistritzetal.,2007),continuummodels good agreement in terms of both the power spectra (frequency
(Kimetal.,2007),non-linearnon-stationarymodels(Celkaand content)andShannonentropy(informationcontent).
Colditz, 2002; Rankine et al., 2007), random neural networks In the present effort, we improve the model introduced in
(AcedoandMorano,2013),andchaoticphenomenaandstabil- Ghorbanianetal.(2015)bymatchingthesampleentropyofthe
ityaspects(Rodriguesetal.,2007;Dafilisetal.,2009).Stochastic modeloutputandEEGsignaltocaptureitscomplexity.Aglobal
approachesbasedonMarkovchainMonteCarlomethods(Het- optimization routine is employed in order to match the output
tiarachchi et al., 2012) and Markov process amplitude (Wang of with EEG recordings in terms of power spectrum, Shannon
et al., 2011) that take into account the inherent randomness of entropy, and sample entropy. The EEG signals were recorded
theEEGsignalhavealsobeenreported.Inthesamevein,limit underrestingECandEOconditionsinanearlierpilotstudyof
cycleoscillators(Hernandezetal.,1996;BurkeandPaor,2004) Alzheimer’sdisease(AD)patientsvs.age-matchedhealthycon-
aswellasstochasticsynchronization(BressloffandLai,2011)and trol(CTL)subjects(Ghorbanianetal.,2013).Themodelparam-
stochasticapproximation(Felletal.,2000;Sunetal.,2008)meth- eters obtained for the oscillators representing EC and EO EEG
odshavebeenconsideredinEEGmodeling.Notably,limitcycle signalsforCTLandADpatientsarecomparedinordertoestab-
behaviorateachofthebrainfrequencybandsappearstoprovide lishstatisticallysignificant,distinctmodelsforADandCTLsub-
amoreaccuraterepresentationoftheEEGsignalthanonebased jects under each condition. In addition, we present new results
onchaoticphenomena. fromaphasespacereconstructionanalysisofthemodeloutputto
Some of the most important features in non-linear dynamic matchtheactualEEGsignal.Theresultsindicatethattheanalyt-
and stochastic approaches are signal information content and icalmodeleffectivelycapturesthefrequencyspectrumandnon-
complexity as measured using various forms of information linearcharacteristicsoftheEEGsignalintermsofcomplexityand
entropy. Measures such as Shannon entropy (Shannon, 1948) informationcontent.Furthermore,itisshownthattheaddition
characterize the information content in a signal and higher ofsampleentropysignificantlyenhancesthemodelperformance
entropycorrespondstoincreasedrandomnessandchaoticbehav- intermsofcomplexityandnon-lineardynamiccharacteristics,as
ior (Abasolo et al., 2006). Importantly, one observes that, with demonstratedbyphasespacereconstruction.Theresultssuggest
respecttotheEEGsignal,higherinformationcontentcorrelates excitingnewpathwaystodevelopbettertoolsfordistinguishing
with better brain function (Shin et al., 2006). Furthermore, it pathological and normal brain states in AD and perhaps other
has been reported thatvariations in information entropic mea- neurologicaldiseasesanddisorders.
sures may be used to detect functional abnormalities in the The rest of the article is set as follows. Details of the EEG
brain caused by disorders or injuries (Slobounov et al., 2009). recordings, the analytical model, the optimization scheme and
Hence,informationcontentoftheEEGsignal,characterizedby thephasespacereconstructiontechniqueareprovidedinSection
information-entropicmeasures,maybeexpectedtobeimportant 2.TheresultsarepresentedinSection3anddiscussedinSection
in identifying distinct states of the brain. This is further rein- 4.Thearticlesconcludeswithcommentsonfurtherresearchin
forced by the recent results of McBride and colleagues on the Section5.
roleofinformationentropicandspectralanalysisinthestudyof
theearlystagesofAlzheimer’sdiseaseandmildTraumaticBrain 2. Materials and Methods
InjuryMcBrideetal.(2013a,b,2014).
Entropy may also be utilized to measure signal complex- 2.1.EEGRecordingBlocks
ity. For instance, embedding entropy provides information Twenty six AD patients and healthy age-matched CTL subjects
about how the EEG signal fluctuates in time by compar- were selected for this study (“A Brain-Computer Interface for
ing the time series with a delayed version of itself (Abasolo DiagnosingBrainFunction,”AspireIRB,HumanSubjectProto-
etal.,2006).Moreover,theconceptofapproximateentropywas colNumberPDMC-001,approvedonOctober7,2010).Ofthe
introduced as a measure of system complexity (Pincus, 1991) 26subjectsselected,onewithdrewandonedidnotqualifyasAD
and has been applied to brain wave signals (Quiroga et al., or CTL. Subjects were asked to relax and wear an EEG record-
2001).However,theapproximateentropymeasuresuffersfrom ingheadsetduringalternatingblocksofECandEOfollowedbya
drawbacks such as bias and inconsistency (Xu et al., 2010). varietyofcognitiveandauditorytasksandafinalEC-EOresting
Hence,thenotionofsampleentropywasintroduced(Richman period.
FrontiersinComputationalNeuroscience|www.frontiersin.org 2 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
FIGURE1|SchematicofthestochasticcoupledDuffing—vanderPoloscillators.
TABLE1|OptimalparametersoftheDuffing—vanderPoloscillatorfor subjects(10ECand10EO).Notethat,thesmallernumberofAD
ECandEOofCTLsubjects(N=40)andthep-valuesfromunpairedt-test, patientsalongwithsmallernumberofADpatientrecordingses-
Wilcoxonranksumtest,andBonferonnicorrection.
sions that were were not dominated by artifacts resulted in the
selectionofsmallerADsamplesize.
Parameter Eyes-Closed Eyes-Open t-test Wilcoxon Bonferroni
(EC) (EO)
2.2.EEGFeatures
k1 7286.5±192.4 2427.2±448.91 1e-15 1e-8 1e-7 The time-varying power spectrum in each of the major brain
k2 4523.5±282.3 499.92±84.04 1e-15 1e-8 1e-7 EEGfrequencybandswascalculatedusingshorttimefastFourier
b1 232.05±18.3 95.61±24.20 1e-15 1e-8 1e-7 transform(FFT)withslidingwindow,sinceagoodmodelmust
b2 10.78±2.3 103.36±9.22 1e-15 1e-8 1e-7 producesignalsthatcanmatchEEG’sfrequencycontent.Specif-
ǫ1 33.60±5.4 48.89±9.49 1e-15 1e-8 1e-7 ically,thepowerspectrumwascomputedinsevenranges:lower
ǫ2 0.97±0.19 28.75±1.74 1e-15 1e-8 1e-7 δ (1–2 Hz), upper δ (2–4Hz), θ (4–8Hz), α (8–13 Hz), lower
β (13–20Hz),upperβ (20–30Hz),andγ (30–60Hz).However,
µ 2.34±0.47 1.82±0.78 0.01 0.06 0.06
lowerδ andγ bandpowers,whichhappentohavelittlepower,
wereignoredduetounreliabilityofthedeviceinthosefrequency
TheEEGsignalswererecordedthroughasingle-dryelectrode ranges.
device at position Fp1 (based on a 10–20 electrode placement Shannon entropy was measured based on a sliding tempo-
system)withaBluetoothenabledtelemetricheadset.Thehead- ralwindowtechnique.Atemporalwindowwasdefinedtoslide
set’seffectivesamplerateis125Hz.Frequenciesbelow1Hzand alongthesignaltimerepresentationwithaslidingstep(interval
above 60 Hz (near Nyquist frequency) were filtered out by the or bin) to sample a part of the signal. A discrete entropy esti-
devicehardware.Oncomparisonof theEEG recordingsbythe matorwasapplied,inwhich10uniformintervalsequallydivided
devicewiththosefromotherwidelyaccepteddevices,frequencies therangeofthenormalizedobservedsignal.Thentheprobability
within2–30Hzweredeemedtobeveryaccurate. thatthesampledsignalbelongstotheintervalistheratiobetween
Therecordingdeviceeliminatedfrequentlyobservedartifacts the number of the samples found within each interval and the
includinglinenoise.Otherartifactsweremainlyduetoeye-and total number of samples of the signal. The Shannon entropy is
muscle-movements,whicharecommonatFp1locationandcan then calculated based on these probabilities (Shin et al., 2006),
be clearly identified by their high amplitudes compared to true separatelyforeach40-sEEGrecordingblock(5000samples).
EEGsignalrecordingsduringrestingstates.Theseartifactswere Sampleentropy(SE)isthenegativenaturallogarithmofthe
removedusingasimpleartifactdetectionthateliminatedanypart conditionalprobabilitythattwosequencesofatimeseries,simi-
of the signal greater than 4.5σ (standard deviation). The algo- larformpoints,remainsimilaratthenextpoint.ForgivenNdata
rithmalsoreconstructedthenulledsamplesusingFFTinterpo- points from a time series, [x(1),x(2),··· ,x(N)], we calculated
lationofthetrailingandsubsequentrecordeddata(Ghorbanian SE of each 40-s EEG recording block (5000 samples) by the
etal.,2013). statistic(Abasoloetal.,2006):
The EEG recordings in this study were obtained from sub-
Um+1(r)
jectsinanADpilotstudywith14control(CTL)subjectsand10 SE(m,r,N)= −ln , (1)
Alzheimer’sDisease(AD)patientspresentedinourearlierwork (cid:26) (cid:20) Um(r) (cid:21)(cid:27)
(Ghorbanian et al., 2013). Recording blocks of 40-s duration
wheremistherunlength,risthetolerancewindowsize,and
(approximately 5000 sample signals) from resting eyes-closed
(EC)andeyes-open(EO)conditionswereselected.Inall,60ran-
N−m
dom blocks were selected from the pilot study: 40 blocks from Um(r)= 1 U. (2)
i
controlCTLsubjects(20ECand20EO)and20blocksfromAD (N−m)(N−m−1)
Xi=1
FrontiersinComputationalNeuroscience|www.frontiersin.org 3 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
Intheaboveequation,U indicatesthenumberofk’s(1 ≤ k ≤ majorbrainfrequencyspectraandvanderPolnon-linearitypro-
i
N − m) such that the Euclidean distance between X (i) and videsself-excitedlimitcyclebehaviorwhichhavebeenpreviously
m
X (k), k 6= i, is less than or equal r and X (i) = [x(i),x(i + reportedforeachmajorbrainfrequencybands(BurkeandPaor,
m m
1),··· ,x(i+m−1)]. 2004).
Generally, large m or small r values result in number of WeconsideraphenomenologicalmodeloftheEEGbasedon
matches being too small for confident estimation of the condi- a coupled system of Duffing—van der Pol oscillators subject to
tionalprobabilityandviceversa(LakeandMoorman,2011).In whitenoiseexcitation,asshowninFigure1.Theequationsfor
thisstudy,weusedm = 2andr = 0.25σ basedontheconsis- themodelmaybewrittenas:
tencyoftheresultsandrecommendedrangesinpreviousstudies
(RichmanandMoorman,2000;Xuetal.,2010).
x¨ +(k +k )x −k x =−b x3−b (x −x )3
1 1 2 1 2 2 1 1 2 1 2
2.3.StochasticCoupledNon-linearOscillators  +ǫ1x˙1(1−x12), (3)
WerecallthattheEEGhasbeenmodeledintheliteraturetaking  x¨2−k2x1+k2x2 =b2(x1−x2)3
+ǫ x˙ (1−x2)+µdW,
intoaccountcharacteristicsincludingnon-linearity(bothchaotic 2 2 2
and non-chaotic), non-stationarity, and randomness of the sig- 
nal(Felletal.,2000;Rankineetal.,2007;Sunetal.,2008).The where x,x˙,x¨,i = 1,2 are positions, velocities, and accelera-
i i i
EEG has also been studied as the manifestation of underlying tionsofthetwooscillators,respectively.Parametersk,b,ǫ,i=
i i i
limitcycleoscillationsatagivenfrequencyandothersuchperi- 1,2 are, respectively, linear stiffness, cubic stiffness, and van
odic solutions (Hernandez et al., 1996; Burke and Paor, 2004). der Pol damping coefficient of the two oscillators. Parameters
While inspired by the above, the authors were fundamentally bs indicate the strength of the Duffing non-linearity resulting
i
motivatedtodevelopmodelsthatcanbetterreproducethesig- in multiple resonant frequencies while ǫs indicate the strength
i
nificant linear and non-linear characteristicsof actual EEG sig- of van der Pol non-linearity and determine the extent of self-
nals.Hence,weproposedaphenomenologicalmodeloftheEEG excitation and the shape of the resulting limit cycle. Parameter
based on a coupled system of Duffing—van der Pol oscillators µ represents the intensity of white noise and dW is a Wiener
subjected to white noise excitation (Ghorbanian et al., 2015). process (Gardiner, 1985; Higham, 2001) representing the addi-
ThisparticularoscillatorwasselectedbecausetheDuffingnon- tive noise in the stochastic differential system. The input exci-
linearity allows a system with only two oscillators capture the tation to the system is provided through µdW. The output
0.7
EEG Signal−−EC
0.6 Generated Signal−−EC
0.5
er 0.4
w
o
P 0.3
0.2
0.1
0
delta_L delta_U theta alpha beta_L beta_U gamma
Frequency (Hz)
0.4
EEG Signal−−EO
Generated Signal−−EO
0.3
er
w 0.2
o
P
0.1
0
delta_L delta_U theta alpha beta_L beta_U gamma
Frequency (Hz)
FIGURE2|ComparisonofmajorbrainfrequencybandmeanpowersofCTLEEGsignalsandoptimaloscillatormodeloutput;EC(top),EO(bottom).
FrontiersinComputationalNeuroscience|www.frontiersin.org 4 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
maybeselectedasanycombinationofthepositionsandveloc- where J is the objective function, p = [k ,k ,b ,b ,ǫ ,ǫ ,µ]
1 2 1 2 1 2
ities to mimic an EEG signal. Note that, the Euler-Maruyama thedecisionvariables,P andP thepowersinthemajorbrain
Ej Oj
method(Higham,2001)wasselectedtointegratethestochastic frequency bands for the normalized EEG signal and the model
differential equations in Equation (3) since standard numerical output,respectively,misnumberoffrequencybands(m = 7),
integrationmethodsarenotapplicable. S and S the Shannon entropies of the EEG signal and the
E O
modeloutput,respectively,SP andSP thesampleentropiesof
E O
the EEG signal and the model output, respectively, w and w
1 2
2.4.OptimizationFormulation
areweightingfactorforabsoluteShannonandsampleentropies,
Wehaveselectedthevelocityofthesecondoscillatorasthesys- respectively,and||representsabsolutevalue.Theweightingfac-
tem output approximating the EEG signal since it is directly torsw andw arerequiredtogiveequalimportancetopower
1 2
impacted by the noise. A global optimization search method spectrumandentropycharacteristicsofthesignal.Notethatthe
basedonamulti-startalgorithm(Ugrayetal.,2007)wasadopted magnitudeoftheoutputsignalsarematchedthroughnormaliza-
to determine the oscillator model parameters that can pro- tionofboththemodeloutputandtheEEGsignalwithrespectto
duce the output matching various EEG signals. The optimiza- theirstandarddeviations.
tionobjectivefunctionwasselectedastherootmeansquaredof The objective function minimization is subject to equality
the errors in power spectrum of each selected brain frequency constraintsrepresentedbythestate(Equation3)andinequality
bands plus weighted values of the errors in absolute Shannon constraintsrepresentedbythedecisionvariablelowerandupper
and sample entropies. Hence, the optimization goal is error bounds:
minimization:
0<k ≤1e4, 0<b ≤ 1k, 0<ǫ ≤ 1k,
i i 2 i i 3 i (5)
i=1,2, 0≤µ≤2.
m
minJ =v (PEj−POj)2+w1|SE−SO|+w2|SPE−SPO|, (4) Theconstraintsforbi’sandǫi’swereimposedtoavoidthechaotic
p uuXj=1 regime(Lietal.,2006)andprovideaperiodicstochasticresponse.
t
100
al
n
g
Si
G
E
E
0 10 20 30 40 50 60 70
el
od 100
M
d
Ol
− 10−2
al
n
g
Si
d 10−4
e
at
er
n
Ge 0 10 20 30 40 50 60 70
100
al
n
g
Si
d
e
at
er
n
e
G
0 10 20 30 40 50 60 70
Frequency (Hz)
FIGURE3|PowerspectrumofasampleCTLEC(top)EEGsignal;(middle)outputofstochasticoscillatormodelusingShannonentropy;(bottom)
outputofstochasticoscillatormodelusingShannonandsampleentropies.
FrontiersinComputationalNeuroscience|www.frontiersin.org 5 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
Noiseintensityisalsoconstrainedtoavoidaresponsedominated 2.6.PhaseSpaceReconstruction
by random noise. The initial guesses for the global optimiza- In addition to matching Shannon and sample entropies of the
tion search are randomly generated within the bounds defined modeloutputandEEGsignalthroughtheoptimizationprocess,
inEquation(5). it is of interest to investigate matching other features such as
The stochastic component was introduced as white noise, the phase plot which plays a significant role in non-linear time
which was generated through a normally distributed random seriesanalysis(KantzandSchreiber,2004).Itisknownthatany
variable and applied to the model via Wiener process. A new dynamicsystemcanbecompletelyrecoveredinthephasespace,
randomprocesswasgeneratedandappliedtothemodelduring which maybe reconstructed from the measured time domain
integrationoftheequations,ateachiterationoftheoptimization responseofthesystem(Nieetal.,2013).Whilephasespacecon-
algorithm. sistsofvelocityandpositionvariablesforamechanicalsystem,in
thecasewherejustthetimerepresentationofasignalisavailable,
aphasespacereconstructiontechniquebasedonthemethodof
2.5.StatisticalAnalysis
delaysisused(KantzandSchreiber,2004).
Akeyobjectiveofthephenomenologicalmodelinginthiswork
The main idea is that one does not need the derivatives to
is the ability to establish a correspondence between variations
form a coordinate system in which to capture the structure
in model parameters and the variations in the data obtained
of phase space, but instead one could directly use the lagged
from different physiological conditions. Hence, the parametric
variables:
unpaired t-test and non-parametric Wilcoxon rank sum statis-
tical testing methods were employed to determine the relative
significance of the model parameters. Furthermore, Bonferroni x(n+T)=x(t0+(n+T)1τs), (6)
correction was applied due to multiple comparisons problem
andadequacyofsamplesizesforstatisticaltestswereestablished where x(n) is the nth sample of the time series, 1τ the time
s
usingpoweranalysis. step, and T the delay integer to be determined. Then, a vector
100
al
n
g
Si
G
E
E
0 10 20 30 40 50 60 70
el 100
d
o
M
d
Ol
− 10−2
al
n
g
Si
ed 10−4
at
er
n
Ge 100 0 10 20 30 40 50 60 70
al
n
g
Si
d
e
at
er
n
e
G
0 10 20 30 40 50 60 70
Frequency (Hz)
FIGURE4|PowerspectrumofasampleCTLEO(top)EEGsignal;(middle)outputofstochasticoscillatormodelusingShannonentropy;(bottom)
outputofstochasticoscillatormodelusingShannonandsampleentropies.
FrontiersinComputationalNeuroscience|www.frontiersin.org 6 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
of(embedding)dimensiondmaybeconstructedusingthetime ameasurementatonetimefromameasurementtakenatanother
lagsas: time. Consider the time series nth sample x(n) and its value
aftertimedelayTwiththeassociatedprobabilitydistributionsof
[x(n),x(n+T),x(n+2T),··· ,x(n+(d−1)T)]. (7) P(x(n))andP(x(n+T)),respectively.Theaverageinformation
whichcanbeobtainedaboutx(n+T)fromx(n)isgivenbythe
Time-delay embedding is probably one of the best systematic mutualinformationofthetwomeasurements(Abarbaneletal.,
methods for converting scalar data to multidimensional phase 1993;Mizrach,1996):
space (Abarbanel et al., 1993; Burke and Paor, 2004; Nie et al.,
2013).Anappropriateandsuccessfulreconstructiondependson P(x(n),x(n+T))
I(x(n),x(n+T))=log , (8)
thechoiceofbothtimedelayTandtheembeddingdimensiond 2(cid:20)P(x(n))P(x(n+T))(cid:21)
(Nieetal.,2013).
Inthisstudy,theappropriatevalueofthetimelagwasdeter- where P(x(n),x(n + T)) is the joint probability of the mea-
minedusingtheaveragemutualinformationmethodappliedto surements x(n) and x(n+T) calculated using a binning-based
eachEEGrecordingblock.Theideabehindmutualinformation method,inwhich20uniformintervalsdividedtherangeofthe
istoidentifytheamountofinformationthatcanbelearnedabout measurementsequally.Theaveragemutualinformationbetween
0.49
n
o
ati 0.485
m
o
nfr 0.48
al I
u
ut
M 0.475
e
g
a
er 0.47
v
A
0.465
0 5 10 15 20 25
0.605
n
o
ati 0.6
m
o
nfr 0.595
ual I 0.59
ut
M
e 0.585
g
a
er 0.58
v
A
0 5 10 15 20 25
0.595
n
o
ati 0.59
m
o
nfr 0.585
ual I 0.58
ut
M
e 0.575
g
a
er 0.57
v
A
0 5 10 15 20 25
lag
FIGURE5|AveragemutualinformationforasampleCTLEC(top)EEGsignal;(middle)outputofstochasticoscillatormodelusingShannonentropy;
(bottom)outputofstochasticoscillatormodelusingShannonandsampleentropies.
FrontiersinComputationalNeuroscience|www.frontiersin.org 7 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
measurementsofanyvaluex(n)andx(n+T)istheaverageover After specifying the correct time delay T, an appropriate
allpossiblemeasurementsofI(x(n),x(n+T))(Abarbaneletal., embedding dimension, d, should also be found for the phase
1993): space reconstruction. If d is too small, the trajectories will not
beunique.Ontheotherhand,toolargead willresultinaddi-
tional computational cost by requiring extra dimensions (Nie
I(T)= P(x(n),x(n+T))I(x(n),x(n+T)). (9)
etal.,2013).
x(n)X,x(n+T)
IfTistoosmall,themeasurementsx(n)andx(n+T)willhavetoo 3. Results
muchoverlap.However,ifTistoolarge,thenI(T)willapproach
zeroandnothingrelatesx(n)tox(n+T).Itissuggestedthatthe Theoptimizationalgorithmwasseparatelyappliedtodetermine
proper T can be chosen as the first minimum of I(T) which is the model parameters (decision variables) for each of the 60
notnecessarilyoptimalbuthasbeenshowntoworkwell(Abar- selectedEEGsignalsusingtheweightingfactorsw =w =0.35.
1 2
baneletal.,1993;Nieetal.,2013).Ifinacase,nominimaexists These weighting factors give equal importance to the entropy
forI(T),thechoiceofT =1or2hasbeensuggested(Abarbanel measuresandpowerspectrum.Wethencategorizedtheresult-
etal.,1993). ing60setofmodelparametersintofourgroupsbasedrecording
1
0.5
)
T
+ 0
n
(
x
−0.5
−1
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
x(n)
1
0.5
)
T
+ 0
n
(
x
−0.5
−1
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
x(n)
1
0.5
)
T
+ 0
n
(
x
−0.5
−1
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
x(n)
FIGURE6|ReconstructedphaseplotofasampleCTLEC(top)EEGsignal;(middle)outputofstochasticoscillatormodelusingShannonentropy;
(bottom)outputofstochasticoscillatormodelusingShannonandsampleentropies.
FrontiersinComputationalNeuroscience|www.frontiersin.org 8 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
conditionsandsubjectdiagnosis:EC-CTL,EO-CTL,EC-AD,and amoreflatfrequencydistributionfromupperδ tolowerβ fre-
EO-AD. quencybands.Furthermore,ShannonandSEvaluesoftheEEG
signals and the model outputs for the EC and EO cases show
3.1.HealthyEyes-ClosedandEyes-OpenResults closeagreement.Shannonentropyvalueswere1.80±0.08and
Initially,westudiedthemodelsderivedfortheECandEOEEG 1.92 ± 0.08 for EC EEG and model output, respectively, and
signalsofCTLsubjectsforvalidationpurposes.Themeansand 1.71 ± 0.11 and 1.57 ± 0.15 for EO EEG and model output,
standarddeviationsoftheoptimalvaluesofthemodelparame- respectively.While,SEvalueswere1.04±0.20and1.17±0.22
tersforEC-CTLandEO-CTLEEGsignalsarelistedinTable1. forECEEGandmodeloutput,respectively,and0.97±0.20and
The p-values from the two statistical tests and non-parametric 1.20±0.18forEOEEGandmodeloutput,respectively.These
methodafterBonferronicorrectionsindicatethatthedifferences resultsshowasignificantimprovementoverourpreviousmodel
betweenofallparametersofthetwomodelsarestronglystatisti- whereonlyShannonentropywasused(Ghorbanianetal.,2015).
callysignificantwiththeexceptionofnoiseintensity.Notethat, Theimprovementisclearlyobservedinthethepowerspectraof
µisalsofoundstatisticallysignificantusingt-testbutisslightly sampleECandEOEEGsignalsandtheircorrespondingoptimal
offwhenthenon-parametricmethodisused. modeloutputs,respectivelyshowninFigures3,4.Bothfigures
In order to ensure that adequate sample sizes are used, the demonstratemoredistributedspectraofthemodeloutputswith
minimum required difference between means of two groups of similarnoisecomplexitiestotheactualEEGsignalswhenSEis
dataforeachparameterarecomputed.Asexpectedduetovery addedtotheobjectivefunction;i.e.,powerspectraofthesignals
small p-values, the sample size for statistical testing is found to withoutmatchingofSEhaveverydiscretepeaksunliketheEEG.
be sufficient with more than 99.9% power for all parameters TheimpactofSEtomatchsignalcomplexityisfurtherdemon-
exceptnoiseintensityµ,whichwasnotfoundtobestatistically stratedthroughphaseplotreconstructionofthetimeseries.Aver-
significantusingthenon-parametricmethod. agemutualinformationforasampleECEEGsignalandoutputs
Power spectrums of the optimal stochastic oscillator model oftheoptimalstochasticoscillatormodelsareshowninFigure5
outputandEEGsignalsfortheECandEOcasesofCTLsubjects asafunctionoflagtime.ThefirstminimumoccursatT = 5lag
arepresentedinFigure2whereθ,α,andβ bandpowersshow samplesforboththeEEGsignalandtheoptimalmodelderived
excellentagreement.Thecomparisonrevealedthat,asexpected, with both Shannon and sample entropies while T = 3 for the
theoptimalmodeliscloselyfollowingtheα-banddominancein output of the model derived solely based on Shannon entropy.
theECcases.While,intheEOcases,theoptimalmodelfollows TheresultingreconstructedphaseplotsoftheECEEGsignaland
theoutputsofthetwooptimalmodelsarepresentedinFigure6.
Clearly,thereconstructedphaseplotsoftheEEGandtheoutput
ofthemodelderivedusingbothShannonandsampleentropies,
TABLE2|OptimalparametersoftheDuffing—vanderPoloscillator displaysimilarbehavior.Whiletheoutputofthemodelderived
modelforECandEOofADsubjects(N=20)andthep-valuesfrom using only Shannon entropy is qualitatively different form the
unpairedt-test,Wilcoxonranksumtest,andBonferonnicorrection. EEGsignalintermsofcomplexityandnoise.Indeedthisresult
providesfurtheraffirmationthatthestochasticDuffing—vander
Parameter Eyes-Closed Eyes-Open t-test WilcoxonBonferroni
PolmodelyieldsanoutputthatmatchestheactualEEGdatain
(EC) (EO)
termsofnon-linearcharacteristicsobservedinthephasespace.
k1 1742.1±197.91 3139.9±1040.9 0.0005 0.0025 0.009
3.2.Alzheimer’sDiseasevs.ControlResults
k2 1270.8±277.13 650.32±175.76 1e-5 0.0005 0.002
Next, we studied the models derived for the EC and EO EEG
b1 771.99±126.81 101.1±27.86 1e-12 0.0001 0.001
signalsofADsubjects.Themeanandstandarddeviationofthe
b2 1.91±0.22 81.3±9.76 1e-15 0.0001 0.001
optimalvaluesofthemodelparametersforEC-ADandEO-AD
ǫ1 63.7±11.64 56.3±5.75 0.0884 0.021 0.063
EEG signals are listed in Table2 along with the p-values from
ǫ2 20.7±5.64 19.12±2.87 0.4234 0.879 0.95
thetwostatisticaltestsandthenon-parametrictestafterBonfer-
µ 1.78±0.8 1.74±0.67 0.905 0.879 0.95
ronicorrectionsindicatingthatthedifferencesbetweenonlythe
TABLE3|Thep-valuesfromunpairedt-test,Wilcoxonranksumtest,andBonferonnicorrectionforcomparisonofmodelparametersbetweenAD
(N=20)andCTL(N=40)subjects.
Parameter t-test(EC) Wilcoxon(EC) Bonf.(EC) t-test(EO) Wilcoxon(EO) Bonf.(EO)
k1 1e-30 1e-5 4e-5 0.013 0.027 0.08
k2 1e-23 1e-5 3e-5 0.0034 0.0015 0.007
b1 1e-17 1e-5 5e-5 0.58 0.027 0.08
b2 1e-12 1e-5 6e-5 1e-6 1e-5 9e-5
ǫ1 1e-10 6e-5 1e-5 0.031 0.0018 0.007
ǫ2 1e-15 5e-5 7e-6 4e-12 1e-5 7e-5
µ 0.02 0.06 0.06 0.80 0.70 0.7
FrontiersinComputationalNeuroscience|www.frontiersin.org 9 April2015|Volume9|Article48
Ghorbanianetal. Stochasticnon-linearoscillatormodelsofEEG
first four parameters of the two models are statistically signifi- the model parameters of CTL and AD subjects under EO con-
cant.Next,weseparatelycomparedthemodelparametersofEC dition are not, however, as strong, though they are still mostly
andEOEEGsignalsofCTLsubjectswiththoseADpatients. statisticallysignificant.IntheEOcase,parameterµisnotstatis-
Table3liststhep-valuesfromthetwostatisticaltestingmeth- ticallysignificantusingeithermethodandt-testdoesnotfindb
1
odsandthenon-parametricmethodafterBonferronicorrections tobestatisticallysignificanteither.
comparingCTLvs.ADsubjectsunderseparateECandEOcon- The power analysis results for 90%, 95%, 99%, and 99.9%
ditions.Theresultsindicatethatthedifferencebetweenallmodel fortwostatisticalarelistedinTables4,5forECandEOcases,
parameters of CTL and AD subjects under EC condition are respectively. The actual difference between means are given
stronglystatisticallysignificantexceptfornoiseintensity.Again, within parentheses following each parameter. The results indi-
µisalsofoundstatisticallysignificantusingt-testbutisslightly cated that our sample size for statistical testing in EC case
offwhennon-parametricmethodisused.Thedifferencebetween between AD and CTL subjects was sufficient for all parameters
TABLE4|Minimumrequireddifferencebetweenmodelparametermean TABLE5|Minimumrequireddifferencebetweenmodelparametermean
valuesofECADvs.ECCTLforvariousdesiredpowersofstatisticaltests. valuesofEOADvs.EOCTLforvariousdesiredpowersofstatisticaltests.
Parameter 90% 95% 99% 99.9% Parameter 90% 95% 99% 99.9%
1k1(5544.3) 281.39 313.53 371.72 439.46 1k1(712.78) 1009.12 1124.36 1333.04 1575.97
1k2(3252.7) 406.65 453.09 537.18 635.08 1k2(175.81) 175.81 195.89 232.25 274.57
1b1(539.94) 106.44 118.60 140.61 166.24 1b1(5.49) 36.86 41.07 48.69 57.56
1b2(8.87) 2.82 3.14 3.72 4.40 1b2(22.02) 13.61 15.17 17.98 21.26
1e1(30.09) 11.54 12.86 15.25 18.03 1e1(7.41) 12.27 13.68 16.22 19.17
1e2(19.79) 4.64 5.17 6.13 7.25 1e2(9.62) 3.14 3.50 4.15 4.91
1µ(0.56) 0.87 0.97 1.15 1.36 1µ(0.07) 1.08 1.21 1.43 1.70
0.5
EEG Signal−−AD−−EC
Generated Signal−−AD−−EC
0.4
0.3
er
w
o
P
0.2
0.1
0
delta_L delta_U theta alpha beta_L beta_U gamma
Frequency (Hz)
0.5
EEG Signal−−AD−−EO
Generated Signal−−AD−−EO
0.4
0.3
er
w
o
P
0.2
0.1
0
delta_L delta_U theta alpha beta_L beta_U gamma
Frequency (Hz)
FIGURE7|ComparisonofmajorbrainfrequencybandmeanpowersofADEEGsignalsandoptimaloscillatormodeloutput;EC(top),EO(bottom).
FrontiersinComputationalNeuroscience|www.frontiersin.org 10 April2015|Volume9|Article48