´ ´I OPTIMALIZACN ˇ ´ICH ALGORITMU ˚ PRO SROVNAN ´ ´ ´ ULTRAZVUKOVOU PERFUZNI ANALYZU METODOU BOLUS & BURST M. M´ezl 1 , R. Jiˇr´ık 2 1 Ustav ´
´ ˇ v.v.i. biomedic´ınsk´eho inˇzen´ yrstv´ı, VUT v Brnˇe; 2 Ustav pˇr´ıstrojov´e techniky AV CR, Abstrakt
Pˇ redloˇ zen´ a studie pojedn´ av´ a o srovn´ an´ı optimalizaˇ cn´ıch algoritm˚ u v Optimization toolboxu v Matlabu. Na simulovan´ ych datech pro ultrazvukovou metodu perf´ uzn´ı anal´ yzy bolus & burst s r˚ uzn´ ymi u ´ rovnˇ emi ˇ sumu je testov´ ana pˇ resnost odhadu arteri´ aln´ı vstupn´ı funkce. Optimalizaˇ cn´ı probl´ em je definov´ an jako parametrick´ a slep´ a dekonvoluce, na kterou je pohl´ıˇ zeno jako na minimalizaci stˇ redn´ı kvadratick´ e chyby. Na z´ akladˇ e v´ ysledk˚ u t´ eto studie jsou stanovena doporuˇ cen´ı pro anal´ yzu dat z preklinick´ ych a klinick´ ych studi´ı.
1
´ Uvod
Metoda bolus & burst (B&B) je pomˇernˇe nov´a metoda ultrazvukov´e perf´ uzn´ı anal´ yzy (UPA), kter´a kombinuje dva z´ akladn´ı a pouˇz´ıvan´e pˇr´ıstupy pro perf´ uzn´ı anal´ yzu – tzv. bolus tracking a burst replenishment [1, 2]. Pomoc´ı t´eto metody je moˇzn´e prov´est plnˇe kvantifikaˇcn´ı anal´ yzu perf´ uzn´ıch parametr˚ u - toku krve, objemu krve na mikrokapil´ arn´ı u ´rovni nebo stˇredn´ı doby pr˚ uchodu. Metody oznaˇcovan´e jako bolus tracking vyuˇz´ıvaj´ı aplikaci ultrazvukov´e kontrastn´ı l´atky (KL) injekc´ı a n´ızkoenergetick´e zobrazov´an´ı oblasti z´ajmu. V r´ amci perf´ uzn´ı anal´ yzy jsou vyhodnocov´any kˇrivky ultrazvukov´e intenzity (z´ıskan´e jako stˇredn´ı hodnota v oblasti z´ajmu), kter´e jsou v line´ arn´ım vztahu ke koncentraci KL. Toto vyhodnocen´ı je nejˇcastˇeji zaloˇzeno na proloˇzen´ı namˇeˇren´ ych dat statistick´ ym rozdˇelen´ım (lognorm´aln´ı, gamma a dalˇs´ı) [8]. Tento pˇr´ıstup nejˇcastˇeji vede ke odhadu parametr˚ u, kter´e maj´ı relaci ke skuteˇcn´ ym perf´ uzn´ım parametr˚ um, ale nen´ı moˇzn´e je fyzik´alnˇe kvantifikovat. Plnˇe kvantifikaˇcn´ı anal´ yza je moˇzn´ a pomoc´ı aplikace perf´ uzn´ıho modelu zaloˇzen´eho na konvoluci arteri´ aln´ı vstupn´ı funkci (AIF) a rezidu´aln´ı funkce tk´anˇe (TRF) [5]. Druh´ a skupina metod oznaˇcovan´ a jako burst replenishment vyuˇz´ıv´a aplikace KL pomoc´ı intraven´ ozn´ı inf´ uze a po dosaˇzen´ı ust´ alen´e koncentrace kontrastn´ı l´atky v krevn´ım obˇehu dojde k aplikaci sekvence vysokoenergetick´ ych pulz˚ u (tzv. burst), kter´e v oblasti z´ajmu zniˇc´ı mikrobubliny KL. V r´ amci perf´ uzn´ı anal´ yzy je pot´e vyhodnocv´ano znovuzaplnˇen´ı oblasti z´ajmu kontrastn´ı l´atkou (tzv. replenishment) [10]. Plnˇe kvantifikatˇcn´ı anal´ yza je u t´eto skupiny metod moˇzn´ a s vyuˇzit´ım informace o intenzitˇe odraˇzen´eho ultrazvuku v krevn´ım ˇreˇciˇsti. Metoda bolus & burst kombinuje oba v´ yˇse uveden´e principy do jednoho [1, 2] Akvizice dat prob´ıh´ a podle n´ asleduj´ıc´ıho sch´ematu. V prvn´ı f´azi (bolus f´aze)je aplikov´ana KL injekc´ı (bolus) s n´ asledn´ ym sn´ım´ an´ım s n´ızk´ ym mechanick´ ym indexem (MI) oblasti z´ajmu. Po dosaˇzen´ı u ´st´ alen´e hladiny kontrastn´ı l´ atky v obˇehu je aplikov´an burst, kter´ y v oblasti z´ajmu zniˇc´ı mikrobubliny KL a d´ ale je provedeno sn´ım´ an´ı opˇetovn´eho zaplnˇen´ı oblasti z´ajmu KL pomoc´ı zobrazov´an´ı s n´ızk´ ym ˇ MI (replenishment f´aze). Casov´ y pr˚ ubˇeh koncentrace je modelov´an jako konvoluce AIF a TRF, kter´e jsou formulov´any parametricky. Odhad parametr˚ u je proveden slepou dekonvoluc´ı. Metoda byla testov´ana na klinick´ ych datech pacient˚ u s Crohnovou chorobou [2, 7] a preklinick´ ych datech n´ adorov´ ych myˇs´ı [3]. Odhad perf´ uzn´ıch parametr˚ u u metody B&B je zaloˇzen na slep´e dekonvoluci, kter´a je formulov´ana jako minimalizace stˇredn´ı kvadratick´e chyby mezi mˇeˇrenou a modelovanou kˇrivkou ´ koncentrace. Ukolem t´eto pr´ ace je srovn´an´ı bˇeˇznˇe dostupn´ ych optimalizaˇcn´ıch metod v Optimization Toolbox v Matlabu na simulovan´ ych datech.
2 2.1
Materi´ aly a metody Optimalizaˇ cn´ı probl´ em
ˇ Casov´ y pr˚ ubˇeh koncentrace kontrastn´ı l´atky v oblasti z´ajmu (ROI), CR (t), je pops´ an konvoluc´ı arteri´ aln´ı vstupn´ı funkce (AIF) a rezidu´aln´ı funkce tk´anˇe (TRF) jako [6] Z ∞ AIF (τ )R(t − τ )dτ, (1) CR (t) = Fb 0
kde parametr Fb je tok krve v ml/min/ml tk´anˇe. Parametrick´e vyj´adˇren´ı TRF, R(t), je v nejjednoduˇsˇs´ı podobˇe definov´ano jako klesaj´ıc´ı exponenci´ aln´ı funkce 0 t≤0 R(t) = −βt (2) e t > 0, kde β je jedin´ y odhadovan´ y parametr. Perf´ uzn´ı parametr stˇredn´ı doba pr˚ uchodu (MTT) je potom definov´an jako pˇrevr´ acen´ a hodnota tohoto parametru (M T T = 1/β). Parametrick´e vyj´adˇren´ı AIF je definov´ano jako souˇcet tˇr´ı gamma funkc´ı, AIF3gam , [3] AIF3gam (t) = t
α
3 X
Ai e−τi t ,
(3)
i=1
kde α, Ai a τi jsou parametry modelu (celkem 7 parametr˚ u). Odhad parametr˚ u modelu prob´ıh´ a pomoc´ı slep´e dekonvoluce, kter´a je definov´ana jako minimalizace stˇredn´ı kvadratick´e chyby mezi mˇeˇrenou kˇrivkou koncentrace v ROI a jej´ım modelem podle rovnice (1) s parametrick´ ymi vyj´adˇren´ımi TRF, resp. AIF, podle (2), resp. (3). Miminalizace je provedena s diskretizovanou ˇcasovou osou, tzn. nez´ avisl´a promˇenn´ a t je nahrazena indexy n a m. Minimalizace je realizov´ana pro obˇe ˇc´asti sign´ alu - f´azi bolusu i f´azi replenishmentu jako [3]: w1 +w2
Xb1
XN −1
n=0
n=b2
Cm (n) − Fb Cm (n) − Fb
Xb1
m=0
XN −1
2 AIF (m)R(n − m) +
m=b2
2 AIF (m)R(n − m) .
(4)
ˇ Casov´ e indexy b1 a b2 vymezuj´ı aplikaci burstu do oblasti z´ajmu, Cm (n) odpov´ıd´ a mˇeˇren´e kˇrivce koncentrace v ROI. V´ahy w1 a w2 upravuj´ı vliv jednotliv´ ych ˇc´asti kˇrivky.
2.2
Simulovan´ a data
Pro u ´ˇcely t´eto studie byl vygenerov´an soubor simulovan´ ych dat zaloˇzen´ y na re´ aln´ ych mˇeˇren´ıch [6]. Uk´ azkov´ y namˇeˇren´ y pr˚ ubˇeh AIF byl parametrizov´an pomoc´ı rovnice (3) a takto z´ıskan´ a referenˇcn´ı AIF byla pouˇzita pro generov´an´ı simulovan´eho pr˚ ubˇehu podle rovnice (1) s definovan´ ym M T T = 5 s. Takto z´ıskan´ y pr˚ ubˇeh byl zaˇsumˇen aditivn´ım ˇsumem s Gaussovsk´ ym rozloˇzen´ım hodnot tak a nulovou stˇredn´ı hodnotou. Byl stanoven pomˇer sign´ al–ˇsum jako pomˇer mezi v´ ykonem nezaˇsumˇen´eho sign´ alu a v´ ykonem ˇsumu. Soubor simulovan´ ych dat obsahoval pˇet u ´rovn´ı ˇsumu odpov´ıdaj´ıc´ı SNR 10, 20, 30, 40 a 50 dB. Pro kaˇzdou u ´roveˇ n bylo generov´ano 100 realizac´ı n´ ahodn´eho ˇsumu. Celkem tak datov´ y soubor obsahuje 500 unik´ atn´ıch pr˚ ubˇeh˚ u. Uk´ azka simulovan´ ych dat je na Obr. 1.
Obr´azek 1: Pˇr´ıklad simulovan´ ych dat - referenˇcn´ı AIF (vlevo nahoˇre), TRF s nastaven´ ym para−1 metrem β = 0, 2 s (vpravo nahoˇre), kˇrivka CR s u ´rovn´ı ˇsumu 40 dB (vlevo dole) a kˇrivka CR su ´rovn´ı ˇsumu 10 dB (vpravo dole).
2.3
Optimalizaˇ cn´ı metody
Pro minimalizaci krit´eria (4) byl pouˇzit Optimization Toolbox v Matlabu ve verzi R2015a. Na samotnou optimalizaci bylo pohl´ıˇzeno jako na nepodm´ınˇenou optimalizaˇcn´ı u ´lohu (unconstrained optimization) i na u ´lohu podm´ınˇenou (constrained optimization). Omezen´ı (podm´ınky) byla definov´ana jako minim´aln´ı a maxim´ aln´ı hodnoty, kter´ ych mohou jednotliv´e parametry nab´ yvat. Tato omezen´ı byla volena heuristicky na z´akladˇe pˇredchoz´ıch simulac´ı. Pro srovn´an´ı bylo testov´ano celkem pˇet algoritm˚ u - tˇri pro podm´ınˇenou optimalizaci (funkce fmincon) a dvˇe pro nepodm´ınˇenou optimalizaci (funkce fminunc a fminsearch). U funkce fmincon byly pouˇzity algoritmy - active-set, interior-point a sqp. Jednotliv´e algoritmy se liˇs´ı ve zp˚ usobu ˇreˇsen´ı tzv. Kuhn–Tuckerov´ ych podm´ınek. U nepodm´ınˇen´e optimalizace pomoc´ı funkce fminunc byl pouˇzit algoritmus quasi-newton, kter´ y vyuˇz´ıv´a BFGS (Broyden-FletcherGoldfarb-Shanno) algoritmus pro aproximaci Hessovy matice. Nepodm´ınˇen´ a optimalizace pomoc´ı funkce fminsearch vyuˇz´ıv´a nederivaˇcn´ı Nelder˚ uv-Mead˚ uv algoritmus. [9] Vˇsechny algoritmy byly testov´any s p˚ uvodn´ım nastaven´ım, mˇenil se pouze pouˇzit´ y algoritmus. Startovac´ı pozice byla volena n´ ahodnˇe z okol´ı skuteˇcn´eho ˇreˇsen´ı. Pro zv´ yˇsen´ı robustnosti byla optimalizace opakov´ana ze 100 n´ ahodn´ ych pozic a za v´ ysledn´ y odhad byl vybr´an nejlepˇs´ı ze vˇsech odhad˚ u. Celkem tak bylo tedy provedeno 100 simulac´ı pro 100 pr˚ ubˇeh˚ u s pˇeti u ´rovnˇemi ˇsumu pomoc´ı pˇeti metod, celkem tedy 250 000 bˇeh˚ u optimalizace. Pro posouzen´ı spr´avnosti odhadu byla vyhodnocena souhrnn´ a relativn´ı chyba odhadu AIF [4], kter´a je stanovena pro cel´ y pr˚ ubˇeh AIF a vˇsechny realizace ˇsumu na dan´e hladinˇe ve srovn´an´ı s referenˇcn´ım pr˚ ubˇehem AIF v procentech. Pro srovn´an´ı ˇcasov´e n´ aroˇcnosti jednotliv´ ych algoritm˚ u byla vyhodnocena tak´e pr˚ umˇern´ a v´ ypoˇcetn´ı doba pro odhad parametr˚ u u jedn´e realizace ˇsumu.
3
V´ ysledky
V´ ysledky souhrnn´e relativn´ı chyby odhadu AIF jsou uveedeny v Tab. 1 pro vˇsechny kombinace u ´rvon´ı ˇsumu a pouˇzit´ ych algoritm˚ u. Pro velmi n´ızk´e u ´rovnˇe ˇsumu (SNR 40 a 50 dB) byly obdrˇzeny obdobn´e v´ ysledky se souhrnnou relativn´ı chybou odhadu AIF menˇs´ı neˇz 0,5 %. Nejlepˇs´ıch v´ ysledk˚ u pro u ´roveˇ n ˇsumu 50 dB dosahoval algoritmus interior-point, kter´ y je ale z´aroveˇ n i nejpomalejˇs´ım algoritmem. Na stˇredn´ı u ´rovni ˇsumu (SNR 30 dB) bylo dosaˇzeno srovnalteln´ ych v´ ysledk˚ u pomoc´ı algoritm˚ u sqp, active-set a quasi-newton. Pro vyˇsˇs´ı u ´rvovnˇe ˇsumu (20 a 10 dB) poskytuj´ı nelepˇs´ı v´ ysledky algoritmy sqp a quasi-newton, u kter´ ych je souhrnn´ a relativn´ı chyba odhadu AIF menˇs´ı neˇz 7 %. ˇ Casov´ a n´ aroˇcnost byla srovn´av´ana na z´akladˇe pr˚ umˇern´eho v´ ypoˇcetn´ıho ˇcasu na perf´ uzn´ı anal´ yzu jedn´e kˇrivky (tedy jedn´e realizace ˇsumu). Nejlepˇs´ıch v´ ysledk˚ u bylo dosaˇzenopomoc´ı algoritmu active-set ve funkci fmincon a pot´e algoritmy sqp a quasi-newton. Nejhorˇs´ıch v´ ypoˇcetn´ıch ˇcas˚ u bylo dosaˇzeno pomoc´ı algoritmu interior-point, coˇz je d´ ano faktem, ˇze tento algoritmus umoˇzn ˇuje tzv. large-scale podm´ınˇenou optimalizaci. ´ relativn´ı chyba odhadu AIF pro kombinace u ´ rovn´ı ˇ Tabulka 1: Souhrnna sumu a ´ ch algoritm˚ ´ho vy ´ poc ˇetn´ıho c ˇasu v sekunda ´ ch. pouˇ zity u v procentech a pr˚ u mˇ erne SNR [dB] v´ yp. ˇcas 50 40 30 20 10 [s] funkce algoritmus fmincon active-set 0,093 0,27 0,73 2,17 7,1 15 60 fmincon interior-point 0,087 0,27 0,78 2,42 7,64 fmincon sqp 0,15 0,29 0,71 2,12 6,87 27 50 fminsearch nonderivative 0,11 0,27 0,79 2,31 7,52 fminunc quasi-newton 0,26 0,32 0,72 2,19 6,71 31
4
Z´ avˇ er
Pro n´ızk´e u ´rovnˇe ˇsumu (SNR ≥ 30 dB) byla souhrnn´ a relativn´ı chyba odhadu AIF menˇs´ı neˇz 1 % a tedy odhad parametr˚ u modelu velmi pˇresn´ y. Se zvyˇsuj´ıc´ı se hodnotou ˇsumu doch´ azelo k vˇetˇs´ım odchylk´ am mezi referenˇcn´ı a odhadnutou kˇrivkou AIF a t´ım i n´ arustu vyhodnocovan´e chyby. Pro perf´ uzn´ı anal´ yzu metodou bolus & burst je ˇsum v datech kritick´ y pro spr´avn´ y odhad perf´ uzn´ıch parametr˚ u. Jak uk´azala prvotn´ı studie [6] na preklinick´ ych datech, dosaˇziteln´ a hodnota ˇsumu je mezi 10 a 20 dB. Pro tyto hodnoty se jako nejvhodnˇejˇs´ı ukazuje algoritmus sqp ze skupiny metod pro podm´ınˇenou optimalizaci, popˇr. algoritmus quasi-newton pro nepodm´ınˇenou optimalizaci. Obˇe metody maj´ı tak´e pˇrijateln´ y v´ ypoˇcetn´ı ˇcas.
Reference ˇ ´IK, Radovan, Kim NYLUND, Odd Helge GILJA, Martin MEZL, ´ ˇ Radim [1] JIR Vratislav HARABIS, ´ R, ˇ Michal STANDARA a Torfinn TAXT. Ultrasound perfusion analysis combining bolusKOLA tracking and burst-replenishment. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control. 2013, 60(2): 310-319. ˇ ´IK, Radovan, Kim NYLUND, Torfinn TAXT, Martin MEZL, ´ [2] JIR Trygve HAUSKEN, Vratislav ˇ ´ ˇ HARABIS, Radim KOLAR, Michal STANDARA a Odd Helge GILJA. Parametric ultrasound perfusion analysis combining bolus tracking and replenishment. 2012 IEEE International Ultrasonics Symposium. IEEE, 2012, : 1323-1326. ˇ ´IK, Radovan, Karel SOUCEK, ˇ ´ ˇ Eva DRAZANOV ˇ ´ Frantisek [3] JIR Martin MEZL, Michal BARTOS, A, ´ ´ ´ ´ ˇ DRAFI, Lucie GROSSOVA, Jiˇr´ı KRATOCHVILA, Ondˇrej MACICEK, et al. Blind deconvolution in dynamic contrast-enhanced MRI and ultrasound. 2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE, 2014, 60(2): 4276-4279.
ˇ ´IK, Michal BARTOS, ˇ Michal STANDARA, Zenon STARCUK ˇ [4] KRATOCHV´ILA, Jiˇr´ı, Radovan JIR a Torfinn TAXT. Distributed capillary adiabatic tissue homogeneity model in parametric multichannel blind AIF estimation using DCE-MRI. Magnetic Resonance in Medicine. 2015, : n/a-n/a. ´ ˇ ´IK, Vratislav HARABIS, ˇ Radim KOLA ´ R, ˇ Michal STANDARA, Kim [5] MEZL, Martin, Radovan JIR NYLUND, Odd Helge GILJA, Torfinn TAXT, Odd Helge GILJA, et al. Absolute ultrasound perfusion parameter quantification of a tissue-mimicking phantom using bolus tracking [Correspondence]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control. IEEE, 2015, 62(5): 983-987. ´ ˇ ´IK, Karel SOUCEK ˇ ´ R. ˇ Evaluation of Accuracy of [6] MEZL, Martin, Radovan JIR a Radim KOLA Bolus and Burst Method for Quantitative Ultrasound Perfusion Analysis with Various Arterial Input Function Models. 2015 IEEE International Ultrasound Symposium. 2015, (in press). ˇ ´IK, Martin MEZL, ´ [7] NYLUND, Kim, Radovan JIR Sabine LEH, Trygve HAUSKEN, Frank PFEFFER, Svein ØDEGAARD, Torfinn TAXT, Odd Helge GILJA, et al. Quantitative Contrast-Enhanced Ultrasound Comparison Between Inflammatory and Fibrotic Lesions in Patients with Crohn’s Disease. Ultrasound in Medicine. IEEE, 2013, 39(7): 1197-1206. [8] STROUTHOS, Costas, Marios LAMPASKIS, Vassilis SBOROS, Alan MCNEILLY, Michalakis AVERKIOU, Kim NYLUND, Odd Helge GILJA, Torfinn TAXT, Odd Helge GILJA, et al. Indicator dilution models for the quantification of microvascular blood flow with bolus administration of ultrasound contrast agents. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control. IEEE, 2010, 57(6): 1296-1310. [9] VENKATARAMAN, P. Applied optimization with MATLAB programming. 2nd ed. Hoboken, N.J.: John Wiley, 2009, xvi, 526 p. ISBN 04-700-8488-X. [10] WEI, K., A. R. JAYAWEERA, S. FIROOZAN, A. LINKA, D. M. SKYBA, S. KAUL, Odd Helge GILJA, Torfinn TAXT, Odd Helge GILJA, et al. Quantification of Myocardial Blood Flow With Ultrasound-Induced Destruction of Microbubbles Administered as a Constant Venous Infusion. Circulation. IEEE, 1998, 97(5): 473-483.
Martin M´ezl ´ Ustav biomedic´ınsk´eho inˇzen´ yrstv´ı, Vysok´e uˇcen´ı technick´e v Brnˇe. e-mail:
[email protected] Radovan Jiˇr´ık ´ ˇ v.v.i.. Ustav pˇr´ıstrojov´e techniky AV CR, e-mail:
[email protected]
