extended smoother Kalman Joint and state parameter estimation of the hemodynamic model byiterative Biomedical Signal Processing and Control

(1)

ContentslistsavailableatScienceDirect

Biomedical Signal Processing and Control

j ou rn a l h o m e p a g e :w w w . e l s e v i e r . c o m / l o c a t e / b s p c

Joint parameter and state estimation of the hemodynamic model by iterative extended Kalman smoother

Serdar Aslan

^a,∗

, Ali Taylan Cemgil

^b

, Murat S¸ amil Aslan

^c

, Behc¸ et U˘gur Töreyin

^d

, Ata Akın

^e,∗∗

aInstituteofBiomedicalEngineering,Bo˘gazic¸iUniversity,Kandilli, ˙Istanbul,Turkey

bDepartmentofComputerEngineering,Bo˘gazic¸iUniversity, ˙Istanbul,Turkey

cTübitakBilgem ˙IltarenAdvancedTechnologiesResearchInstitute,Ankara,Turkey

dC¸ankayaUniversity,FacultyofEngineering,DepartmentofElectricalandElectronicsEngineering,Ankara,Turkey

eDepartmentofMedicalEngineering,AcıbademUniversity,Istanbul,Turkey

a r t i c l e i n f o

Articlehistory:

Received30March2015

Receivedinrevisedform7September2015 Accepted9September2015

Availableonline27September2015

Keywords:

Hemodynamicmodel

ExtentedKalmanﬁlter/smoother CubatureKalmanﬁlter/smoother

a b s t r a c t

Thejointestimationoftheparametersandthestatesofthehemodynamicmodelfromthebloodoxygen leveldependent(BOLD)signalisachallengingproblem.Inthefunctionalmagneticresonanceimaging (fMRI)literature,quiteinterestingly,manyproposedalgorithmsworkonlyasafilteringmethod.This makestheestimationofhiddenstatesandparameterslessreliablecomparedwiththealgorithmsthat usesmoothing.Instandardimplementations,smoothingisperformedonlyonce.However,jointstate andparameterestimationcanbeimprovedsubstantiallybyiteratingsmoothingschemessuchasthe extendedKalmansmoother(IEKS).InthefMRIliterature,extendedKalmanfilteringisthoughttobeless accuratethanstandardparticlefiltering(PF).WecomparedEKFwithPFandobservedthatthecontraryis true.WeimprovedtheEKFperformancebyaddingsmoother.Byiterativeschemejointhemodynamicand parameterestimationisimprovedsubstantially.WecomparedIEKSperformancewiththesquare-root cubatureKalmansmoother(SCKS)algorithm.Weshowthatitsaccuracyforthestateandtheparameter estimationisbetterandmuchfasterthaniterativeSCKS.SCKSwasfoundtobeabetterestimatorthanthe dynamicexpectationmaximization(DEM),EKF,locallinearizationfilter(LLF)andPFmethods.Weshow inthispaperthatIEKSisabetterestimatorthaniterativeSCKSunderdifferentprocessandmeasurement noiseconditions.Asaresult,IEKSseemstobethebestmethodweevaluatedinallaspects.

1. Introduction

Inestimationtheory,therearetwomaincategoriesforthestate estimation.Thesearethefilteringandsmoothing.Inthefilteringof thestateestimate,theobservationsequencey₁,y₂,···,y_kory_1:kis usedtofindthestatisticalinformationaboutthestatesx₁,x₂,···, x_k.Weareinterestedtofindingtheconditionalprobabilitydensity function(pdf)p(x_k|y1:k)forthefiltering.Sincebythenewobser- vationy_k+1wecanfindtheconditionalpdfp(x_k+1|y1:k+1),filtering is alsoknownas online estimation, whereas, inthe smoothing approach,alltheobservationsequenceincludingthefuturedata (y_1:N,N>k)areusedtofindp(x_k|y1:N)[1].Sinceforeverytimestep weneedthecompleteobservationsequence,thiskindofprocessing

∗ Correspondingauthor.Tel.:+905545694249.

∗∗ Correspondingauthor.Tel.:+902123117683.

E-mailaddresses:serdar.aslan@gmail.com(S.Aslan),ata.akin@acibadem.edu.tr (A.Akın).

isalsoknownasofﬂineestimation. Theuseoffuture datacon- tainsinformationaboutthepastdataand,asaresult,improvesthe stateandparameterestimatesofthesystem.Thechoiceofthestate estimationtechniquealsohasanimportanteffectontheparame- terestimation.Havingamoreaccuratestateestimationalgorithm alsoimprovestheparameterestimationofthesystem.InthefMRI modelinversionliterature,mostofthestateestimationtechniques workonlyintheﬁlteringsense.However,inmosthemodynamic dataanalysiscases,thecomputationisnotperformedinrealtime.

Byusingsmoothingtechniques,theestimationofboththestateand theparametersofthehemodynamicsystemsisimproved.Trade- offisthatincorporatingsmoothingintoexistingalgorithmsleads toanincreasedcomputationtime.

In this paper, we worked on thejoint parameter and state estimationof thehemodynamicnonlineardynamic systemrep- resentation.The bloodoxygenlevel dependent(BOLD)signalis a measure for the localized hemodynamic response, which is observed after neuronal activation [2]. Thisresponse is a non- linearfunctionof theblood ﬂowand theblood oxygencontent http://dx.doi.org/10.1016/j.bspc.2015.09.006

(2)

[3].Hemodynamicresponsetypicallylastsforseveralsecondsto respondanddecay[4].Thisresponsechangesfromsubjecttosub- ject,evenindifferenttimesduringthedayanddifferentregions ofthebrain[4,5].Thegeneralshapeoftheresponseisatypically skewedbell-shapedfunctionand iscommonlyrepresentedasa mixtureofgammafunctions[6].Alternatively,thetotalsystemrep- resentingthehemodynamicresponsecanbemodeledbynonlinear differentialequations[7–9].

WeestimatetheparametersandthestatesbyusingtheBOLD signalandthefunctionalrepresentationofthesystem.Inthispaper, weusethefollowingfunctionalmodel:

x_k+1=f(,x_k)+w_k (1)

y_k=h(,x_k)+vk (2)

Here,fandharenonlinearfunctionsofthestatex_k,wherekis thetimeindex.Thestateinourcaseisavectorx_k∈R⁴,because therearefourhemodynamicstatevariables.TheobservedBOLD signalattimekisdenotedbyy_k∈R,becausetheBOLDsignalis onedimensional.Thestatetransitionnoisew_k isGaussianwith N(0,Qk),andtheobservationnoisevkisGaussianwithN(0,Rk).The setofparametersofthemodelisdenotedby.GivenNobservations y₁,y₂,···,y_Nory_1:N,ouraimistoﬁndtheparametersetandthe statesx1,x2,···,xNorx1:N.

InthefMRImodelinversionliterature,Rieraetal.[10]useda kindofextendedKalmanfilterwherethediscretizationwasnot basedontheEuler-Maruyamamethod.Theyusedthemethodpro- posedbyJimenezetal.[11].Huetal.,ontheotherhand,usedthe unscentedKalmanfilter(UKF)[12]intheirworkfortheestima- tionofthehemodynamicmodelstatesandparameters[13].Riera etal.andHuetal.usedthesetechniquesinaforwardpassmanner withoutsmoothing.TheearlyattemptsinthefMRImodelinversion literaturedisregardedthestatenoise.Forexample,in[7]Friston etal.modelstheinputandtheoutputrelationviaVolterraker- nels.Later,FristonusedaBayesianestimationtechnique[14]to findtheposterioroftheparametersagainunderzerostatenoise assumption.Later,Fristonet al.suggestedadvanced techniques basedonvariationalfilteringinwhichhealsoincludedthenonzero statenoisesinthemodel[15–17].Variationalfilteringtechniques assumefactorizationoftheparametersandstates.Dynamicexpec- tationmaximization(DEM)isthemostpopularofthem.AsHavlicek etal.pointedout,DEMalsoworksintheforwardpassmanner[3].

Johnstonetal.[2]usedparticleﬁlters,withoutusinganysmooth- ing.However,MurrayandStorkey[18]usedparticlesmoothers.

Mostrecently,Havlicek etal.[3]usedthesquare-rootcubature Kalmansmoother(SCKS)andobtainedquitearemarkablesuccess comparedwithDEMundercertainnoiseconditions.Mostofthe signalanalysistechniquesinthefMRIliteraturecanberegardedas akindofdeterministicmodelinversiontechnique.Finally,thereisa classofvariationalﬁlteringschemesknownasgeneralizedﬁltering thatoperateingeneralizedcoordinatesofmotion[17].Crucially, fromourperspective,theseeffectivelyperformonlinesmoothing.

Thisisbecausethefilterestimatesnotjustthecurrentstatebut highorder derivativesthatare basedon(local)pastand future states.Inotherwords,generalizedfilteringcanberegardedasa localBayesiansmoother.Generalizedfilteringis,however,arecent additiontotheportfolioofBayesiandeconvolutionschemesand willnotbeconsideredfurtherinthispaper.

Inthispaper,theiterativeextendedKalmansmoother(IEKS) methodisusedasamodelinversiontechniqueanditisshownthat theparameterandthestateestimationperformancesareimproved comparedwiththestandard extendedKalmantype algorithms, whichareused intheliterature.Quiteinterestingly,EKF inthe hemodynamicliteratureis usedonlyin theforwardpasswith- outsmoothing[2,10].Furthermore,Johnstonetal.[2]concluded thatEKFisnotrobustandpoorinperformancecomparedwiththe

particlefiltersandLLFfilters.TheassertionthatEKFispoorerinper- formancewasrepeatedinthefMRIliteratureinseverallandmark paperswithoutexaminingspecificallyforthehemodynamicmod- elingcase[3,16].Inthispaper,weshowedthatthecontraryistrue.

WefollowedasimilarapproachwithHavliceketal.’smethodand showedthatIEKSisapowerfulalternativeofformermodelinver- sionalgorithms.WecomparedIEKSwithEKSandEKFandnoticed theremarkablestateestimationperformanceunderdifferentpro- cess/measurementnoiseconditions.WecomparedIEKSwithSCKS forjointstateandparameterestimationperformance.Inallcondi- tions,IEKSisshowntoberobustandmoreaccuratethaniterative SCKS.Also,IEKSwasmorethantwofoldfasterthaniterativeSCKS.

1.1. Outline

Thepaperisorganizedasfollows: inSection 2,we givethe detailsofthealgorithmsused.ThestepsofEKF,EKSandIEKSare detailed,includingthejointparameterandthestateestimation parts.

HemodynamicmodelsystemidentiﬁcationispresentedinSec- tion 3;weworkonthedetailsofthefMRIhemodynamicmodel statespaceformulation.Weﬁrstrepresentthesysteminthedeter- ministicnonlineardynamicaldifferentialequationsform.Weuse theEuler-Maruyamamethodwithtimestept=0.1todiscretize thesystem.Weassumethatthereisanassociatingmeasurement foreveryonesecond.Weperturbthesystemwithprocessandmea- surementnoises.Weworkondifferentprocessandmeasurement noiseconditions,including[16,3].

In the Results and Discussion section, we first comparethe hemodynamic state estimation performance of EKF with the standardparticlefilters(PF).Johnstonetal.[2]assertedthatPF outperformsEKFandLLFfilters.Inthispaper,wereexaminedthis assertionandconcludedthatthemostbasicusageofparticlefilters donotoutperformotherfilteringtechniques.Wecheckedthestate estimationperformanceofPFwithEKFunderavarietyofprocess andmeasurementnoiseconditions.Subsequently,wecheckedthe numericalimprovementinstateestimation,achievedbyapplying smoothertotheEKFstateestimationforawiderangeofnoisecon- ditions.Subsequently,forthejointparameterandstateestimation, wecheckedtheperformancecomparisonofEKSandIEKS.Wecon- cludethatIEKSsubstantiallyimprovestheestimationperformance.

Subsequently,forthejointparameterandstateestimation,wegive thesimulationresultsofthealgorithmsIEKSanditerativeSCKSfor ﬁvedifferentscenarios.SCKSwastestedintheformerhemody- namicmodelinversionworkforlowprocessnoiseconditions[3].

WetestIEKSandSCKSformorechallengingprocessandmeasure- mentnoiselevels.TheperformanceimprovementofEKSoverSCKS forallprocessandmeasurementnoiseconditionsisshowninterms ofthestateandparameteraccuracy.Thebiasimprovementforthe parametersforIEKSisstressedcomparedwithiterativeSCKS.In Section 5,weapplyIEKSmethodtoa realdataset.Wediscuss furtheraspectsofthemethodindetail.

The main ﬁndings and contributions for the hemodynamic modelinversionforPF,EKF,EKSandIEKSareasfollows:

-TheassertionmadeintheliteraturethatPFoutperformsGauss- ian approximated model inversion methods seems not valid for the hemodynamic model tested for a wide range of process/measurementnoiseconditions.

-ThenotionthatEKFisnotanefﬁcientmodelinversionscheme canberefuted.Itappearstooperatewellforvariouslevelsof measurementnoiseanddifferentprocesses.

-EKF’shemodynamicstateestimationperformanceisevenbetter thanstandardparticleﬁltersunderawiderangeofprocessand measurementnoiseconditions.

(3)

-EKSisrobustandimprovesEKFespeciallyinthehigherhemoyd- namicprocess/measurementnoiseconditions.

-IEKSsubstantiallyimproveshemodynamicstateestimationper- formancebyutilizingtheEKSalgorithmiteratively.

-IEKS hasbetterhemodynamic stateestimationaccuracy com- paredwithiterativeSCKSfordifferentprocessandmeasurement noiseconditions.

-IEKShasbetterparameterestimationaccuracycomparedwith iterative SCKS for differenthemodynamic model process and measurement noise conditions. The biasof the parameters is lower.

-IEKSismorethantwicefasterthaniterativeSCKSforthejoint hemodynamicmodelinversion.

IntheformerevaluationofSCKSandDEM[3],SCKSisconcluded tobethebestmethodformodelinversioninthefMRIliterature.In thispaper,itisconcludedthatIEKSseemstobethebestmethod tobe usedfor the jointparameter and stateestimation ofthe hemodynamicdynamicalsystem.Thisalsostresseshowimportant smoothingisinthehemodynamicmodelinversionalgorithms.

2. Materialsandmethods

Inthissection,thedetailsoftheEKF,EKSandIEKSalgorithm isgiven.For theSCKSimplementation,thereaderisreferredto [3].TheextendedKalmanfilter/smootherandcubatureKalmanfil- ter/smootheralgorithmscanbeclassifiedunderthedeterministic inferencealgorithms.Theyapproximatethefilteredstateestimate p(x_k|y1:k)andthesmoothedstateestimatep(x_k|y1:N)byGaussian probabilitydensityfunctions(pdf).Since thepdfsareGaussian, bothofthealgorithmsiterativelyupdatethemeanandcovariance ofthestateestimates.

2.1. ExtendedKalmanﬁlter

FortheEKFalgorithm,thesystemislinearizedintheestimated statevalues.Afterthelinearizationstep,thestandardKalmanfilter algorithmisapplied,whichisknowntooperateoptimallyinlinear systems[19].SotheextendedKalmanfilterisakindofdeterminis- ticapproximationtechnique.Ithasthefurtherassumptionthatthe filteredandthepredictedstatescanbeapproximatedbyGaussian densities:

p(xk|y1:k)=N(ˆxk|k,P_k|k) (3)

p(x_k+1|y1:k)=N(ˆx_k+1|k,P_k+1|k) (4)

Here, ˆx_k|kand ˆx_k+1|karetheﬁlteredandthepredictedmeanvalues ofthedensities,whereasP_k|k andP_k+1|k aretheﬁlteredandthe predictedcovariancesofthestates.

Subsequently, the EKF algorithm recursively updates these statisticsasfollows:

Predictionupdate:

xˆ_k|k−1=f(ˆx_k−1|k−1) (5)

F_ij(x)=∂f(x_i)

∂x_j (7)

Thisconcludesthepredictionstep.Thenextstepisthemeasure- mentupdate.Theﬁrststepistocalculatetheso-calledinnovation vk.

Measurementupdate:

Ifthereisameasurementattimek,thentheupdatesare:

vk=y_k−h(ˆx_k|k−1) (8)

Sk=H(ˆx_k|k−1)P_k|k−1H(ˆx_k|k−1)^T+Rk (9)

K_k=P_k|k−1H(ˆx_k|k−1)^TS⁻¹_k (10)

xˆ_k|k= ˆxk|k−1+Kkvk (11)

P_k|k=P_k|k−1−KkSkK_k^T (12)

Ifthereisnomeasurementavailableattimek,

xˆ_k|k= ˆx_k|k−1 (13)

P_k|k=P_k|k−1 (14)

InordertocalculatetheKalmangain,K_k,themeasurementfunction hislinearizedsimilartothepredictionstep.TheJacobianmatrixof hatthestatevaluexisdenotedbyHwiththecomponentvalues evaluatedas:

Hij(x)=∂h(x_i)

∂xj

(15)

2.2. ExtendedKalmansmoother

ByusingtheEKFalgorithm,thestatex_kispredictedfromthe observedsequencey_1:k.Fromthefuturedata,oncewehavethem, a more accuratestate estimation can be accomplished.This is achievedthroughtheEKSalgorithm[1].TheEKSalgorithmisalso easytoimplement.Therecursionisdonebygoingbackwardin time.Therecursionstepbeginsfromk=N−1.

First,thestateispredictedas:

xˆ_k+1|k=f(ˆx_k|k) (16)

andthecovarianceispredictedbytheformula

P_k+1|k=F(ˆx_k|k)P_k|kF(ˆx_k|k)^T+Qk (17)

TheKalmangainK_kforthesmootheris:

K_k=P_k|kF(ˆx_k|k)^TP_k+1|k⁻¹ (18)

Smootherestimatesforthestatemean ˆx_k|N andstatecovariance P_k|NarefoundbyusingtheKalmangainK_kandtheﬁlterestimates ofthestatemean ˆx_k|kandthestatecovarianceP_k_|k:

xˆ_k|N= ˆxk|k+Kk(ˆx_k+1|N− ˆxk+1|k) (19)

P_k|N=P_k|k+Kk(P_k+1|N−P_k+1|k)K_k^T (20)

2.3. ParameterestimationanditerativeEKS

EKFand EKScanalsobeusedtoestimatetheparametersof thesystem.Itispossibletoaugmenttheparameterstothestate variables.Inthisdescription,theparametersaretreatedastime- varyingvariableswithsmallnoiseperturbationsasin[3].Sothe parameterupdatesbecome:

_k+1,1=_k,1+w_k,1 (21)

_k+1,2=k,2+wk,2 (22)

··· (23)

_k+1,p=_k,p+w_k,p (24)

(4)

where_k,istandsforthei-thparameterattimek,i=1,2,···,p.The newstatex_a,kbecomesanextendedversionoftheoriginalstatex_k.

xa,k= x_k

k,1

_k,2 .. .

_k,p

(25)

ForSCKS,Havliceketal.usedthesimilartechniqueforparame- terestimation[3].Sinceweaugmentedtheparameterstothestate, EKSalsoestimatesfortheparametersbesidesthestates.Attheend oftheEKSalgorithm,weobtainanestimateoftheparameterset.

Usingtheestimatesoftheparametersasthenewinitialestimate fortheparameters,weiteratetheEKSalgorithmuntilconvergence similarto[3].

3. Hemodynamicmodelsystemidentiﬁcation

Inthissection,weﬁrstreviewthedeterministichemodynamic model.FollowingFriston,wethenmakeanonlineartransformation toguaranteepositivevaluesforthehemodynamicvariables[16].

Wethenperturbthesystemwithwhitenoisetoobtainthestochas- tic representation of the hemodynamic model. As a ﬁnal step, wediscretizethesystembyusingtheEuler-Maruyamamethod toobtainthediscretetimenonlinearstatespacemodel.Wealso specifytheparametervaluesandthenoisestatisticsofthehemo- dynamicmodel.

Thehemodynamic model consists of four dimensionalstate transitionequationsandonedimensionalmeasurementequation.

Thehemodynamicmodelrepresentstherelationshipbetweenthe neuronalactivityuandthemetabolicdemands[18].Aftertheneu- ronalactivity,anincreaseinthevasodilatorysignalh₁isobserved.

In proportion to the increase in the vasodilatory signal h1, an increaseinthebloodﬂowh₂isobserved.Subsequently,achange intheblood volumeh₃ andthedeoxyhemoglobincontenth₄ is observed[16].Thesequantitativefactsareexpressedinthefollow- ingnonlineardynamicdifferentialequations:

˙h₁(t)=u(t)−h₁(t)−(h₂(t)−1) (26)

˙h2(t)=h1(t) (27)

˙h₃(t)=(h2(t)−F(h3(t))) (28)

˙h₄(t)=

h₂(t)E(h₂(t)))−F(h₃(t))h₄(t) h₃(t)

(29) where

F(h₃(t))=h₃(t)^1/˛ (30)

E(h₂(t))= 1

ϕ(1−(1−ϕ)^1/h²^(t)) (31)

Here,theparameterϕisdenotedastherestingoxygenextrac- tionfraction,andtheparameter˛iscalledGrubb’sexponent.The parameteriscalledthetransittime.Theparameterisamea- sureforthedecayofthevasodilatorysignalh1.Theparameteris therateconstantforthefeedbackofthebloodﬂowh₂.Finally,the parameterrepresentstheneuronalefﬁcacyfactor.

Themeasurementequationisgivenby y(t)=V0

k1(1−h4(t))+k2

1−h₄(t) h₃(t)

+k3(1−h3(t))

(32) HeretheparameterV₀iscalledtherestingbloodvolumefraction, andk₁,k₂andk₃arecalledtheintravascularcoefficient,theconcen- trationcoefficientandtheextravascularcoefficient,respectively.

Themeasurementequationrelatestheobserved BOLDsignalin termsofthebloodvolumeh3andthedeoxyhemoglobincontent h₄.Forthederivationandjustiﬁcationofthisequation,thereader isreferredto[20,9,21].

For arriving at the ﬁnal discrete time nonlinear state space model,we ﬁrstutilizethe techniqueof variabletransformation proposedbyFriston[22].Subsequently,weperturbthemodelby whitenoises.Afterward,weemploytheEulermethodtoderive thediscretetimenonlinearstatespacerepresentation.Inorderto guaranteepositivevaluesforthehemodynamicstates,wemakethe transformationx1(t)=h1(t)andx_i(t)=log(h_i(t))fori=2,3,4[16,3], andweaddthewhitenoises.Asaresult,weworkonthefollowing transformedsystemofequationsforthestates:

x˙1(t)=u(t)−x1(t)−(e^x²^(t)−1)+ˇ1(t) (33) x˙₂(t)=x₁(t)

e^x²^(t)+ˇ₂(t) (34)

x˙₃(t)=(e^x²^(t)−F(e^x³^(t)))

e^x⁴^(t) +ˇ₃(t) (35)

x˙₄(t)=(e^x²^(t)E(e^x²^(t)))−F(e^x³^(t))(e^x⁴^(t)/e^x³^(t)))

e^x⁴^(t) +ˇ₄(t) (36)

ByutilizingEulerdiscretizationfortheﬁrstderivative,weobtain thefollowingrepresentationofthestochasticnonlinearsystemfor hemodynamicmodelingoffMRIsignals:

x₁(t+t)≈x₁(t)+t(u(t)−x₁(t)−(e^x²^(t)−1))+

tˇ₁(t) (37)

x₂(t+t)≈x₂(t)+t

_x

1(t) e^x²^(t)

+

tˇ₂(t) (38)

x3(t+t)≈x3(t)+t

(e^x²^(t)−F(e^x³^(t))) e^x⁴^(t)

+

tˇ3(t) (39)

x4(t+t)≈x4(t)+t

(e^x²^(t)E(e^x²^(t))−F(e^x³^(t))(e^x⁴^(t)/e^x³^(t))) e^x⁴^(t)

+

tˇ4(t) (40)

Forthediscrete timeinstantst=t_kkt,k=1,2,...,putting x_k+1=x(t_k+t),x_k=x(t_k),u_k=u(t_k),y_k=y(t_k),andw_k=√

tˇ(t_k),¹ wearriveatthediscretetime statespacerepresentationofthe form:

x_k+1=f(,x_k)+w_k (41)

y_k=h(,x_k)+vk (42)

where thestatetransitionfunction f(,x_k)and theobservation functionh(,x_k)are

f(,xk)=

x_k,1+t(u_k−x_k,1−(e^x^k,2−1)) x_k,2+t

_x

k,1

e^x^k,2

x_k,3+t

(e^x^k,2−F(e^x^k,3)) e^x^k,4

xk,4+t

(e^x^k,2E(e^x^k,2)−F(e^x^k,3)(e^x^k,4/e^x^k,3)) e^x^k,4

(43)

1Here,thescalingfactor√

tonthenoisestandarddeviationcomesfromthe discretizationofthestochasticintegral

dˇt.

(5)

Table1

Parametervaluesofthehemodynamicstatetransitionequations.

Parameter Description Value

Rateofsignaldecay 0.65

Hemodynamictransittime 1.02

Rateofﬂowdependentelimination 0.41

˛ Grubb’sexponent 0.32

ϕ Restingoxygenextractionfraction 0.34

Neuralefﬁciency 0.5

Table2

Parametervaluesofthehemodynamicobservationequation.

Parameter Description Value

V0 Restingbloodvolumefraction 0.04

k1 Intravascularcoefﬁcient 7ϕ

k2 Concentrationcoefﬁcient 2

k3 Extravascularcoefﬁcient 2ϕ−2

and h(,x_k)=V₀

k₁(1−e^x^k,4)+k₂

1−e^x^k,4 e^x^k,3

+k₃(1−e^x^k,3)

(44)

respectively.Intheequationabove,x_k,ifori=1,2,3,4denotethe componentsofthestatevectorx_k.

Inthesimulations,wetooktheparametervaluesofthehemo- dynamicmodelandtheBOLDsignalasshowninTables1and2as in[16].Theinputistakenasin[3]asGaussianbumpswithdifferent amplitudescenteredinthetimepoints(10,15,39,48).Thelength oftheinputsignalis64s.tistakenas0.1sasin[2].Thesampling intervalforthemeasurementistakenas1s.

ThenoisestatisticscanbefoundinTable3.Forthemeasurement noisev_k,weﬁrstusedthenumericalvaluee⁻⁶asthestandarddevi- ationv.Notethatthiscorrespondstoanoisevarianceof_v²=e⁻¹², anditisanequivalentchoicetothenoisecharacteristicsgiveninthe histogramplotsin[2,10].Weincludedfourmorescenariostotest therobustnessandaccuracyoftheEKSandSCKSmethodsunder differentprocessandmeasurementnoiseconditionsascanbeseen fromTable3.Fortheparameterestimation,weestimatedthesame parameters,,asin[3].Asindicatedby[3],Grubb’sexponent

˛andtherestingoxygenextractionfractionϕwerealsofixedin theworkof[10].Thereasonforfocusingonthethreeparameters aboveisthatthereisconditionaldependencyamongtheparame- tersofthehemodynamicmodel.Thismeansthattheparameters consideredabovecanaccountfornormalvariationsinthehemo- dynamic responsefunction. In effect,this correspondsto fixing Grubb’sexponentandoxygenextractionfractiontotheirphysio- logicallyexpectedvalues[23,3].Furthermore,wealsoassumedthat theneuronalefficiencyisknownandtakenas0.5.IntheMonte

Table3

Noisestatisticsandinitialvalues.

Variable Value

InitialstatevalueX0 [0000]

Initialstatepdf N(0,0.01)foreachcomponent Initialparameterpdf N(true,1/12)foreachcomponent Parameternoisew N(0,10⁻⁵)foreachcomponent

Scenario1:wk,vk N(0,_w²=te⁻¹⁶)foreachstatecomponent, N(0,_v²=e⁻¹²)

Scenario2:w_k,vk N(0,_w²=te⁻¹²)foreachstatecomponent, N(0,_v²=e⁻¹²)

Scenario3:w_k,vk N(0,_w²=te⁻⁸)foreachstatecomponent, N(0,_v²=e⁻¹²)

Scenario4:wk,vk N(0,_w²=te⁻⁸)foreachstatecomponent, N(0,_v²=e⁻¹¹)

Scenario5:w_k,vk N(0,_w²=te⁻⁸)foreachstatecomponent, N(0,_v²=e⁻¹⁰)

Table4

RMSstateerrors:samplemeanandstandarddeviation.

Scenarionumber EKF PF EKS

Scenario1 0.0070±0.0031 0.0075±0.0030 0.0066±0.0029 Scenario2 0.0095±0.0026 0.0098±0.0026 0.0092±0.0023 Scenario3 0.0408±0.0034 0.0411±0.0035 0.0344±0.0028 Scenario4 0.0433±0.0041 0.0434±0.0042 0.0381±0.0036 Scenario5 0.0454±0.0051 0.0455±0.0050 0.0423±0.0048

Carloanalysis,wedisturbedtheinitialestimationoftheparameters

,,bythevarianceN(0,1/12)asdonein[3].

Basedontheaboveproblemformulation,inthenextsection,we performadetailedcomparisonofPF,EKS,IEKS,anditerativeSCKS underﬁvedifferentsimulationscenarios.

4. Results

4.1. Simulationresults

Inthissection,weperformMonteCarlosimulationsunderfive differentscenarioswithdifferentprocessandmeasurementnoise conditions.Theground-truthvaluesforalltheparametersofthe hemodynamicmodeldefinedinEqs.(41)–(44)aretakenasshown inTables1and2[22].Foreachscenario,100MonteCarlorunsare performed.Inthefirstsetofscenarios1,2,and3,wefixthemea- surementnoisestandarddeviationto_v=e⁻⁶asin[3]andchange processnoiselevels.Inscenario1,theprocessnoisestandarddevi- ationisasin[3](i.e.,w=√

te⁻⁸).Inscenarios2and3,wetry morechallengingprocessnoiseconditionswithw=√

te⁻⁶and w=√

te⁻⁴.Inthesecondsetofscenarios,3,4,and 5,weﬁx theprocess noise standard deviationtow=√

te⁻⁴,which is thevalueinscenario3,andchangethemeasurementnoiselev- els.Weincreasethemeasurementnoisevarianceto²_v =e⁻¹¹and _v²=e⁻¹⁰forscenarios4and5,respectively.

A prioriinformation aboutthestates are takenasshown in Table3similarto[3].Thestatecomponentsareallinitializedwith 0.

Inthesimulation,forsomerarecases,weobservedthatthesec- ondstatex₂cangoto−∞.Thiscorrespondsto0fortheoriginal untransformedvariableh₂.Sincethetransformationisoftheexpo- nentialtype,wesetalowerlimitx2=−4(h2=0.0183)foreachtime stepattheﬁlteringstepofEKS.Thesamelimitisalsoappliedtothe othervariables.ThesamethingisalsodoneforPF.PFisevenmore pronetodivergetheselimits.Thereasonisthefollowing:EKFtracks themeanofthestate,whichisexpectedtobeinthestableregion.

ButthePFtriestorepresenttheconditionalpdfp(x_k|y1:k),including theextremeconditionsalso.Forthatreason,theyaremoreprone todivergeforthoseparticles.ThesamelimitsarealsoputforSCKS.

4.1.1. Comparisonoftheﬁlteringalgorithms

Inthissection,wewanttoshowthatthecommonlyheldview thatEKFisnotanappropriatefilteringalgorithmisnotright.John- stonetal.[2]comparedEKFstateestimationwithparticlefilter andconcludedthatPFoutperformsEKF,anditisstatedthatEKF divergesmostofthetime,andasaresult,itisnotrobustandpoor inperformance.Weseeinthissectionthat,onthecontrary,not onlystateestimationbutalsojointestimationofparameterand stateisrobustandperformsbetterthanPF.Inthissubsection,we firstcomparetheEKFandPFstateestimationandshowthatEKF performsbetter.InTable4,theMonteCarloresultsofEKFwithPF arecompared.Inallscenarios,EKFperformsbetterthanPF,andit isrobust.Theparticlenumberischosenas500,whichisfivefold morethanthecasewheretheoriginalcomparisonwasmade[2].

(6)

0 2 4 6 8 10 12 14 16 18 20 0

0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02

RMS State Error

Iteration Number

RMS State Error vs Iteration Number

IEKS

Fig.1.IEKSandimportanceofiteration.

4.1.2. PerformanceimprovementbyEKSoverEKF

In the fMRI state estimation literature, extended Kalman- type estimation algorithmsare only used inthe filteringsense [2,10]. Eventhe parameterestimation algorithms which rested onextendedKalmanfilteralgorithms,relyonthefilteringalgo- rithm[10].Inthissubsection,assumingthatwehaveknownand fixedparameters,ouraimistoestimatethehiddenhemodynamic statesbyusingextendedKalmanfilterandsmootheralgorithms.

WeshowtheperformanceimprovementinTable4.Wesummarize theEKFandEKSstateestimationerrorsinRMS.Forhighprocessand measurementnoiseconditions,theimprovementismoreapparent.

4.1.3. JointparameterandstateestimationwithiterativeEKS Aprioriinformationabouttheparametersaretakenasshownin Table3similarto[3].Thestatecomponentsareallinitializedwith 0.Notethattheinitialvaluesofalltheparametersareassumedto beGaussiandistributedaroundtheirtruevalueswithavarianceof 1/12asin[3].Intheliterature,Riera[10]usedextendedKalmanﬁl- terwiththediscretizationproposedbyJimenezetal.[11]without anyiteration.ThisisthestandardusageofEKFandEKS.Forthat reason,inthissectionwealsocomparetheEKSstateRMSerror withIEKSutilizedinthispaper.Fig.1alsoshowstheimportanceof theiteration.Thisﬁgureisatypicalrepresentativeforshowingthe decreaseofRMSstateerrorwithrespecttotheiterationnumber.

TheexampleplotisaparticularMCmeanestimateforscenario2.By increasingtheiterationnumber,theRMSstateerrorsaredecreased substantially.Thereasonthatthestateestimationimproveswith subsequentiterationsisthattheparameterestimatesbecomepro- gressivelymoreaccurate;therebyenablingabetterestimationof thestates.Obviously,thisperformanceshowsthatiterativeEKS outperformsthedirectusageofEKSwithoutiteration.

We alsosummarized theMonte Carloparameter estimation resultsofthealgorithmsIEKSanditerativeSCKSinTables5–7.The truevaluefortheparameterswere=0.65,=1.0204and=0.41.

In all ﬁve different process/measurement noise levels, for all parametervalues,theaccuracyoftheestimationofEKSwasbetter thanthatofSCKS.Onlyinonecaseweretheyequal.Wenotethat thebiasoftheEKSMonteCarloestimateislessthanthebiasofthe SCKScase.Similarly,forthestateestimation,theerrorinRMSis summarizedinTable8.Althoughthedifferencewasnotbig,forall cases,IEKSwasbetterthantheSCKSmethod.Furthermore,itwas muchfasterthanSCKS.Bothalgorithmsarerobustunderdifferent

Table5

Simulationresultsfortheparameterestimates.

Scenario EKSEstimate SCKSEstimate EKSBias SCKSBias Scenario1 0.6489±0.0282 0.6511±0.0280 0.0011 0.0011 Scenario2 0.6494±0.0289 0.6517±0.0288 0.0006 0.0017 Scenario3 0.6545±0.0556 0.6580±0.0556 0.0045 0.0080 Scenario4 0.6561±0.0627 0.6588±0.0630 0.0061 0.0088 Scenario5 0.6560±0.0748 0.6571±0.0752 0.0060 0.0071

Table6

Table7

measurement noise conditions. As expected, by increasing the measurementnoise,thegradualdecreaseoftheperformanceof theestimatesisobserved.Thebiasandsamplevarianceestimates areincreasedbyincreasingtheprocessand measurementnoise variance.

Table8

RMSstateerrors:samplemeanandstandarddeviation.

Scenarionumber EKS SCKS

Scenario1 0.0128±0.0038 0.0130±0.0039

Scenario2 0.0140±0.0035 0.0143±0.0036

Scenario3 0.0374±0.0046 0.0376±0.0047

Scenario4 0.0418±0.0053 0.0420±0.0054

Scenario5 0.0483±0.0071 0.0486±0.0074

(7)

0 50 100 150 200 250 300 350

−10

−8

−6

−4

−2 0 2 4 6 8 10

time bin with Δ t =3.22 sec.

BOLD Signal

BOLD Signal vs. time

BOLD Signal

Fig.2.RealBOLDdataiscollectedfromthemotionsensitiveareaV5.TimeintervaltisinthedurationofthesamplingtimeTR=3.22s.Theﬁrst256sampleswereusedfor modelinversion.

4.2. Discussion

InthefMRIinversionliteratureakeypaper[2]hasassertedthat EKFis notarobustmethodforstateestimation,and itsperfor- mancewaspoorcomparedwithPF.Inthispaper,weexamined EKFfor a varietyofnoise conditionsand concludedthatthis is notthecase.Itevenperformedbetterthanparticlefilter.Wetake fivefoldmoreparticlesthanthepaper,whichmakesthatasser- tion.Standardparticlefilterswereallegedtoperformbetterthan theotherhemodynamicstateestimationalgorithms.Inthispaper, weconcludedthat,onthecontrary,EKFperformedbetterthanthe otherfilteringalgorithms.Thereasonisthatstandardparticlefilters useanonoptimalproposalfunction,whichdegradestheperfor- mance.Thereisstillroomforparticlefilters,whichmayusemore sophisticatedproposalfunctions.

In thispaper,the importanceofthesmoothing for boththe hemodynamic state and the parameter estimation is empha- sized.TheextendedKalmanfilterismodifiedbyincorporatingthe smoother.For fixed hemodynamicparameters, we checkedthe stateestimationimprovementforvariousnoiseconditions.Espe- ciallyforhigherconditions,theimprovementwasmoreapparent.

WenotethatextendedKalman-typealgorithmsinthefMRIlitera- turewerenotusedwithsmoothinguptonow[10,2].

Thestandardapplicationoftheﬁlteringandsmoothingisper- formedonlyonce.ByiterativelycallingEKFandEKSaftereachnew parameterestimate,wenoticedthehugeimprovementovernon- iterativeusage.ThisimprovementisshowninFig.1.EKSisjust thespecialcaseofIEKSwithjustoneiteration.LikeHavliceketal.

[3],treatingtheparametersastime-varyingvariablesbyadding artiﬁcialdynamics,thejointstateandparameterestimationofthe hemodynamicmodelisperformedbyEKSinthispaper.Havlicek etal.compared theiterativeusageofSCKSwithDEMand con- cludedthatSCKSismoreaccurateintheirpaper[3].Inthisreport, IEKSiscomparedwithiterativeSCKSunderdifferentnoisecondi- tions.Theﬁrstscenarioweusedhasthesamenoiseconditionsas in[3].Thismodelhasaverylowprocessnoisecovariance.Evenin thiscondition,ourmethodwasmoreaccurateforboththestate andtheparameterestimation.Overall,bothmethodsweregood.

Forlowprocessandmeasurementnoise,itisquitereasonableto approximatethedensitiesbyGaussianpdf,becausethesystemis almostdeterministic.Asa result,weexpectthat EKSand SCKS performwell.However,innonlinearsystems, eventheadditive

noiseisGaussian;theconditionalpdfsp(x_k|y1:k)andp(x_k+1|y1:k) actually deviate fromthe Gaussian assumption.The higher the covariancesofthenoises,themoreseveretheseconditionalpdfs deviatefromtheGaussianassumption.Wetestedthealgorithms underﬁvedifferentprocess/measurementnoiseconditions.IEKS wasmoreaccuratethaniterativeSCKSforbothparameterandstate estimationforallcases.

BesidestheeaseofimplementationofEKS,itwasalsomuch fasterthanSCKS.EKSisaround2.3timesfasterthanSCKS.

Sooverall,asasserted,EKFisnotanunreliablemodelinversion method.ItisbetterthantheallegedmostsuccessfulalgorithmPF.

TheiterativeEKSissubstantiallymoreaccuratethantheformerly usedEKF-typealgorithmsinthefMRImodelinversionliterature.

EKSisrobustunderawiderangeofnoiseconditions.EKSisfaster andhaslowerparameterbiasandmoreaccuratestateestimation thantheSCKSmethod,whichseemstoworkbestamongthecur- rentfMRImodelinversionmethods[3]forallthenoiselevelswe worked.

5. Applicationtorealdata

5.1. Backgroundinformation:testsetupandBOLDdata

Inthissection,wewanttoshowthevalidityoftheIEKSmethod onarealBOLDdataasshowninFig.2.Thisdatawasalsousedby Fristonetal.toshowthepracticalapplicabilityoftheirmodelinver- sionmethodsDEM,variationalﬁlter(VF)andgeneralizedﬁlter(GF) [15–17].²

TheBOLDdatawascollectedfromthemotionsensitiveareaV5 [15–17].Duringthetest,thesubjectwasexposedto5differentcon- ditions.Theﬁrstonewastoallowformagneticsaturationeffects [15–17].Inthesecondonecalled“Fixation”,subjectwasviewing aﬁxationpointinthescreen.Inthethirdonecalled“Attention”, subjectwasviewing250dots.Thosedotsweremovingradiallyby 4.7^◦fromthecenter.Atthesametime,thesubjectwassupposed todetectthechangesinradialvelocity.Inthefourthonecalled

“NoAttention”,thesubjectjustviewedthemovingdots.Intheﬁfth one,thepersonwassupposedtoviewthestationarypoints.During

2WethankProfessorKarlFristonforhiskindnesstogiveuspermissiontousethe BOLDdata.

(8)

0 1000 2000 3000 4000 5000 6000

−1

−0.5 0 0.5 1 1.5

Neuronal Input

time bin with Δt=0.2

Neuronal Input vs. time

Neuronal Input − Vision Neuronal Input − Motion Neuronal Input − Attention

Fig.3.Inthissetup,thepersonwasexposedto3differentneuronalinputs.TheyareVision,MotionandAttention.Thoseneuronalinputsaremodeledasboxcarfunctions.

Duringthetest,thepersonwassubjectto5differentconditions.ThoseareSaturation,Fixation,Attention,NoAttentionandMotionView.Dependingonthecondition,for eachneuronalinput,wehavethevalueofeither1or0.However,followingFristonetal.[15–17],DCvaluesoftheinputsareremoved.

Fig.4. Inthisﬁgure,weshowthesimulationswhichconvergetoglobaloptimag.Speciﬁcally,weestimatetheparameterset=[₁,₂,₃,,,].Theconvergedglobal optimaisg=[0.1024,0.2102,0.0175,0.7285,0.4981,0.6460].

theentiretest,thesubjectwasexposedtotheﬁxationandvisual stimulusconditionsinanalternatively.

Fristonetal.[15–17]modeledthesefivedifferentconditionsas acombinationof3differentneuronalinputs.Thefirstoneisvisual stimulus,thesecondoneismotionstimulusandthethirdoneisthe attentionstimulusrepresentedbyu1,u2andu3,respectively.Each neuronalinputismodeledbyboxcarfunctions.Forexamplefor thethirdcondition“Attention”,inthosetimeinstantsu₁,u₂andu₃ haveallthevalue1.Thereasonisthatinthe“Attention”condition thepersonissupposedtoview,trackthemotionandshallattend thechangesinradialvelocities.Inthe“Noattention”conditionu₁, u₂ andu₃ are1,1and0,respectively.Thereasonisthatinthis conditionthepersonviews,tracksthemotionbutthereisnothing forattention.Inthefifthconditionu₁,u₂andu₃are1,0and0.The reasonisthatthereisvisualstimulus,butthedotsarenotmoving andthereisnothingfortheattention.Similarly,inthe“Fixation”

conditionu₁,u₂andu₃ areall0.OnefurthernoteisthatFriston etal.[15–17]havetherealBOLDsignalwithDCvalue0,forthat

reasonheremovestheDCvaluesfromu₁,u₂andu₃.Asaresult, theyusedtheneuronalinputsequencesdepictedinFig.3[15–17].

5.2. Extendedhemodynamicmodelformulti-inputsystem

Sincewehaveinthissetuparathercomplicatedmulti-input system,weextendthenonlinearhemodynamicsystemtothefol- lowingnonlineardifferentialequationsystem.

˙h₁(t)=1u₁(t)+2u₂(t)+3u₃(t)−h₁(t)−(h₂(t)−1) (45)

˙h2(t)=h1(t) (46)

˙h₃(t)=(h₂(t)−F(h₃(t))) (47)

˙h₄(t)=

h₂(t)E(h₂(t))−F(h₃(t))h4(t) h3(t)

(48)

Here1,2and3 aretheneuronalefﬁciencyfactorsforthe neuronalinputsu1,u2andu3,respectively.Inthissetup,1,2and

(9)

0 20 40 60 80 100 120 140 160 180 200 0.18

0.185 0.19 0.195 0.2 0.205 0.21 0.215

Iteration number

RMS Error

RMS Error vs Iteration Number, Global Optimum Case

RMS Error

Fig.5. Inthisﬁgure,weshowthesimulationswhichconvergetoglobaloptima.Foreveryiteration,theRMSerrorofthepredictionisshown.Ateverystep,theRMSerrors aredecreasing.Afteraround25iterations,thealgorithmconverges.

Fig.6. Inthisﬁgure,weshowthesimulationswhichconvergetolocaloptima.Speciﬁcally,weestimatetheparameterset=[1,2,3,,,].Theconvergedlocaloptima isl=[0.0329,0.0776,0.0091,0.6349,2.2990,0.2381].

3 giveinformation howmuchtheselectedbrainregionissen- sitivetovisioninputu₁,motioninputu₂andattentioninputu₃, respectively.

Inordertohaveabestestimatefortheparameterset=[1,2,

3,,,].WeperformedaMonteCarlosimulationandbeganwith differentinitialparametersets.Theinitialparametersetandnoise conditionsarechosenaccordingtoTable9.

Table9

Noisestatisticsandinitialvalues.

Variable Value

InitialstatevalueX0 [0000]

Initialstatepdf N(0,0.01)foreachcomponent Initialparameterpdf N(0,1/12)for₁,₂,₃

Initialparameterpdf N(0.65,1/12),N(1.02,1/12),N(0.41,1/12) for,,respectively

Parameternoisew N(0,t10⁻⁸)foreachcomponent Processnoise:w_k N(0,_w²=te⁻⁸)foreachstatecomponent Measurementnoise:vk N(0,_v²=e⁻¹²)

5.3. Results:IEKSmodelinversion

Inthis section,weapplytheIEKS methodtotheBOLDdata.

WeperformMonteCarlosimulationwith100runstoestimatethe parameters.Theinitialparameterset=[1,2,3,,,]ischosen accordingtothepriorstatisticsgiveninTable9.Weiterateuntil thealgorithmconverges.TheIEKSalgorithmconvergedto2values whichareglobalandlocaloptima.Wecallthetwoparameterset asgand_l.

g=[0.1024,0.2102,0.0175,0.7285,0.4981,0.6460] (49)

l=[0.0329,0.0776,0.0091,0.6349,2.2990,0.2381] (50) Whenwenotetheneuronalefﬁciencyfactors1,2 and3,we observecloseresemblanceoftheconvergedsequencestotheones reportedbyFristonetal.[15–17].Bothgand_lshowstheneuronal activationmostfor2whichcorrespondsforthe“Motion”condition.Somewhatfor1whichcorrespondstothe“Vision”condition andalmostnothingforthe“Attention”condition.Thedifferenceis inthestrength.

(10)

0 20 40 60 80 100 120 140 160 180 200 0.195

0.2 0.205 0.21 0.215 0.22 0.225 0.23

Iteration number

RMS Error

RMS Error vs Iteration Number, Local Optimum Case

RMS Error

Fig.7. Inthisﬁgure,weshowthesimulationswhichconvergetolocaloptima.Foreveryiteration,theRMSerrorofthepredictionisshown.Ateverystep,theRMSerrors areincreasing.Afteraround200iterations,thealgorithmconverges.

0 100 200 300 400 500 600

−5 0 5 10

time index

BOLD Signal

BOLD Signal vs time

input BOLD Signal

Fig.8.BOLDsignalchangewithrespecttochange.ischangedintheinterval[0.35;1.05].

0 100 200 300 400 500 600

−2 0 2 4 6 8 10

time index

BOLD Signal

BOLD Signal vs time

input BOLD Signal

Fig.9.BOLDsignalchangewithrespecttochange.ischangedintheinterval[0.32;2.32].

(11)

0 100 200 300 400 500 600

−2 0 2 4 6 8 10 12 14

time index

BOLD Signal

BOLD Signal vs time

input BOLD Signal

Fig.10.BOLDsignalchangewithrespecttochange.ischangedintheinterval[0.21;1.01].

0 100 200 300 400 500 600

−5 0 5 10 15

time index

BOLD Signal

BOLD Signal vs time

input BOLD Signal

Fig.11.BOLDsignalchangewithrespecttochange.ischangedintheinterval[0.01;0.51].

Foreachsimulationandateveryiteration,RMSerrorintheBOLD signalpredictionisevaluatedaccordingtotheformula.

erms=

¹

N

k=1

y_k−y^p_k

² (51)

wherey_kandy^p_kistherealandpredictedBOLDsignal,respectively.

RMSerrorsforglobaloptimumvaluearelowerthantheones forlocaloptimum.Furthermore,ateveryiteration,theglobalopti- mumRMSvaluesareimproved,whereasthelocaloptimumRMS errorsaregettingworse.Theparameterestimateswhichconverge toglobaloptimaandlocaloptimaareshowninFigs.4and6.Simi- larlytheRMSerrorsinpredictioncanbeseeninFigs.5and7.

5.3.1. ParametersweepandtheBOLDsignal

In this section, we want to observe the changes in the hemodynamicstatevariables and theBOLDsignalwithrespect

to the parameter changes. We let the parameter values as in Tables1and2.Wesweeptheparametersforintheinterval[0.35 1.05],intheinterval[0.322.32],intheinterval[0.211.01], intheinterval[0.010.51].Thoseintervalscontainthewidelyused parameterregimesusedintheliterature[15–18,10,13].Figs.8–11 theBOLDsignalchangewithrespecttotheparameters,,,, respectively.

5.4. Discussion

Inthissection,wediscusstheresultsindetail.Fortheglobal optima,theIEKSmethodconvergedtotheparameterset_g.When wechecktowhichinputs theselectedregionismostsensitive, weobservethat2is,asexpected,signiﬁcantlybiggerthan1and

3.ThisresultisincloseagreementwithFristonetal.[15].They usedDEMformodelinversion.Ashere,theyhaverelativelyhigh

2 estimatedregion.WealsovisualizetheresultsinFig.12.We

extended smoother Kalman Joint and state parameter estimation of the hemodynamic model byiterative Biomedical Signal Processing and Control

Biomedical Signal Processing and Control

Joint parameter and state estimation of the hemodynamic model by iterative extended Kalman smoother

Serdar Aslan

, Ali Taylan Cemgil

, Murat S¸ amil Aslan

, Behc¸ et U˘gur Töreyin

, Ata Akın



























