HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele...

21
HEART :: High PErformance Computing of PersonAlized CaRdio ComponenT Models Newsletter Numărul 2 Octombrie 2016 http://heart.unitbv.ro/ HEART High PErformance Computing of PersonAlized CaRdio ComponenT Models Newsletter Numărul 2 Octombrie 2016 În cadrul acestui newsletter se prezintă o serie de activități care au fost desfășurate în vederea atingerii obiectivelor formulate la începutul proiectului. Astfel, s-a continuat dezvoltarea platformei utilizate pentru reconstrucţia tridimensională a arterelor coronariene, introducându-se posibilitatea prelucrării imaginilor angiografice care nu au EKG ataşat. În ceea ce priveşte modelarea multiscalară a sistemului cardiovascular şi a circulaţiei coronariene, activităţile au fost împărţite în trei mari categorii. În primul rând s-a dezvoltat un nou algoritm de estimare a debitului coronarian. Acesta ia în considerare toate secţiunile transversale sănătoase ale tuturor ramurilor modelului anatomic reconstruit. Aspectele cheie sunt: folosirea unei metodologii locală-globală-locală care asigură că principiul de conservare a masei este respectat şi definirea unor coeficienţi de încredere atât pentru ramuri cât şi pentru generaţiile de ramuri. În al doilea rând, s-a dezvoltat un framework ierarhic de personalizare a condiţiilor de frontieră în simulări hemodinamice care folosesc arbori structuraţi ca şi condiţii de frontieră de ieşire. La primul nivel se estimează proprietăţi hemodinamice precum rezistenţa şi complianţa, iar la cel de-al doilea nivel se estimează parametrii arborilor structuraţi astfel încât proprietăţile hemodinamice estimate la primul nivel să fie respectate. Un aspect cheie al abordării propuse este faptul că pentru fiecare proprietate hemodinamică se estimează valorile a doi parametrii diferiţi ai arborilor structuraţi. Framework-ul propus a fost validat cu succes folosind un model aortic cu coarctaţie. În al treilea rând, s-a dezvoltat o metodologie pentru simularea interacţiunii fluid-solid, bazată pe metoda Lattice Boltzmann. În particular, s-a dezvoltat o abordare eficientă pentru înglobarea geometriei variabile în timp. Metodologia a fost validată printr-o serie de experimente sintetice: un cilindru care se expandează şi se contractă periodic şi curgerea peristaltică tridimensională (deformarea periodică a peretelui conduce la curgerea fluidului din interiorul vasului). Soluţiile numerice obţinute au fost de fiecare dată foarte apropiate de cele analitice. O altă activitate importantă a fost achiziţia unui volum mare de date clinice. Astfel, la Spitalul Clinic de Urgenţă București au fost incluşi până în prezent 48 de pacienţi în studiu şi s-au realizat în total 120 de măsurători intracoronariene, dintre care 92 înainte de stentare şi 28 după stentare. În continuare s-au dezvoltat algoritmi pentru postprocesarea datelor achiziţionate în vederea determinării indicatorilor diagnostici IFR şi FFR. Folosind platforma menţionată mai sus s-a realizat reconstrucţia tridimensională a arterelor coronariene, s-au rulat simulările hemodinamice multiscalare şi s-au estimat indicatorii IFR şi FFR prin 1. Reconstrucția tridimensională a arterelor coronariene ........................................................ 2 2. Modelarea multiscalară a sistemului cardiovascular .................................................... 2 3. Achiziția datelor de la pacienți ......................... 11 4. Procesarea datelor achiziționate și validarea modelului multiscalar al circulației coronariene ......................................................................... 11 5. Model cu parametrii distribuiți al sistemului cardiovascular .................................................. 13 6. Utilizarea dimensiunii fractale în identificarea stenozelor coronariene .................................... 16 7. Utilizarea procesării paralele în accelerarea timpilor de execuție al algoritmilor dezvoltați . 18

Transcript of HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele...

Page 1: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART :: High PErformance Computing of PersonAlized CaRdio ComponenT Models

Newsletter

Numărul 2

Octombrie 2016

http://heart.unitbv.ro/

HEART

High PErformance Computing of

PersonAlized CaRdio ComponenT Models

Newsletter

Numărul 2

Octombrie 2016

În cadrul acestui newsletter se prezintă o serie de activități care au fost desfășurate în

vederea atingerii obiectivelor formulate la începutul proiectului. Astfel, s-a continuat

dezvoltarea platformei utilizate pentru reconstrucţia tridimensională a arterelor

coronariene, introducându-se posibilitatea prelucrării imaginilor angiografice care nu au

EKG ataşat.

În ceea ce priveşte modelarea multiscalară a sistemului cardiovascular şi a circulaţiei

coronariene, activităţile au fost împărţite în trei mari categorii. În primul rând s-a

dezvoltat un nou algoritm de estimare a debitului coronarian. Acesta ia în considerare

toate secţiunile transversale sănătoase ale tuturor ramurilor modelului anatomic

reconstruit. Aspectele cheie sunt: folosirea unei metodologii locală-globală-locală care

asigură că principiul de conservare a masei este respectat şi definirea unor coeficienţi de

încredere atât pentru ramuri cât şi pentru generaţiile de ramuri. În al doilea rând, s-a

dezvoltat un framework ierarhic de personalizare a condiţiilor de frontieră în simulări

hemodinamice care folosesc arbori structuraţi ca şi condiţii de frontieră de ieşire. La

primul nivel se estimează proprietăţi hemodinamice precum rezistenţa şi complianţa, iar

la cel de-al doilea nivel se estimează parametrii arborilor structuraţi astfel încât

proprietăţile hemodinamice estimate la primul nivel să fie respectate. Un aspect cheie al

abordării propuse este faptul că pentru fiecare proprietate hemodinamică se estimează

valorile a doi parametrii diferiţi ai arborilor structuraţi. Framework-ul propus a fost

validat cu succes folosind un model aortic cu coarctaţie. În al treilea rând, s-a dezvoltat o

metodologie pentru simularea interacţiunii fluid-solid, bazată pe metoda Lattice

Boltzmann. În particular, s-a dezvoltat o abordare eficientă pentru înglobarea geometriei

variabile în timp. Metodologia a fost validată printr-o serie de experimente sintetice: un

cilindru care se expandează şi se contractă periodic şi curgerea peristaltică

tridimensională (deformarea periodică a peretelui conduce la curgerea fluidului din

interiorul vasului). Soluţiile numerice obţinute au fost de fiecare dată foarte apropiate de

cele analitice.

O altă activitate importantă a fost achiziţia unui volum mare de date clinice. Astfel, la

Spitalul Clinic de Urgenţă București au fost incluşi până în prezent 48 de pacienţi în

studiu şi s-au realizat în total 120 de măsurători intracoronariene, dintre care 92 înainte

de stentare şi 28 după stentare.

În continuare s-au dezvoltat algoritmi pentru postprocesarea datelor achiziţionate în

vederea determinării indicatorilor diagnostici IFR şi FFR. Folosind platforma menţionată

mai sus s-a realizat reconstrucţia tridimensională a arterelor coronariene, s-au rulat

simulările hemodinamice multiscalare şi s-au estimat indicatorii IFR şi FFR prin

1.

Reconstrucția tridimensională a arterelor

coronariene ........................................................ 2

2.

Modelarea multiscalară a sistemului

cardiovascular .................................................... 2

3.

Achiziția datelor de la pacienți ......................... 11

4.

Procesarea datelor achiziționate și validarea

modelului multiscalar al circulației coronariene

......................................................................... 11

5.

Model cu parametrii distribuiți al sistemului

cardiovascular .................................................. 13

6.

Utilizarea dimensiunii fractale în identificarea

stenozelor coronariene .................................... 16

7.

Utilizarea procesării paralele în accelerarea

timpilor de execuție al algoritmilor dezvoltați . 18

Page 2: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

2

http://heart.unitbv.ro/

intermediul simulării.

Acurateţea clasificării realizate de mărimile estimate din simulări a fost de 74.4% şi

respectiv 86.5% pentru cei doi indicatori diagnostici. De asemenea, s-a dezvoltat o

strategie hibridă de diagnosticare, în urma căreia doar 34% dintre leziuni ar necesita o

măsurătoare invazivă (66% dintre leziuni ar fi clasificate pe baza valorilor estimate din

simulare). Acurateţea clasificării acestei strategii este de 96.5%.

De asemenea, s-a continuat şi dezvoltarea modelului cu parametrii distribuiţi al

întregului circuit cardiovascular. S-a dezvoltat o nouă strategie de personalizare bazată

pe 20 de obiective de personalizare formulate pe baza presiunilor ventriculare şi aortice

şi pe baza volumelor ventriculare.

O altă activitatea s-a concentrat pe utilizarea dimensiunii fractale în identificarea

stenozelor coronariene. Metoda a fost testată cu succes la şase pacienţi, demonstrându-

se că prin utilizarea dimensiunii fractale se pot identifica cu acurateţe stenozele

coronariene.

Pentru a reduce timpii de execuție necesari pentru rularea simulărilor hemodinamice s-

au desfăşurat două activităţi. În primul rând s-a realizat o implementare multicore a

modelului utilizat pentru simulările coronariene. Ca urmare, timpul de execuţie mediu

pentru un model anatomic a fost redus de la 230.9 secunde la 38.1 secunde. În al doilea

rând, s-a realizat o implementare pe GPU a algoritmului multigrid geometric pentru

rezolvarea sistemelor de ecuaţii liniare sparse. Acesta a fost comparat cu un algoritm

PCG optimizat şi s-a obţinut un timp de execuţie de 7.1 – 9.2 ori mai mic.

1.

Reconstrucția tridimensională a arterelor coronariene

Platforma utilizată pentru reconstrucţia tridimensională a vaselor,

care a fost descrisă în newsletter-ul anterior a fost dezvoltată în

continuare. Una din modificările importante a fost introducerea

posibilităţii prelucrării imaginilor angiografice care nu au EKG ataşat.

EKG-ul este folosit pentru a identifica finalul diastolei (reconstrucţia

geometriei tridimensionale foloseşte cadre de la finalul diastolei

deoarece mişcarea miocardului este redusă în această fază a ciclului

cardiac).

Deoarece însă datele achiziţionate de la o serie de pacienţi nu

conţineau şi EKG-ul, algoritmul de reconstrucţie a fost adaptat astfel

încât să se poată folosi orice cadru angiografic.

Alte îmbunătăţiri aduse platformei sunt:

editarea mai simplă şi mai rapidă a segmentării;

posibilitatea de a marca startul şi finalul stenozei de către

utilizator.

2.

Modelarea multiscalară a sistemului cardiovascular

2.1. Algoritm de estimare a debitului coronarian

S-a dezvoltat un nou algoritm de estimare a debitului coronarian

pornind de la geometria coronariană. Acest algoritm este format din

trei paşi principali şi foloseşte o metodologie locală-globală-locală,

prezentată în figura 2.1. Iniţial se estimează pentru fiecare ramură o

valoare a debitului (pasul 1, local). Deoarece toate debitele ramurilor

sunt calculate independent condiţia conform căreia debitul unei

ramuri părinte trebuie să fie egal cu suma debitelor ramurilor copil

nu este îndeplinită. De aceea, se estimează o singură valoare

(globală) pentru întreg arborele coronarian prin medierea valorilor

ramurilor de la diferite generaţii (pasul 2, global). În final, debitul

total este distribuit către ramurile individuale astfel încât să fie

satisfăcută condiţia descrisă mai sus (pasul 3, local).

La primul pas (local) debitul fiecărei ramuri este estimat

independent. Deoarece raza variază de-a lungul axei centrale a unui

segment, se realizează o mediere a tuturor razelor care aparţin unor

segmente coronariene sănătoase (fără stenoze).

La cel de-al doilea pas se calculează un debit total pentru întreg

arborele coronarian (stâng sau drept) pe baza debitelor estimate la

primul pas pentru fiecare ramură. În acest sens se va calcula mai întâi

o valoare globală pentru fiecare generaţie. Figura 2.2 prezintă un

arbore coronarian, în care fiecărei ramuri i-a fost asociat un număr

de generaţie. Ramura rădăcină are numărul de generaţie egal cu 0,

care apoi creşte cu unu la fiecare bifurcaţie. În continuare se va face

referire la o generaţie prin indexul g.

Page 3: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

3

http://heart.unitbv.ro/

Figura 2.1. Schema bloc a algoritmului de estimare a debitelor segmentelor

coronariene.

Figura 2.2. Exemplu de arbore coronarian în care fiecărei ramuri i-a fost

asociat un număr de generaţie.

Înainte de a estima debitul global, fiecărei ramuri îi este asociat şi un

coeficient care reflectă încrederea în corectitudinea valorii debitului

estimat pentru acea ramură. Ramuri foarte scurte sau ramuri

complet stenozate vor avea valori de încredere mici, în timp ce

ramuri lungi fără variaţii mari ale razelor vor avea valori mari de

încredere.

În continuare valoarea finală a debitului total este calculată din

valorile individuale estimate separat pentru fiecare generaţie. Şi de

această dată se foloseşte un coeficient de încredere care este ataşat

fiecărei generaţii.

La cel de-al treilea pas (local), se estimează valoarea finală a debitului

fiecărei ramuri. În final, debitele celorlalte ramuri sunt calculate ca

sumă a debitelor ramurilor terminale aflate în aval.

2.2. Dezvoltarea unui framework ierarhic de personalizare a

condiţiilor de frontieră în simulări hemodinamice

Sistemul cardiovascular este compus din aproximativ 10 miliarde de

vase, ale căror dimensiuni variază într-un interval care acoperă mai

multe ordine de mărime. Prin urmare modelarea 3D sau chiar de

ordin redus (1D) a întregului sistem nu este fezabilă din punctul de

vedere al resurselor de calcul necesare. De aceea, doar zona de

interes este modelată spaţial, în timp ce restul componentelor sunt

modelate cu parametri distribuiţi, reprezentând condiţii de frontieră

artificiale pentru zona de interes. În funcţie de disponibilitatea

măsurătorilor in-vivo şi de ipotezele modelului, cercetările anterioare

au folosit unul din următoarele tipuri de condiţii de frontieră: viteză

(debit) variabilă în timp [Olufsen et al., 2000], [LaDisa et al., 2011],

sau un model cu parametri distribuiţi al inimii [Formaggia et al.,

2006], [Coogan et al., 2011]. Impunerea condiţiilor de frontieră de

ieşire este mai dificilă, deoarece:

vasculatura distală (microvasculatura) generează cea mai mare

parte a rezistenţei şi prin urmare este responsabilă de distribuţia

debitului şi de nivelul global de presiune din zona de interes;

undele de debit şi presiune se propagă dincolo de frontierele de

ieşire. Atunci când vasele îşi schimbă geometria şi proprietăţile,

sau se bifurcă, undele se modifică şi ele la rândul lor.

Suplimentar, undele sunt reflectate şi se propagă în sens invers,

revenind în zona de interes.

Diferite abordări au fost propuse în literatura de specialitate pentru

impunerea condiţiilor de frontieră de ieşire, acestea acoperind o

gamă largă de la condiţii de tip presiune sau debit până la modele cu

parametrii distribuiţi. S-a concluzionat că pentru a realiza simulări

personalizate trebuie folosite condiţii fiziologice. Prin cuplarea

acestor condiţii de frontieră cu zona de interes se obţin modele

geometrice multiscalare. Condiţiile de frontieră de ieşire sunt în

general dezvoltate cu scopul suprinderii unuia sau mai multor efecte

dintre următoarele: (i) rezistenţa totală, (ii) complianţa totală şi (iii)

efectele propagării şi reflexiei undelor în vasculatura distală. Cea mai

des utilizată condiţie de frontieră este modelul windkessel cu 3

elemente, care este simplu (are doar 3 parametrii) şi este capabil să

suprindă primele două dintre cele trei efecte menţionate mai sus.

Numărul mic de parametrii simplifică personalizarea modelului, dar

dezavantajul este că nu poate surprinde fenomenele de propagare şi

reflexie a undelor hemodinamice din vasculatura distală.

O altă condiţie de frontieră, dezvoltată special pentru a surprinde

fenomenele de propagare, este modelul bazat pe un arbore

structurat. A fost iniţial prezentat în [Taylor, 1966] şi apoi dezvoltat

de Olufsen et al. [Olufsen et al., 1999]. Vasculatura distală este

modelată ca o structură geometrică simplă iar în urma aplicării unor

ipoteze simplificatoare se obţine o expresie analitică bazată pe debit

şi presiune. Datorită caracteristicilor sale, modelele bazate pe arbori

structuraţi surprind toate cele trei efecte care sunt de interes pentru

o condiţie de frontieră de ieşire. Pe de altă parte, deoarece rezistenţa

şi complianţa nu sunt parametrizate în mod explicit în expresia

analitică a arborelui structurat, este mult mai dificil să se dezvolte

proceduri de estimare a parametrilor pentru simulările

hemodinamice care folosesc astfel de condiţii de frontieră. Din acest

motiv, arborele structurat a fost folosit mai rar pentru simulări

hemodinamice personalizate.

Cu toate acestea, condiţia de frontieră de tip arbore structurat a fost

studiată intens în ultimii ani. Astfel, Cousins et al., au introdus o

metodologie mai simplă de obţinere a expresiei analitice şi au realizat

o analiză de sensibilitate în raport cu parametrii arborelui structurat

[Cousins et al., 2012]. Într-un alt studiu, s-a introdus o formulare

diferită care poate fi utilizată pentru a modela curgeri tranzitorii

[Cousins et al., 2013]. Recent, arborele structurat a fost utilizat cu

succes în simulările hemodinamice ale arterelor şi venelor retinei

[Malek et al., 2015].

Page 4: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

4

http://heart.unitbv.ro/

Diferite proceduri de estimare a parametrilor au fost propuse în

ultimii ani [Spilker et al., 2010], [Itu et al., 2015], [Ismail et al., 2015],

[Blanco et al., 2012], [Bertoglio et al., 2012], în special pentru

condiţiile de frontieră de tip windekessel. Recent, s-a introdus şi o

metodologie mai riguroasă de calibrare a parametrilor arborelui

stucturat [Cousins et al., 2014], unde metodă de tip trust-region a

fost aplicată pentru a adapta raportul lungime-rază, astfel încât să se

obţină un rezultat hemodinamic care să corespundă distribuţiei

măsurate a debitelor.

În continuare se va prezenta un framework ierarhic iterativ de

estimare a parametrilor, dezvoltat pentru a fi utilizat în simulări

hemodinamice care aplică modele bazate pe arbori structuraţi la

frontieră de ieşire. Parametrii arborilor structuraţi sunt adaptaţi

astfel încât să minimizeze erorile dintre valorile simulate şi măsurate

de debit şi presiune. Se foloseşte o metodă ierarhică deoarece

proprietăţile vaselor (rezistenţa şi complianţa) nu sunt parametrizate

în mod explicit în expresia matematică a arborelui structurat. La

primul nivel se estimzează rezistenţele şi complianţele astfel încât

rezultatle simulării să corespundă măsurătorilor efectuate, iar la cel

de-al doilea nivel se estimează parametrii arborilor structuraţi pentru

a obţine valorile de rezistenţă şi complianţă estimate la primul nivel.

Iniţial s-a evaluat cel de-al doilea nivel pe baza valorilor parametrilor

modelelor windkessel utilizate ca şi condiţii de frontieră într-un

model al circulaţiei sistemice arteriale [Stergiopulos et al., 1992].

Apoi, geometria unui pacient cu coarctaţie aortică este utilizată

pentru a testa şi valida întreg framework-ul. Rezultatele sunt

comparate cu o configuraţie în care se folosesc modele windkessel ca

şi condiţii de frontieră.

2.2.2. Metode

2.2.2.1. Condiţia de frontieră de tip arbore structurat

În figura 2.3 se prezintă structura generică a unui arbore structurat.

În continuare, se prezintă doar o scurtă descriere a principalelor sale

aspecte.

Figura 2.3. Arbore structurat asimetric [Olufsen, 1998].

Arborele structurat este un arbore asimetric, binar, în interiorul

căruia fiecare vas este simetric în raport cu axa centrală şi are o rază

constantă.

În cazul condiţiei de frontieră de tip arbore structurat, bifurcaţiile

sunt considerate a fi asimetrice. Prin urmare, razele celor două vase

copil pot fi calculate pe baza razei vasului părinte și a unui coeficient

de asimetrie:

Pornind de la vasul rădăcină, cu o anumită rază, arborele structurat

se bifurcă până când raza vaselor devine mai mică decât o anumită

rază minimă. Un alt parametru important al arborelui structurat

determină lungimea vaselor. Pe baza rezultatelor obţinute în studii

anterioare [Iberall, 1967], lungimea poate fi exprimată prin

intermediul razei fiecărui vas, adoptând un raport lungime-rază.

În final, trebuie specificate şi proprietăţile pereţilor vaselor de sânge

care compun arborele structurat. După cum a fost prezentat în

literatură [Olufsen et al., 2000], deoarece arterele mici sunt compuse

din aceleaşi materiale ca şi arterele mari, expresia dedusă pentru

arterele mari poate fi aplicată și pentru arterele mici [Olufsen et al.,

2000].

Ecuaţiile care modelează curgerea sângelui în arborele structurat

sunt obţinute din ecuaţiile Navier-Stokes, cu supoziţia de simetrie

axială [Olufsen, 1998]. Deoarece în arterele mici efectul vâscozităţii

este mult mai important decât efectul inerţiei, termenii neliniari,

inerţiali, pot fi neglijaţi. Dacă debitul şi presiunea sunt periodice, se

poate determina o soluţie analitică în domeniul frecvenţei.

Impedanţa la rădăcina arborelui structurat este calculată recursiv.

Impedanţa obţinută la rădăcina arborelui este apoi aplicată ca şi

condiţie de frontieră de ieşire.

Deoarece condiţia de frontieră la intrarea arborelui structurat este

periodică, presiunea şi debitul pot fi exprimate prin serii Fourier

complexe. Prin această abordare caracteristicile răspunsului

sistemului pot fi determinate separat pentru fiecare termen.

În cadrul acestei activităţi arborele structurat a fost folosit ca şi

condiţie de frontieră periodică. Pentru a aplica o condiţie de frontieră

Page 5: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

5

http://heart.unitbv.ro/

periodică, istoricul debitului de-a lungul ultimei perioade trebuie

memorat şi trebuie realizată o operaţie de scanare de înmulţire-

adunare la fiecare moment de timp discret. Acest aspect conduce la

timpi de execuţie mult mai mari decât pentru o condiţie de frontieră

aperiodică.

Se subliniază faptul că modul de aplicare a celor două tipuri de

condiţii de frontieră nu este impus prin natura acestora. Modelul

windkessel poate fi aplicat ca şi condiţie de frontieră periodică şi,

după cum a fost descris recent în lucrarea [Cousins et al., 2012], şi

arborele structurat poate fi aplicat ca şi condiţie de frontieră

aperiodică.

2.2.2.2. Framework-ul de estimare a parametrilor

În continuare este prezentat framework-ul ierarhic de estimare a

parametrilor simulărilor hemodinamice multiscalare care utilizează

arbori structuraţi ca şi condiţii de frontieră de ieşire (figura 2.4).

Pentru a evalua framework-ul se foloseşte un model multiscalar de

ordin redus care combină un model cvasi unu-dimensional cu

modelul arborelui structurat. Ca şi condiţii de frontieră de intrare se

utilizează debite variabile în timp. Ambele nivele ale framework-ului

implică metode de calibrare complet automatizate pentru estimarea

valorilor parametrilor.

Metoda de calibrare a rezistenţelor şi complianţelor

Metoda de calibrare utilizată la primul nivel (figura 2.4a) estimează în

mod automat rezistenţa totală şi complianţa arborilor structuraţi

astfel încât diferenţele dintre valorile măsurate şi simulate de

presiune şi debit să fie cât mai mici. Procedura de estimare utilizată

la primul nivel este o versiune modificată a unui framework introdus

anterior [Itu et al., 2015]. Diferenţa apare la pasul 9, unde este apelat

cel de-al doilea nivel al framework-ului pentru adapta parametrii

arborilor structuraţi.

Cel de-al doilea nivel al framework-ului (figura 2.4b) estimează

valorile parametrilor arborilor structuraţi astfel încât să se obţină

rezistenţele şi complianțele determinate la primul nivel (la pasul 8 în

figura 2.4a). Trei abordări diferite au fost propuse în trecut pentru

adaptarea rezistenţei totale a arborelui structurat:

impunerea unei rezistenţele la ieşirea fiecărui vas terminal a

arborelui structurat [Olufsen, 1998];

adaptarea razei minime, până la care se generează arborele

structurat [Cousins et al., 2012];

adaptarea raportului lungime-rază, care determină lungimea

fiecărui vas din arborele structurat [Cousins et al., 2014].

Figura 2.4: Framework ierarhic de estimare a parametrilor utilizat în simulări hemodinamice care folosesc arbori structuraţi la frontierele de ieşire: (a) Nivelul 1 –

Metoda de calibrare care estimează rezistenţele totale şi complianţa arborilor structuraţi; (b) Nivelul 2 – Metoda de calibrare care estimează valorile parametrilor

arborilor structuraţi astfel încât să se obţină rezistenţele şi complianţele totale estimate la nivelul 1.

Page 6: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

6

http://heart.unitbv.ro/

Metoda de calibrare a parametrilor arborilor structuraţi

Pentru a putea simula o gamă largă de stări pentru acelaşi pacient

(de exemplu: repaus, exerciţiu uşor, exerciţiu intens, etc.) rezistenţa

trebuie să poată fi variată într-un interval mare de valori. Prin

urmare, în acest studiu s-a folosit o combinaţie între raza minimă şi

rezistenţa impusă la vasele terminale ale arborelui structurat pentru

a adapta rezistenţa totală a arborelui structurat. Deoarece vasele

terminale au aproximativ aceeaşi rază, se presupune că rezistenţele

terminale sunt egale. Raza minimă se foloseşte pentru o calibrare

grosieră, iar rezistenţa terminală pentru calibrarea fină.

În ceea ce priveşte complianţa totală, în literatura de specialitate nu

au fost raportate metode pentru a realiza o adaptare a acestei

caracteristici. Pentru ca valoarea complianţei să poată fi variată într-o

gamă largă de valori corespunzătoare mai multor stări ale

pacientului, parametrii modelului arterial sunt adaptaţi pentru

calibrarea grosieră și pentru calibrarea fină.

Metoda de calibrare fină este similară cu cea de la primul nivel:

parametrii arborelui structurat sunt determinaţi ca soluţie a unui

sistem de ecuaţii neliniare, care are o rădăcină în punctul în care

proprietăţile calculate şi de referinţă ale arborelui structurat sunt

egale.

Deoarece rezistenţele şi complianţele calculate pot fi determinate

rapid, iacobianul este recalculat la fiecare iteraţie. Prin urmare, la

acest nivel se foloseşte o metoda dogleg trust region în locul unei

metode cvasi-Newton.

Complianţa volumetrică este calculată analitic prin adunarea

complianţelor volumetrice ale tuturor vaselor din arborele structurat.

2.2.3. Rezultate

Pentru a evalua performanţa framework-ului propus, s-a folosit un

model anatomic de coarctaţie reconstruit din imagini MRI.

Framework-ul de calibrare asigură faptul că simularea este

personalizată şi că, prin urmare, valorile de debit şi presiune calculate

sunt apropiate de măsurătorile clinice. Atât arbori structuraţi cât şi

modele windkessel au fost folosite pentru a impune condiţiile de

frontieră de ieşire, şi rezultatele au fost comparate pentru a sublinia

avantajele furnizate de condiţia de frontieră de tip arbore structurat.

Iniţial însă sunt prezentate rezultatele obţinute în cadrul unui test în

care nivelul al doilea al framework-ului de estimare a fost evaluat

separat pe baza condiţiilor de frontieră de ieşire ale unui model

sistemic arterial al întregului corp.

2.2.3.1. Model sistemic arterial

S-a folosit modelul arterial propus anterior în [Stergiopulos et al.,

1992]. Astfel, rezistenţele şi complianţele totale ale condiţiilor de

frontieră de ieşire au fost folosite ca valori de referinţă, iar metoda

de calibrare prezentată în figura 2.4b a fost rulată separat pentru

fiecare frontieră.

Metoda de calibrare a convers pentru fiecare arbore structurat (doar

2-10 iteraţii au fost necesare). Raza minimă este mai mare de 0.05cm

pentru anumite artere: în aceste cazuri rezistenţa totală iniţială a fost

mai mare decât rezistenţa de referinţă, chiar şi atunci când rezistenţa

terminală a fost setată la valoarea 0. În general, rezistenţele

terminale impuse la finalul arborilor structuraţi sunt cu trei până la

cinci ordine de magnitudine mai mari decât rezistenţa totală a

arborelui structurat.

Metoda de calibrare este eficientă din punct de vedere al timpului de

calcul: deoarece numărul de iteraţii este mic şi atât rezistenţa totală

cât şi complianţa unui arbore structurat sunt determinate analitic,

timpul de execuţie necesar pentru a calibra un arbore structurat este

mai mic de 3 secunde pe un procesor Intel i7 CPU cu 3.4 GHz.

Trebuie remarcat faptul că atât rezistenţele cât şi complianţele de

referinţă iau valori fiziologice. Este posibil ca sistemul de ecuaţii

neliniare să nu aibă o soluţie dacă una din valorile de referinţă se află

înafara intervalului valorilor fiziologice obţinute.

2.2.3.2. Simulare hemodinamică personalizată a unui model

anatomic cu coarctaţie aortică

Metode bazate pe dinamica fluidelor au fost propuse în trecut pentru

evaluarea non-invazivă a căderii de presiune de-a lungul stenozelor

[Itu et al., 2012], [Keshavarz-Motamed et al., 2011], [LaDisa et al.,

2011], [Ismail et al., 2013]. Pentru a estima corect căderea de

presiune, rezultatele modelului hemodinamic trebuie să corespundă

valorilor măsurate de presiune şi debit.

Modelul anatomic utilizat în cadrul acestei activităţi [***CFD

challenge, 2014] conţine aorta ascendentă, trei ramuri supra-aortice,

arcul aortic şi aorta descendentă cu coarctaţie (Figura 2.5a). Figura

2.5b prezintă modelul multiscalar corespunzător.

Debitul măsurat în aorta ascendentă este utilizat ca şi condiţie de

frontieră de intrare în timp ce rapoartele medii de divizare a

debitelor la frontierele de ieşire sunt cunoscute. Scopul final este de

a determina căderea de presiune de-a lungul coarctaţiei. Obiectivele

de la primul nivel al framework-ului de estimare a parametrilor sunt

formulate pe baza rapoartelor de divizare a debitelor şi a presiunilor

sistolice şi diastolice. Modelul vasculaturii distale utilizat la primul pas

al procedurii de calibrare pentru a găsi o soluţie iniţială apropiata de

cea finală este afişat în figura 2.5c.

Pentru a crea modelul multiscalar (pasul 5 din figura 2.4) s-a folosit

biblioteca vmtk [***vmtk, 2014] în vederea determinării axei

centrale şi ariilor secţiunilor transversale. În continuare, pentru

fiecare ramură a modelului anatomic s-au folosit mai multe

segmente unu-dimensionale cu arie transversală care variază de-a

lungul direcţiei longitudinale.

Parametrii estimaţi la primul nivel al framework-ului de estimare a

parametrilor sunt rezistenţele totale ale vaselor supra-aortice şi ale

aortei descendente şi complianţa totală.

Sistemul de ecuaţii neliniare a fost rezolvat în două configuraţii:

configuraţia din figura 2.5, în care arborii structuraţi sunt folosiţi

ca şi condiţii de frontieră de ieşire şi framework-ul de estimare a

parametrilor din figura 2.4 a fost aplicat pentru personalizarea

parametrilor;

o configuraţie în care modelele windkessel cu trei elemente sunt

folosite ca şi condiţii de frontieră de ieşire şi doar primul nivel al

framework-ului de estimare a parametrilor a fost aplicat pentru

personalizarea parametrilor (pasul 9 este eliminat). Rezistenţa

proximă a fiecărui model windkessel este setată la o valoare

egală cu rezistenţa caracteristică şi este menţinută constantă de-

a lungul procedurii de calibrare a parametrilor.

Page 7: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

7

http://heart.unitbv.ro/

Figura 2.5. (a) Geometria aortei proxime cu coarctaţie; (b) Model multiscalar utilizat pentru determinarea valorilor parametrilor arborilor structuraţi specifici

pacientului; (c) Model cu parametrii distribuiţi utilizat la primul pas al algoritmului de personalizare pentru găsirea unei soluţii iniţiale a problemei de calibrare.

Deoarece doar complianţa totală este utilizată ca parametru în

sistemul de ecuaţii neliniare, în ambele configuraţii aceasta este

distribuită către modelele utilizate ca şi condiţii de frontieră de ieşire.

Un aspect important al simulărilor hemodinamice cu pereţi

complianţi este estimarea proprietăţilor mecanice ale peretelui

arterial. Pentru a calcula proprietăţile segmentelor aortice se

foloseşte o metodă bazată pe estimarea vitezei de undă [Olufsen et

al., 2000]. Pentru a estima proprietăţile vaselor supra-aortice, se

foloseşte o metodă diferită în baza căreia proprietăţile fiecărei artere

supra-aortice sunt estimate separat [Itu et al., 2012].

În continuare se prezintă rezultatele procedurii de calibrare obţinute

prin aplicarea framework-ului ierarhic şi arbori structuraţi ca şi

condiţii de frontieră în modelul multiscalar. Valorile parametrilor şi

obiectivelor sunt prezentate în figura 2.6: şase iteraţii sunt necesare

pentru a obţine convergenţă. Iteraţia zero se referă la rezultatele

obţinute prin rularea simulării cu valorile parametrilor obţinute cu

modelul cu parametrii distribuiţi. Modelul vasculaturii distale

reprezintă o bună aproximare a modelului multiscalar, deoarece

toate obiectivele se află în interiorul unui interval de 5% în jurul

valorii de referinţă corespunzătoare. Totuşi, pentru ca valorile

calculate să fie perfect aliniate cu cele de referinţă, valorile

parametrilor trebuie adaptate semnificativ. Acest aspect este cauzat

pe de o parte de faptul că geometria aortei conduce la o creştere

considerabilă a complianţei totale a modelului multiscalar, iar pe de

altă parte de faptul că segmentul coarctaţei duce la o creştere a

rezistenţei totale a modelului multiscalar. În continuare se compară

valorile simulate şi măsurate ale mărimilor variabile în timp din aorta

ascendentă şi descendentă.

Figura 2.7 prezintă o comparaţie a presiunilor de intrare, şi a

presiunilor şi debitelor din aorta descendentă, obţinute prin

măsurători respectiv prin simulări, în cele două configuraţii descrise

mai sus (cu arbore structurat sau cu model windkessel la frontiera de

ieşire). Presiunea din aorta descendentă obţinută prin utilizarea

condiţiei de frontieră de tip arbore structurat este mai apropiată de

cea măsurată decât presiunea obţinută cu modelul windkessel:

atunci când se folosesc arbori structuraţi presiunea simulată de la

începutul sistolei este identică cu cea măsurată. Suplimentar, atunci

când se foloseşte modelul windkessel se pot observa oscilaţii

pronunţate de-a lungul diastolei. Aceste oscilaţii sunt absente sau

semnificativ amortizate în cazul utilizării arborilor structuraţi. În cazul

debitului din aorta descendentă, se obţin rezultate mai bune cu

modelul windessel.

Pentru a compara în mod cantitativ rezultatele, s-au calculat

diferenţele medii şi maxime între mărimile simulate şi calculate.

Acestea sunt prezentate în tabelul 2.1 şi confirmă observaţiile de mai

sus. La presiuni diferenţele sunt mai mici în cazul utilizării arborilor

structuraţi, în timp de modelul windkessel conduce la diferenţe mai

mici pentru debitul aortei descendente.

În final, în tabelul 2.2 se prezintă valorile măsurate şi simulate ale

căderii de presiune de-a lungul coarctaţiei. Rezultatele

hemodinamice confirmă evaluarea vizuală a anatomiei, şi anume

coarctaţia nu este semnificativă, conducând la o cădere de presiune

de doar 6.5 mmHg. Ambele valori simulate sunt apropiate de

valoarea măsurată, dar o eroare mai mică se obţine în cazul utilizării

arborilor structuraţi.

Page 8: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

8

http://heart.unitbv.ro/

Figura 2.6. Evoluţia procedurii de estimare a parametrilor modelului aortei proxime.

Figura 2.7. Comparaţie între presiunile şi debitele variabile în timp măsurate şi simulate, obişnuite pentru configuraţiile cu arbori structuraţi şi cu modele windkessel.

Page 9: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

9

http://heart.unitbv.ro/

Tabelul 2.1. Diferenţe medii şi maxime între mărimile variabile în timp măsurate şi simulate,

obţinute pentru configuraţiile cu arbore structurat şi cu model windkessel.

Diferenţă Condiţie de frontieră de ieşire Presiunea aortei ascendente

[mmHg]

Presiunea aortei descendente

[mmHg]

Debitul aortei descendente

[ml/s]

Medie Arbore structurat 3.35 3.00 20.19

Windkessel 3.73 4.36 11.47

Maximă Arbore structurat 7.60 9.29 46.24

Windkessel 8.14 11.69 37.11

Tabelul 2.2. Comparaţie între căderile de presiune medii şi maxime de-a lungul coarctaţiei, obţinute prin măsurare invazivă, respectiv din simulările hemodinamice.

Configuraţie ΔP Mediu

[mmHg]

ΔP Maxim

[mmHg]

Simulat – Condiţie de frontieră de tip windkessel 1.17 8.53

Computed – Condiţie de frontieră de tip arbore structurat 1.36 7.20

Măsurat 1.23 6.50

2.3. Simularea tridimensională a interacţiunii fluid-solid

2.3.1. Metodologie

S-a dezvoltat o metodologie pentru simulări hemodinamice

tridimensionale cu pereţi complianţi. În prima fază s-a considerat o

interacţiune fluid-solid unidirecţională: mişcarea peretelui arterial

este cunoscută şi impusă în simularea hemodinamică. Pentru

simularea curgerii fluidului s-a folosit metoda Lattice Boltzmann

[Aidun et al, 2010], [Chen et al., 1998], [Yu et al., 2003].

Geometria variabilă în timp este dată sub forma unui set de

suprafeţe poligonale, care la fiecare moment de timp au acelaşi

număr de noduri. Pentru fiecare nod se calculează tranformata

Fourier discretă, care este apoi folosită pentru a determina poziţia şi

viteza peretelui la orice moment de timp. Deoarece mişcarea

peretelui este de obicei periodică şi lină, spectrul Foruier va conţine

un număr mic de armonice. Transformata Fourier se calculează într-o

etapă de pre-procesare, iar în timpul simulării doar transformata

inversă este evaluată pentru a determina geometria la un anumit

moment de timp. Viteza fiecărui nod al peretelui este calculată

analitic ca derivată în timp a transformatei inverse.

Când geometria este actualizată anumite noduri de fluid devin noduri

de solid şi invers. Pentru a evita reconfigurarea gridului, toate

nodurile care sunt etichetate ca fiind noduri fluid la cel puţin un

moment de timp discret sunt identificate pe parcursul pre-procesării

(acest aspect este important pentru alocarea corectă a memoriei).

Funcţia level-set, care este o reprezentare implicită a suprafeţei

trebuie şi ea actualizată la fiecare moment de timp. Deoarece

valoarea exactă a funcţiei este importantă doar în apropierea

peretelui, doar câteva straturi de noduri din jurul peretelui sunt

considerate. Semnul funcţiei este determinat pe baza direcţiei

vectorului normal la suprafaţa peretelui.

Deoarece actualizarea geometriei este o operaţie care necesită un

timp de execuţie semnificativ, aceasta nu este realizată la fiecare

moment de timp. Intervalul de timp după care geometria este

actualizată este ales astfel încât pentru orice nod deplasarea maximă

să fie 0.5Δx. În acest sens, la pre-procesare se determină viteza

maximă (în timp şi spaţiu), iar intervalul de timp după care trebuie

actualizată geometria este calculat cu 0.5 Δx / umax. Această abordare

conduce la o scădere semnificativă a timpului total de execuţie

deoarece pasul de timp al simulării LBM (10-5 secunde) este mult mai

mic decât intervalul de timp la care trebuie actualizată geometria.

Deoarece viteza peretelui este descrisă prin spectrul Fourier,

determinarea analitică a vitezei maxime nu este simplă. De aceea, se

foloseşte o metodă numerică, care determină valoarea maximă

printr-o abordare iterativ. Când geometria este actualizată, viteza

peretelui up este asociată fiecărui nod apropiat de perete împreună

cu noua valoare a funcţiei level set. Vitezele peretelui sunt impuse la

pasul de streaming al simulării LBM.

2.3.2. Rezultate

Pentru a valida metodologia descrisă mai sus s-a considerat iniţial un

experiment sintetic pentru care rezultatul analitic este cunoscut:

curgerea unui fluid printr-un cilindru care se expandează şi contractă

periodic. Geometria a fost generată sintetic prin aplicarea unei funcţii

de deformare a unui cilindru drept.

La frontierele stângă şi dreapta ale domeniului au fost impuse

condiţii de frontieră de tip presiune constantă cu o constrângere de

tip valvă: debitul este setat la valoarea zero când diferenţa de

presiune îşi schimbă semnul. Astfel, fluidul poate să intre în domeniu

doar prin frontiera stângă şi să-l părăsească doar prin frontiera

dreaptă. Figura 2.8 prezintă configuraţia geometriei împreună cu

vectorii de viteză pentru fazele de expandare şi contracţie.

Figura 2.8. Vectorii de viteză de-a lungul unui secţiuni longitudinale a unui

cilindru cu volum variabil.

S-au efectuat simulări în configuraţia descrisă mai sus şi debitele au

fost comparate cu soluţia analitică. Deoarece fluidul este

incompresibil, debitul Q trebuie să fie egal cu rata de modificare a

volumului dV/dt. Debitele au fost calculate pe planele localizate la

frontierele stânga şi dreapta şi rezultatele sunt prezentate în figura

2.9. Se poate observa că rezultatele obţinute sunt aproape identice

Page 10: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

10

http://heart.unitbv.ro/

cu soluţia analitică. Diferenţa medie între rezultatul simulării şi

soluţia a fost de 0.11 mm3/s iar diferenţa maximă de 0.34 mm

3/s,

ceea ce reprezintă o eroare normalizată de 1.45% şi respectiv 4.25%.

Figura 2.9. Cilindru cu volum variabil: comparaţie între debitele măsurate şi

soluţia analitică.

În continuare modelul a fost aplicat pentru studiul curgerii

peristaltice tridimensionale: o deformare periodică a pereţilor

generează o mişcare a fluidului. În acest caz, geometria este dată de

un vas cilindric, asupra căruia se aplică o funcţie de deformare

periodică. În figura 2.10 se prezintă geometria împreună cu vectorii

de viteză.

Figura 2.10. Contururile componentei axiale a vitezei: debit pozitiv în regiunea

aflată în faza expandare şi debit negativ în regiunea aflată în faza de

contracţie.

Frontierele de intrare şi de ieşire sunt periodice. Mişcarea fluidului

este generată doar de mişcarea peretelui. După cum a fost descris în

literatura de specialitate [Shapiro et al., 1969], debitul volumetric

pentru număr Reynolds mic poate fi determinat analitic. Debitul

mediu din simulare este evaluat numeric prin integrarea câmpului de

viteze pe o secţiune transversală şi medierea rezultatului de-a lungul

unui ciclu.

Pentru a realiza o validare detaliată a metodologiei s-au efectuat

experimente cu diferite valori ale parametrului Φ (0.16, 0.27, 0.38,

0.49, 0.6) şi debitele obţinute prin simulare s-au comparat cu cele

calculate analitici. Ceilalţi parametrii au fost setaţi la valorile a

0.5mm, c 8 mm/s şi λ 2mm. Vâscozitatea cinematică a fost ν 1

mm2/s, care conduce la un număr Reynolds egal cu 1.

Figura 2.11 prezintă rezultatele obţinute: debitele rezultate în

simulare sunt foarte apropiate de cele analitice, eroarea maximă

absolută a fost de 0.17 mm3/s, iar eroare medie abosultă de 0.11

mm3/s.

Figura 2.11. Simularea curgerii peristaltice. Comparaţie între debitul mediu

măsurat şi soluţia analitică.

De asemenea, au fost efectuate experimente cu diferite valori ale

numărului Reynolds. Numărul Reynolds a fost modificat prin

adaptarea vitezei de undă c, iar ceilalţi parametrii au fost menţinuţi la

valori constante. În figura 2.12 se prezintă liniile de undă ale debitului

pentru curgeri cu număr Reynolds variabil între 1 şi 100. Se poate

observa că cele două vortexuri, care se rotesc în direcţii opuse, devin

mai mici atunci când numărul Reynolds creşte. De asemenea,

eficienţa procesului de mişcare a fluidului scade: în comparaţie cu

debitul analtici, debitul măsurat este cu 13% mai mic pentru Re = 10

şi cu 35% mai mic pentru Re = 100.

Figura 2.12. Liniile de undă ale curgerii peristaltice, obţinute pentru diferite valori ale numărului Reynolds.

3.

Achiziția datelor de la pacienți

În vederea validării modelelor multiscalare dezvoltate s-au

achiziționat date de la un număr mare de pacienți cu stenoze

coronariane (pacienți stabili care au indicație de angiografie în

vederea diagnosticării circulației coronariene) la Spitalul Clinic de

Urgență București. Datele achiziționate de la pacienți pot fi

Page 11: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

11

http://heart.unitbv.ro/

împărțite în trei categorii, în funcție de modalitatea de achiziție:

angiografie, ecografie, cateter cu senzor de presiune.

Până în prezent au fost incluşi în studiu 48 de pacienţi şi s-au

efectuat un total de 120 de măsurători de presiune şi FFR, dintre

care 92 înainte de stentare şi 28 după stentare. Cele 120 de

măsurători de presiune sunt distribuite astfel pe cele 3 ramuri

principale:

LAD: 62

LCx: 33

RCA: 25

Caracteristicile pacienţilor de la care s-au achiziţionat date sunt

prezentate în tabelul 3.1.

În vederea stocării anonimizate a datelor s-a dezvoltat un pachet

software, folosit de partenerul clinic, care anonimizează datele

angiografice, ecografice şi de presiune într-un mod automatizat.

Tabelul 3.1: Caracteristicile pacienţilor incluşi in studiu.

Caracteristică Valoare

Sex 76.9%

Vârstă 63.9 ± 8.0 ani

BMI 28.8 ± 3.6 kg/m2

Diabet 30.77%

Hipertensiune 84.6%

Hipercolesterolemie 46.1%

Fumator 57.7%

Istoric CAD în familie 30.7%

Infarct anterior 7.7%

PCI anterior 54%

CABG anterior 3.8%

Fracţie de ejecţie 51.71 ± 7.9%

4.

Procesarea datelor achiziționate și validarea modelului multiscalar al circulației

coronariene

4.1. Postprocesarea presiunilor invazive achiziţionate

Presiunile invazive aortice şi coronariene, la repaus şi hiperemie,

achiziţionate de la fiecare pacient (un exemplu este prezentat în

figura 4.1) au fost post-procesate pentru a extrage diferite mărimi de

interes pentru diagnosticare funcţională a stenozelor coronariene:

IFR (instantaneous wave free ratio) de repaus,

IFR hiperemic,

Pd/Pa la repaus,

FFR (Pd/Pa hiperemic).

Figura 4.1. Presiunea aortic (roşu) şi coronariană (verde) a unui pacient inclus în studiu.

Figura 4.2. Identificarea intervalului diastolic fără unde de debit şi presiune; raportul instantaneu al presiunilor distale şi aortice.

Page 12: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

12

http://heart.unitbv.ro/

Mărimile Pd/Pa la repaus şi FFR sunt calculate ca raport al

presiunilor medii aortice şi distale de-a lungul unui ciclu cardiac.

Pentru a putea determina mărimile IFR trebuie iniţial identificat

intervalul diastolic în care undele de presiune şi de debit sunt

absente. În acest sens s-a folosit o abordare descrisă anterior în

literatura de specialitate [Sen et al., 2012], şi anume: începutul

perioadei este identificat la 25% din durata diastolei iar finalul

perioadei este identificat la 5ms înainte de finalul diastolei (figura

4.2). Tot în figura 4.2 este prezentat raportul instantaneu dintre

presiunea coronariană distală şi presiunea aortică. Se poate observa

că acest raport este relativ constant de-a lungul perioadei diastolice

fără unde. IFR reprezintă valoarea medie a raportului de-a lungul

acestui interval.

4.2. Reconstrucţia tridimensională a arterelor coronariane

În continuare s-a folosit platforma software descrisă în secţiunea 1

pentru a genera modelele anatomice tridimensionale ale arterelor

coronariene. În figura 4.3 se prezintă ca exemplu un astfel de model

anatomic. În acest sens se parcurg următorii paşi:

se identifică două serii angiografice care prezintă artera

coronariană cu stenoza de interes;

pe fiecare serie se identifică un cadru de la finalul diastolei în

care calitatea imaginii este bună (se folosesc cadre de la finalul

diastolei deoarece în acea etapă a ciclului cardiac mişcarea

arterelor este minimă iar calitatea imaginii este maximă);

se rulează aplicaţia descrisă în secţiunea 1 şi se generează

modelul anatomic tridimensional.

Figura 4.3. Model anatomic tridimensional al arterei LAD.

4.3. Estimarea indicatorilor diagnostici din modelul hemodinamic

În continuare s-au efectuat simulările hemodinamice coronariene

folosind modelele hemodinamice descrise în acest raport şi în

rapoartele anterioare şi modelele anatomice tridimensionale

generate conform descrierii de mai sus. Din rezultatele simulării au

fost estimate următoarele mărimi:

IFR de repaus,

FFR.

Rezultatele obţinute pentru toţi pacienţii procesaţi până în acest

moment sunt prezentate în anexa 1. În tabelul 4.1 sunt prezentate

statisticile de diagnosticare ale mărimilor estimate non-invaziv

(acestea au prefixul ’C-‘). Pentru FFR s-a folosit condiţia ≤ 0.8 pentru

identificarea leziunilor pozitive, iar pentru IFR s-a folosit condiţia ≤

0.89 pentru identificarea leziunilor pozitive. Se poate observa că

rezultatele sunt mult mai bune pentru indicatorul FFR. Acest aspect

este dat de faptul că starea de repaus a pacienţilor are o

variabilitate mai mare decât starea de hiperemie (care este o stare

cu debit coronarian maxim). Această variabilitate nu poate fi

surprinsă în totalitate de modelul hemodinamic.

De asemenea se poate observa că indicatorul FFR a fost estimat

pentru 89 de leziuni în timp ce indicatorul IFR a fost estimat pentru

86 de leziuni. Această diferenţă este dată de faptul că pentru 3

leziuni presiunile de repaus nu au fost disponibile şi prin urmare

valoarea IFR invazivă nu a putut fi extrasă.

În continuare a fost evaluată şi o strategie hibridă bazată pe C-IFR şi

FFR invaziv. Aceasta are la bază un studiu anterior care a evaluat o

strategie hibridă IFR-FFR pentru evaluarea funcţională a stenozelor

coronariene [Petraco et al., 2012]. Conform acestei strategii dacă se

obţine IFR > 0.93 leziunea nu este revascularizată, dacă IFR < 0.86

leziunea este revascularizată, iar leziunile cu valori IFR intermediare

sunt evaluate pe baza indicatorului FFR (folosind un pragul de 0.8).

Strategia hibridă evaluată în cadrul acestei activităţi a fost similară,

dar în loc de IFR s-a folosit C-IFR. În tabelul 4.2 sunt prezentate

statisticile de diagnosticare ale acestei strategii. Dintre cele 86 de

leziuni, 23 au fost semnificative hemodinamic (FFR ≤ 0.8) iar

acurateţea a fost de 96.5%, mult superioară celor de mai sus.

Comparativ cu o strategie 100% invazivă, în care indicatorul FFR este

măsurat invaziv pentru toate leziunile, cu strategia hibridă doar 27

de leziuni necesită măsurători invazive. Dintre cele 60 de leziuni care

sunt clasificate exclusiv pe baza indicatorului C-IFR estimat din

modelul hemodinamic, 17 au avut IFR < 0.86, iar 42 au avut IFR >

0.93. Timpul mediu de calcul pentru estimarea indicatorului C-IFR

pentru leziune a fost de 38.1 ± 8.8 secunde pe un desktop cu

procesor Intel i7 cu 8 core-uri @ 3.4 GHz. De asemenea, prin

aplicarea strategiei hibride 46.5% dintre pacienţi nu ar necesita nici

o măsurătoare invazivă.

Tabelul 4.1: Statisticile de diagnosticare ale indicatorilor IFR şi FFR estimaţi

din simulările hemodinamice.

Indicator statistic Mărime estimată

C-IFR C-FFR

Real pozitiv 15 17

Fals pozitiv 5 7

Real negativ 49 60

Fals negativ 17 5

Sensibilitate 46.87% 77.27%

Specificitate 90.74% 89.55%

Valoarea predictiv pozitivă 75% 70.83%

Valoarea predictiv negativă 74.24% 92.30%

Acurateţe 74.41% 86.51%

Tabelul 4.2: Statisticile de diagnosticare ale strategiei hibride C-IFR - FFR.

Indicator statistic Valoare

Real pozitiv 21

Fals pozitiv 2

Real negativ 62

Fals negativ 1

Sensibilitate 95.45%

Specificitate 96.87%

Valoarea predictiv pozitivă 91.30%

Valoarea predictiv negativă 98.41%

Acurateţe 96.51%

Page 13: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

13

http://heart.unitbv.ro/

5.

Obiective formulate pe baza presiunilor

S-a dezvoltat o nouă procedură de personalizare pentru modelul cu

parametrii distribuiţi al întregului circuit cardiovascular. În acest sens

numărul de obiective şi de parametrii personalizaţi a fost crescut la

20:

6 obiective de bază pentru circulaţia sistemică

8 obiective speciale pentru circulaţia sistemică

6 obiective de bază pentru circulaţia pulmonară

Toate obiectivele sunt bazate pe presiuni şi volume ventriculare şi

aortice.

5.1. Obiective formulate pe baza presiunilor

Obiectivele de bază extrase din măsurătorile de presiune sunt:

Presiunea ventriculară maximă

Presiunea aortică medie

Presiunea aortică minimă

Intervalul de timp de-a lungul căruia valva aortică este deschisă.

Acestea sunt formulate atât pentru circulaţia sistemică cât şi pentru

circulaţia pulmonară.

Obiectivele speciale formulate pe baza presiunilor sunt prezentate în

figura 5.1. Acestea sunt:

Panta 1: panta presiunii ventriculare la începutul sistolei (înainte

de deschiderea valvei aortice)

Presiunea ventriculară diastolică

Panta 2: panta presiunii aortice de-a lungul diastolei.

5.2. Obiective formulate pe baza volumelor

Obiectivele de bază extrase din măsurătorile de volum sunt:

Volumul ventricular maxim (la finalul diastolei)

Volumul ventricular minim (la finalul sistolei)

Acestea sunt formulate atât pentru circulaţia sistemică cât şi pentru

circulaţia pulmonară. Obiectivele speciale formulate pe baza

volumelor sunt prezentate în figura 5.2.

În vederea formulării acestor obiective, pacienţii au fost împărţiţi în

trei categorii, care sunt enumerate în continuare împreună cu

obiectivele speciale:

Pacienţi cu platou diastolic:

o Nivelul platoului diastolic, exprimat în procentaje în raport

cu valorile minime şi maxime ale volumului

o Δt1: intervalul de timp dintre momentele la care volumul

este egal cu media dintre volumul minim şi volumul la

platoul diastolic

o Panta 1: panta curbei de volum la cel de-al doilea moment

de timp utilizat pentru determinarea lui Δt1

o Δt2: durata platoului diastolic

o Panta 2: panta curbei de volum când atriul se contractă

(panta maximă a curbei de volum după platoul diastolic)

Pacienţi fără creştere a volumului în ultima parte a diastolei:

o Δt1: intervalul de timp dintre momentele la care volumul se

află la 30% între valoarea minimă şi maximă

o Panta 1: panta curbei de volum la cel de-al doilea moment

de timp utilizat pentru determinarea lui Δt1

Pacienţi fără platou diastolic, dar cu diferite pante ale volumului:

o Valoarea volumului la care începe contracţia volumului

(punctul de tranziţie de la prima pantă la cea de-a doua

pantă a volumului diastolic)

o Panta 1: panta curbei de volum în prima parte a diastolei,

înainte de contracţia atriuluiâ

o Panta 2: panta curbei de volum în a doua parte a diastolei,

după contracţia atriului

o Δt1: intervalul de timp dintre momentele la care volumul

este egal cu media dintre volumul minim şi volumul la care

se declanşează contracţia atriului.

Figura 5.1. Obiectivele speciale formulate pe baza presiunilor.

Page 14: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

14

http://heart.unitbv.ro/

Figura 5.2. Obiectivele speciale formulate pe baza volumelor.

Deoarece numărul de obiective este foarte mare, la început se

realizează o iniţializare iterativ-individuală a parametrilor. Pentru

acest pas se realizează o asociere între obiectivele definite şi

parametrii care trebuie personalizaţi (tabelul 5.1 prezintă asocierea

pentru circulaţia sistemică). Această asociere a fost stabilită în urma

realizării unei analize de sensibilitate: pentru fiecare obiectiv s-a

determinat parametrul care are cea mai mare influenţă asupra

obiectivului respectiv. Apoi fiecare parametru este adaptat pe rând

pentru a îndeplini obiectivul care îi corespunde. După ce au fost

adaptaţi toţi parametrii, procedura este reluată până când valorile

obiectivelor extrase din simulare se apropie de valorile de referinţă.

Reluarea procedurii este necesară deoarece prin personalizarea

succesivă a parametrilor, obiectivele ale căror valori au fost apropiate

de cele de referinţă se vor îndepărta din nou de aceste valori prin

modificarea altor parametrii care au influenţă asupra acestor

obiective.

Tabelul 5.1. Corespondenţa între obiectivele şi parametrii circulaţiei sistemice în etapa de iniţializare.

Obiectiv Parametrii

Presiunea aortică medie Volumul iniţial al sistemului în circuit închis

Valoare medie a volumului ventriculului Volum mort al ventriculului

Interval de timp valva deschisa Momentul de timp la care se atinge elastanţa maximă

Δt1 Momentul de final al curbei de elastanţă

Presiune diastolică Valoarea minimă a elastanţei ventriculului

Nivelul platoului diastolic Elastanţa minimă a atriului

Panta 1 volum Rezistenţa valvei mitrale

Δt2 – durata platoului Moment de timp la care se declanşează contracţia atriului

Panta 2 volum Elastanţa maximă a atriului

Panta 1 presiune Elastanţa maximă a ventriculului stâng

Presiune de puls aortică Raportul dintre rezistenţa proximă şi rezistenţa distală sistemică

Panta 2 presiune Complianţa sistemică

5.3. Rezultate

Metodologia descrisă mai a fost aplicată cu succes pentru

personalizarea a 52 de pacienţi. Dintre aceştia 36 de pacienţi au avut

platou diastolic, 11 nu au avut creştere a volumului în ultima parte a

diastolei, iar 5 pacienţi au fost fără platou diastolic dar cu diferite

pante ale volumului.

În figura 5.3 se prezintă rezultatele obţinute pentru unul din pacienţii

cu platou diastolic, care reprezintă totodată şi categoria cu numărul

cel mai mare de pacienţi. Se poate observa că rezultatele simulării

sunt foarte apropiate de valorile măsurate invaziv.

În continuare cele 8 obiective speciale formulate pentru circulaţia

sistemică vor fi determinate şi pentru circulaţia pulmonară, astfel

încât numărul de obiective şi parametrii va creşte la 28.

Page 15: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

15

http://heart.unitbv.ro/

Figura 5.3. Comparaţie între presiunile şi volumele măsurate respectiv obţinute din simulare pentru unul din pacienţii cu platou diastolic.

6.

Utilizarea dimensiunii fractale în identificarea stenozelor coronariene

6.1. Introducere

Dimensiunea fractală este o măsură a complexității morfologice a

structurilor biologice. Se urmărește dacă pentru diferite imagini

medicale afectate de diverse boli vasculare au loc schimbări în

magnitudinea dimensiunii fractale. Comparația a fost făcută cu

regiunile de interes care conțin structuri morfologice normale.

Analiza fractală reprezintă o metodă obiectivă de cuantificare a

complexității geometrice și a rugozității conturului unei imagini cu

stenoză.

6.2. Procesarea spațială

Un pas important în îmbunătățirea diagnosticului unor imagini

medicale este procesul de segmentare care poate fi privit ca o

clasificare a pixelilor având ca scop localizarea sau detectarea

conturului obiectelor. Acest obiectiv nu este atins în mod

corespunzător în cazul în care zgomotul nu este filtrat corect în

imaginea originală. Prin urmare, o îmbunătățire a imaginii este o

operațiune crucială pentru angiografii.

Frangi et al. au propus o filtrare neliniară specifică extragerii vaselor

de sânge din imaginii [Frangi et al., 1999]. Această metodă utilizează

matricea Hessiană pentru calcularea vectorilor proprii care în funcție

de caracteristicile lor estimează pixelii susceptibili de a reprezenta

vasele. În continuare, a fost folosită metoda de prag (thresholding)

pentru a genera imagini alb-negru iar pragul global al intensității

pixelilor este stabilit manual sau prin utilizarea metodei Otsu [Otsu,

1979]. În ceea ce privește detectarea conturului, algoritmul lui

Canny calculează gradientul imaginii în pixeli pentru o imagine care

a fost prefiltrată cu filtrul Gaussian. Pentru corecție, operațiile

morfologice sunt aplicate înainte și după această etapa pentru a

elimina pixelii nedoriți.

Există mai multe metode de a calcula dimensiunea fractală, iar cele

mai comune sunt dimensiunea Hausdorff, Minkowski și

dimensiunea box-counting. Cea din urmă presupune umplerea

imaginii cu figuri geometrice, cum ar fi pătrate de dimensiuni

diferite. Pornind de la scala cea mai mare, care este chiar

dimensiunea imaginii, algoritmul este reiterat cu un anumit pas și se

oprește atunci când dimensiunea cutiei devine de 2x2 pixeli,

numărând permanent pătratele care conțin conturul obiectului. În

cele din urma, se generează graficul evoluției numărului de pătrate

în raport cu scala în vederea obținerii dimensiunii fractale care este

panta dreptei de regresie a graficului in scară logaritmică.

6.3. Rezultate

Pentru angiografiile cardiace, două situații diferite au fost

considerate, unul pentru cazul sănătos și altul pentru vasele cu

stenoză. Un număr de 6 pacienți și 8 imagini au fost testate și

rezultatele sunt prezentate în tabelul 6.1.

Tabelul 6.1. Dimensiune fractală în angiografia cardiacă.

Page 16: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

16

http://heart.unitbv.ro/

Comparând media dimensiunilor fractale pentru vasul sănătos cu cel

stenozat se observă o creștere la nivelul sutimilor. Mai mult decât

atât, s-a încercat determinarea unei corelații între dimensiunea

fractală și gradul stenozei calculat din Quantitative Coronary

Angiography (QCA) care este un grafic des utilizat de medici pentru

a analiza evoluția ariei axiale a vasului de-a lungul lungimii

segmentului de vas. Indicele PAS-procetanjul ariei stenozei, calculat

ca diferență între aria maximă și minimă normată apoi la valoarea

minimă este calculată după reconstrucția 3D a conturului vaselor și

comparată cu dimensiunea fractală pentru cele două proiecţii și este

ilustrată în tabelul 6.2.

În concluzie, metoda implementată demonstrează că se pot extrage

caracteristici interesante și măsurători din imagini cu patologii

vasculare. Metoda box counting a fost aleasă în studiul de față

pentru că este o tehnică simplă și robustă pentru a caracteriza în

mod automat stenozele din angiografiile cu raze X și care este

corelată cu rezultatele analizei Quantitative Coronary Angiography.

Totuși, cea mai mare dificultate rămâne procesarea imaginii care

trebuie adaptată tipului de achiziție și necesită o supervizare din

partea unui expert imagistic. De asemenea, rezoluția spațială bună

este o premiză a succesului acestei metode.

Tabelul 6.2. Corelarea dimensiunii fractale cu 3D QCA.

7.

Utilizarea procesării paralele în accelerarea timpilor de executie al algoritmilor

dezvoltati

Activităţile de procesare paralelă din această perioadă de raportare

s-au concentrat pe două activităţi, care sunt prezentate mai jos:

Implementare multi-CPU a modelului hemodinamic utilizat

pentru simulările coronariene

Implementare pe GPU a unui solver multigrid

7.1. Implementare multi-CPU a modelului hemodinamic utilizat

pentru simulările coronariene

În vederea reducerii timpului de execuţie al modelului hemodinamic

utilizat pentru rezultatele prezentate în secţiunea 4 s-a realizat o

implementare multicore a acestuia. S-au parcurs următoarele etape:

identificarea componentelor aplicaţiei care ocupă cea mai mare

parte a timpului de execuţie

identificarea componentelor aplicaţiei care sunt paralelizabile

paralelizarea acestor componente

evaluarea îmbunătăţirilor obţinute din punctul de vedere al

timpului de execuţie

În tabelul 7.1 sunt prezentate principalele operaţii ale aplicaţiei:

iniţializare model arterial: în această etapă se procesează

modelul anatomic, se identifică stenozele, se estimează debitul

de intrare şi valorile initţiale ale parametrilor modelului

microvascular

alocare şi iniţializare memorie: se alocă memorie pentru

tablourile de date care urmează să stocheze mărimile

necunoscute: debit, presiune, arii transversale, etc.

procesare nod frontieră de intrare: determinarea valorilor

mărimilor necunoscute la frontiera de intrare pentru fiecare

moment de timp discret; deoarece debitul este cunoscut doar

presiunea şi aria transversală trebuie actualizate

procesare noduri interioare: determinarea valorilor mărimilor

necunoscute la nodurile interioare ale segmentelor arteriale

pentru fiecare moment de timp discret

procesare noduri bifurcaţie: determinarea valorilor mărimilor

necunoscute la nodurile bifurcaţiilor arteriale pentru fiecare

moment de timp discret

procesare noduri frontieră de ieşire: determinarea valorilor

mărimilor necunoscute la nodurile frontierelor de ieşire pentru

fiecare moment de timp discret

postprocesare: extragerea mărimilor de interes din rezultatele

simulării hemodinamice (FFR, IFR)

alte operaţii: cuprinde toate celelalte operaţii care ocupă un timp

de execuţie nesemnificativ, precum cele de actualizare a valorilor

parametrilor modelului microvascular în vederea personalizării

modelului

Dintre toate aceste operaţii doar cele de procesare a nodurilor

interioare, de bifurcaţie şi de la frontierele de ieşire pot fi paralelizate

eficient. Datorită schemei de discretizare utilizate pentru rezolvarea

numerică a modelului hemodinamic (Lax-Wendroff), toate calculele

efectuate pentru nodurile dintr-o anumită categorie pot fi realizate în

paralel (sunt independente de la un nod la altul). Pe de altă parte,

ordinea corectă de actualizare a nodurilor din diferite categorii este:

nod de intrare, noduri interioare, noduri de bifurcaţie, noduri de la

frontierele de ieşire.

Page 17: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

17

http://heart.unitbv.ro/

În următoarea etapă cele trei operaţii identificate mai sus au fost

paralelizate folosind openMP [***openMP, 2015]. În tabelul 7.1 se

prezintă timpii de execuţie obţinuţi pentru implementarea

secvenţială, respectiv pentru implementarea multicore. Valorile din

tabel reprezintă valorile medii obţinute pe baza a 86 de simulări

hemodinamice. Timpii medii de execuţie au fost de 230.9 ± 59.7

secunde pentru aplicaţia secvenţială şi 38.1 ± 8.8 secunde pentru

aplicaţia multicore. În tabel sunt prezentate doar valorile medii fără

deviaţii standard. Testele au fost rulate pe un procesor Intel i7 cu 8

core-uri @ 3.4 GHz: pentru aplicaţia secvenţială s-a folosit un singur

core al procesorului iar pentru aplicaţia multicore s-au folosit toate

cele opt core-uri.

Pornind de la timpul de execuţie secvential şi de la legea lui Amdahl,

timpul de execuţie minim în cazul unei paralelizări optime ale celor

trei operaţii menţionate mai sus ar fi de 32.69 secunde. Timpul

mediu real obţinut cu implementarea paralelă este de 38.1 secunde,

gradul de accelerare al operaţiilor paralelizate variind între 5.11 şi

7.33. În medie numărul de noduri interioare ale unui model arterial

variază între 300 şi 500, numărul de noduri de bifurcaţie între 40 şi

50 iar numărul de noduri de frontieră de ieşire între 12 şi 20. Prin

urmare, cu cât numărul de operaţii paralele a fost mic cu atât şi

gradul de accelerare a fost mai mic. Gradul de accelerare global, al

întregii aplicaţii, a fost 6.06.

Tabelul 7.1. Analiza timpilor de execuţie pentru implementarea multicore.

Operaţie Secvenţial (1

core)

Multicore

(8 core-uri)

Accelerare

Iniţializare model arterial 1.43 1.43 -

Alocare şi iniţializare memorie 0.78 0.78 -

Procesare nod frontieră de intrare 0.54 0.54 -

Procesare noduri interioare 127.3 17.36 7.33

Procesare noduri de bifurcaţie 72.54 11.16 6.50

Procesare noduri frontieră de ieşire 26.8 5.24 5.11

Postprocesare 0.72 0.72 -

Dealocare memorie 0.33 0.33 -

Alte operaţii 0.56 0.56 -

Total 230.9 38.1 6.06

7.2. Implementare pe GPU a unui solver multigrid

7.2.1. Metodologie

Cei mai populari algoritmi pentru rezolvarea sistemelor de ecuaţii

liniare sparse sunt Preconditioned Conjugate Gradient (PCG) şi

Multigrid (MG) în două variante: multigrid geometric (GMG) şi

multigrid algebraic (AMG). Ambele metode se bazează pe algoritmi

iterativi, care sunt mai rapizi decât algoritmi de rezolvare directă. În

cadrul acestei activităţi ne-am concentrat pe algoritmul GMG. Acesta

se bazează pe o ierarhie de discretizări şi pe principiul conform căruia

corecţiile efectuate la nivelele de discretizazare coarse îmbunătăţesc

rata de convergenţă a soluţiei de la nivelul de discretizare fin.

Studii publicate anterior au demonstrat deja că implementările GMG

bazate pe GPU sunt mai rapide decât cele bazate pe CPU [Williams et

al., 2012], [Geveler et al., 2013], [Anzt et al., 2012]. De aceea în

cadrul acestei activităţi ne-am concentrat asupra unei analize mai

detaliate a algoritmului GMG implementat pe GPU. Pentru a analiza

performanţa, s-a considerat ecuaţia Poisson. Prin discretizare se

obţine un sistem de ecuaţii liniare. În mod particular s-a considerat

ecuaţia staţionară de conducţie a căldurii şi mai multe scheme de

discretizare cu diferenţe finite centrale: 7 puncte, 19 puncte şi 27 de

puncte [Stroia et al., 2015]. Aceste scheme de discretizare sunt

folosite pentru a genera mai multe grid-uri. Toate variantele GMG

sunt bazate pe tranziţii succesive de la griduri fine la griduri coarse şi

inapoi. Prin urmare operaţiile de bază ale algoritmului GMG sunt:

relaxarea: o metodă iterativă de bază precum Jacobi sau Gauss-

Seidel este folosită pentru a diminua erorile de frecvenţă mare

ale soluţiei numerice

restricţie: reziduul de pe un grid fin este folosit pentru a calcula

reziduul pe un grid coarse

interpolare: reziduul de pe un grid fin este determinat prin

interpolare pe baza valorilor de pe gridul coarse

Figura 7.1. Conceptul de bază al metodei multigrid.

În cadrul acestui studiu s-au folosit pentru pasul de relaxare

metodele red black Gauss-Seidel (RBGS) pentru discretizarea cu 7

puncte şi Jacobi pentru discretizarea cu 19 şi 27 de puncte. Metoda

RBGS necesită un singur tablou de valori, dar calculele sunt efectuate

în două etape secvenţiale: iniţial se actualizează nodurile marcate ca

fiind roşii şi apoi nodurile marcate ca fiind negre. Metoda Jacobi

necesită o singură etapă de calcule.

Variantele GMG utilizate în acest studiu sunt prezentate în figura 7.2:

ciclu V, ciclu W, full MG (FMG). În timp ce la schemele V şi W iteraţiile

Page 18: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

18

http://heart.unitbv.ro/

pornesc de la gridul cel mai fin, la schema FMG iteraţiile pornesc de

la nivelul coarse. Pentru comparaţia performanţei s-a folosit o

implementare optimizată bazată pe GPU a algoritmului PCG [Nita et

al., 2014].

Figura 7.2. Variante GMG: (a) ciclu V, (b) ciclu W, (c) Full MG (FMG).

Schemele V şi W folosesc algoritmul cu ciclu µ (algoritmul 7.1), care

este un algoritm recursiv. Prin parametrul µ se selectează tipul

algoritmului: µ = 1 ciclu V, µ = 2 ciclu W. Valorile n1, n2, n3

determină numărul iteraţiilor de relaxare de pe ramura descendentă,

la nivelul coarse şi repsectiv de pe ramura ascendentă. Suplimentar

faţă de operaţia de restricţie, pe ramura ascendentă se execută şi o

operaţie de corectare.

Algoritmul 7.1: Ciclu µ

µ-Cycle(nivel) daca ( nivel este nivelCoarse ) aplică n2 paşi de relaxare altfel aplică n1 paşi de relaxare calcul reziduu restricţie către grid mai coarse µ-Cycle(nivel+1) interpolare corecţie aplică n3 paşi de relaxare final

Metoda multigrid necesită un tablou de stocare a datelor pentru

fiecare nivel de discretizare. Nivelul 0 utilizează cel mai mare tablou

şi ocupă cea mai mare parte a timpului de execuţie total.

7.2.2. Rezultate

Pentru evaluarea performanţelor s-a folosit un GPU NVIDIA GTX Titan

Black şi CUDA toolkit 6.0. Fiecare configuraţie este descrisă pintr-un

număr cu trei cifre care reprezintă valorile celor trei parametrii.

Algoritmul GMG este executat până când valoarea medie a reziduului

nu mai scade de la o iteraţie la alta. Iniţial au fost comparate diferite

scheme GMG cu o configuraţie 313 şi relaxare RBGS. Figura 7.3.

prezintă evoluţia reziduului mediu în raport cu timpul de execuţie.

Schema V conduce la cea mai bună performanţă: deşi necesită mai

multe iteraţii decât celelalte scheme, reziduul mediu scade cel mao

rapid până la 1e-14. Prin urmare în continuare s-a folosit doar

schema V.

Figura 7.3. Comparaţie între performanţele schemelor V, W şi FMG.

În continuare s-a analizat efectul configuraţiei de relaxare. Figura 7.4

prezintă primele patru configuraţii din punctul de vedere al

performanţelor: 2 sau 3 iteraţii de relaxare sunt necesare pe

parcursul operaţiilor de restricţie şi interpolare, şi o singură iteraţie

este necesară pe gridul cel mai coarse. Configuraţiile 212 şi 312

conduc la cea mai bună performanţă, cu uşoare avantaje pentru 212.

Figura 7.4. Efectul configuraţiei de relaxare asupra performanţei algoritmului GMG.

În continuare se analizează timpul de execuţie ocupat de fiecare

operaţie a algoritmului GMG (tabelul 7.2). Pasul de relaxare ocupă cel

mai mult timp (aproape jumătate din timpul total). De asemenea,

tranziţia de la nivelul coarse la nivelul fin necesită mai mult timp

decât tranziţia de la nivelul fin la nivelul coarse. Acest aspect este dat

de faptul că la interpolare nivelul destinaţie al operaţiei este mai fin

decât la restricţie.

În final s-a comparat varianta GMG cu cea mai bună performanţă

(ciclu V, configuraţie de relaxare 312) cu metoda PCG optimizată. S-

au considerat diferite rezoluţii ale gridului fin precum şi cele trei

Page 19: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

19

http://heart.unitbv.ro/

scheme de discretizare. Rezultatele sunt prezentate în tabelul 7.3 şi

indică faptul că pentru un grid fin de 129x129x129 GMG conduce la

un timp de execuţie de 7.1 – 9.2 ori mai mic decât algoritmul PCG şi

în acelaşi timp conduce şi la un reziduu mediu mai mic. Gradul de

accelerare este mai mic pe gridurile mai mici, ceea ce poate fi

justificat prin paralelismul scazut al algoritmului GMG pe aceste

griduri. În practică însă se folosesc în mod tipic griduri cu mai mult de

un milion de noduri.

Trebuie menţionat faptul că într-adevăr GMG oferă aceste avantaje

din punctul de vedere al timpului de execuţie în comparaţie cu PCG,

dar necesită informaţii suplimentare despre ecuaţiile analitice pentru

a putea genera ecuaţiile numerice la diferite nivele de discretizare.

Tabel 7.2. Timpul relativ de execuţie ocupat de fiecare operaţie.

Operaţie % din timpul de execuţie

Relaxare 46.7%

Interpolare 26.7%

Calcul reziduu 14.4%

Corecţie 6.1%

Restricţie 5.9%

Memset(0) 0.2%

Tabelul 7.3. Comparaţie între timpii de execuţie şi reziduurile medii pentru metodele GMG şi PCG.

Rezoluţia gridului fin 129x129x129 65x65x65 33x33x33

Metoda PCG GMG PCG GMG PCG GMG

RBGS

Discret. cu 7 puncte

Eroare medie 5.22e-12 1.44e-14 1.41e-11 1.66e-14 4.21e-11 1.43e-14

Timp [ms] 1118 121 124 50 28 33

Jacobi

Discret. cu 19 puncte

Eroare medie 5.21e-12 7.00e-14 1.25e-11 6.99e-14 3.89e-11 6.79e-14

Timp [ms] 1255 172 127 48 28 32

Jacobi

Discret. cu 27 puncte

Eroare medie 4.30e-12 1.49e-13 1.40e-11 1.17e-13 2.94e-11 1.19e-13

Timp [ms] 1502 211 145 61 29 35

Referințe

[Aidun et al, 2010] Cyrus K Aidun and Jonathan R Clausen, “Lattice-boltzmann method for complex flows,” Annual Review of Fluid Mechanics, vol.

42, pp. 439–472, 2010.

[Anzt et al., 2012] H. Anzt, S. Tomov, M. Gates, J. Dongarra, and V. Heuveline, “Block-asynchronous multigrid smoothers for GPU-accelerated

systems,” Procedia Computer Science, vol. 9, pp. 7-16, 2012.

[Bertoglio et al., 2012] Bertoglio C, Moireau P, Gerbeau JF. Sequential parameter estimation for fluid–structure problems: Application to

hemodynamics. International Journal of Numerical Methods in Biomedical Engineering 2012; 28: 434–455, DOI: 10.1002/cnm.1476.

[Blanco et al., 2012] Blanco PJ, Watanabe SM, Feijo RA. Identification of vascular territory resistances in one-dimensional hemodynamics

simulations. Journal of Biomechanics 2012; 45: 2066–2073, DOI: 10.1016/j.jbiomech.2012.06.002.

[Calmac et al., 2015] Calmac, L., Niculescu, R., Badila, E., Weiss, E., Zamfir, D., Itu, L., Lazar, L., Carp, M., Itu, A., Suciu, C., Passerini, T., Sharma, P.,

Georgescu, B., Comaniciu, D., Image-Based Computation of Instantaneous Wave-free Ratio from Routine Coronary Angiography - Initial Validation

by Invasively Measured Coronary Pressures, TCT 2015, Journal of the Americal College of Cardiology, Vol. 66, 15_S.

[Calmac et al., 2016] Calmac, L., Niculescu, R., Badila, E., Weiss, E., Zamfir, D., Itu, L., Lazar, L., Carp, M., Itu, A., Suciu, C., Passerini, T., Sharma, P.,

Georgescu, B., Comaniciu, D., Image-Based Computation of Instantaneous Wave-free Ratio from Routine Coronary Angiography - Evaluation of a

Hybrid Decision Making Strategy with FFR, ACC 2016.

[Chen et al., 1998] Shiyi Chen and Gary D Doolen, “Lattice boltzmann method for fluid flows,” Annual review of fluid mechanics, vol. 30, no. 1, pp.

329–364, 1998.

[Connington et al., 2009] Kevin Connington, Qinjun Kang, Hari Viswanathan, Amr Abdel-Fattah, and Shiyi Chen, “Peristaltic particle transport using

the lattice boltzmann method,” Physics of Fluids (1994-present), vol. 21, no. 5, pp. 053301, 2009.

[Coogan et al., 2011] Coogan JS, Chan FP, Taylor CA, Feinstein JA. Computational fluid dynamic simulations of aortic coarctation comparing the

effects of surgical- and stent-based treatments on aortic compliance and ventricular workload. Catheterization and Cardiovascular Interventions

2011; 77: 680–691, DOI: 10.1002/ccd.22878.

Page 20: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

20

http://heart.unitbv.ro/

[Cousins et al., 2012] Cousins W, Gremaud PA. Boundary conditions for hemodynamics: The structured tree revisited. Journal of Computational

Physics 2012; 231: 6086–6096, DOI: 10.1016/j.jcp.2012.04.038.

[Cousins et al., 2013] Cousins W, Gremaud PA, Tartakovsky DM. A New Physiological Boundary Condition for Hemodynamics, SIAM Journal on

Applied Mathematics 2013; 73: 1203–1223, DOI: 10.1137/120895470.

[Cousins et al., 2014] Cousins W, Gremaud PA. Impedance boundary conditions for general transient hemodynamics, International Journal of

Numerical Methods in Biomedical Engineering 2014; 30: 1294–1313, DOI: 10.1002/cnm.2658.

[Formaggia et al., 2006] Formaggia L, Lamponi D, Tuveri M, Veneziani A. Numerical modeling of 1D arterial networks coupled with a lumped

parameters description of the heart. Computer Methods in Biomechanics and Biomedical Engineering 2006; 9: 273–88, DOI:

10.1080/10255840600857767.

[Frangi et al., 1999] Frangi, A.F., W.J. Niessen, R.M. Hoogeveen, T. van Walsum, M.A. Viergever, (1999), Model-based quantitation of 3-D magnetic

resonance angiographic images, IEEE Trans. Med. Imaging, 18 (10) pp. 946–956

[Geveler et al., 2013] M. Geveler, D. Ribbrock, D. Göddeke, P. Zajac, and S. Turek, “Towards a complete fem-based simulation toolkit on gpus:

Unstructured grid finite element geometric multigrid solvers with strong smoothers based on sparse approximate inverses,” Computers & Fluids,

vol. 80, pp. 327-332, 2013.

[Hutchins et al., 1976] Hutchins, G. M., Miner, M. M., Boitnott, J. K. “Vessel Caliber and Branch-Angle of Human Coronary Artery Branch-Points”,

Circulation Research, Vol. 38, pp. 572–576, 1976.

[Iberall, 1967] Iberall, A. S. “Anatomy and Steady Flow Characteristics of the Arterial System with an Introduction to its Pulsatile Characteristics”,

Mathematical Biosciences, Vol. 1, pp. 375–385, 1967.

[Ismail et al., 2013] Ismail M, Gee MW, Wall WA. CFD challenge: Hemodynamic simulation of a patient-specific aortic coarctation model with

adjoint-based calibrated windkessel elements. Lecture Notes in Computational Science 2013; 7746: 44–52, DOI: 10.1007/978-3-642-36961-2_6.

[Ismail et al., 2015] Ismail M, Wall WA, Gee MW. Adjoint-based inverse analysis of windkessel parameters for patient-specific vascular models.

Journal of Computational Physics 2015; 244: 113–130, DOI: 10.1016/j.jcp.2012.10.028.

[Itu et al., 2012] Itu L, Sharma P, Ralovich K, Mihalef V, Ionasec R, Everett A, Ringel R, Kamen A, Comaniciu D. Non-invasive hemodynamic

assessment of aortic coarctation: validation with in vivo measurements. Annals of Biomedical Engineering 2012; 41: 669–681, DOI:

10.1007/s10439-012-0715-0.

[Itu et al., 2015] Itu L, Sharma P, Passerini T, Kamen A, Suciu C, Comaniciu D. A parameter estimation framework for patient-specific hemodynamic

computations. Journal of Computational Physics 2015; 281: 316–333, DOI: 10.1016/j.jcp.2014.10.034.

[Kassab et al., 1995] Kassab, G., Fung, Y. “The Pattern of Coronary Arteriolar Bifurcations and the Uniform Shear Hypothesis”, Annals of

Biomechanical Engineering, Vol. 23, pp. 13–20, 1995.

[Keshavarz-Motamed et al., 2011] Keshavarz-Motamed Z, Garcia J, Pibarot P, Larose E, Kadem L. Modeling the impact of concomitant aortic

stenosis and coarctation of the aorta on left ventricular workload. Journal of Biomechanics 2011; 44: 2817–2825, DOI:

10.1016/j.jbiomech.2011.08.00.

[LaDisa et al., 2011] LaDisa JFJ, Figueroa CA, Vignon-Clementel IE, Kim HJ, Xiao N, Ellwein LM, Chan FP, Feinstein JA, Taylor CA. Computational

simulations for aortic coarctation: representative results from a sampling of patients. Journal of Biomechanical Engineering 2011; 133: 091008,

DOI: 10.1115/1.4004996.

[Malek et al., 2015] Malek J, Azar AT, Nasralli B, Tekari M, Kamoun H, Tourki R, Computational analysis of blood flow in the retinal arteries and

veins using fundus image. Computers & Mathematics with Applications 2015; 69: 101–116, DOI: 10.1016/j.camwa.2014.11.017.

[Murray, 1926] Murray, C. D. “The Physiological Principle of Minimum Work: I. The Vascular System and the Cost of Blood Volume”, Proc. of the

National Academy of Sciences of the USA, Vol. 12, pp. 207–214, 1926.

[Nita et al., 2014] C. Nita, Y. Chen, L. Lazar, V. Mihalef, L. M. Itu, M. Viceconti, C.Suciu, GPU Accelerated Finite Element Analysis of Trabecular Bone

Tissue, Proc. of the Virtual Physiological Human Conference, Trondheim, Norway, Sept 9-12, 2014

[Olufsen, 1998] Olufsen, M. S. “Modeling the Arterial System with Reference to an Anesthesia Simulator”, Ph.D. Thesis, 1998.

Page 21: HEARTheart.unitbv.ro/wp-content/uploads/2016/11/HEART-newsletter-octombrie-2016.pdfarborele coronarian (stâng sau drept) pe baza debitelor estimate la primul pas pentru fiecare ramură.

HEART Newsletter

Nr. 2 | Octombrie 2016

21

http://heart.unitbv.ro/

[Olufsen et al., 1999] Olufsen MS. Structured tree outflow condition for blood in the larger systemic arteries. The American Journal of Physiology

1999; 276: 257–268, DOI: 10.1114/1.1326031.

[Olufsen et al., 2000] Olufsen M, Peskin C. Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow

conditions. Annals of Biomedical Engineering 2000; 28: 1281-1299, DOI: 10.1114/1.1326031.

[Otsu, 1979] Otsu, N., (1979), A Threshold Selection Method from Gray-Level Histograms, IEEE Transactions on Systems, Man, and Cybernetics,

Vol. 9, No. 1, pp. 62-66.

[Petraco et al., 2012] R. Petraco et al. “Hybrid iFR-FFR Decision-Making Strategy: Implications for Enhancing Universal Adoption of Physiology-

Guided Coronary Revascularization”, EuroIntervention, vol. 8, pp. 1157-65, 2012.

[Sen et al., 2012] S. Sen et al., “Development and Validation of a New Adenosine-Independent Index of Stenosis Severity From Coronary Wave-

Intensity Analysis Results of the ADVISE Study”, J Am Coll Cardiol, vol. 59, pp. 1392-1402, 2012.

[Shapiro et al., 1969] Ascher H Shapiro, Michel Y Jaffrin, and Steven L Weinberg, “Peristaltic pumping with long wavelengths at low reynolds

number,” Journal of Fluid Mechanics, vol. 37, no. 04, pp. 799–825, 1969.

[Spilker et al., 2010] Spilker R, Taylor CA. Tuning multidomain hemodynamic simulations to match physiological measurements. Annals of

Biomedical Engineering 2010; 38: 2635–2648, DOI: 10.1007/s10439-010-0011-9.

[Stergiopulos et al., 1992] Stergiopulos N, Young DF, Rogge TR. Computer simulation of arterial flow with applications to arterial and aortic

stenosis. Journal of Biomechanics 1992; 25: 1477–1488; DOI: 10.1016/0021-9290(92)90060-E.

[Stroia et al., 2015] Stroia, I., Itu, L., Niţă, C, Lazăr, L., Suciu, C. GPU Accelerated Geometric Multigrid Method: Comparison with Preconditioned

Conjugate Gradient, 19th IEEE High Performance Extreme Computing Conference, Waltham, MA, USA, Sept. 15-17, 2015.

[Taylor, 1966] Taylor MG. Wave transmission through an assembly of randomly branching elastic tubes. Biophysical Journal 1966; 6: 697–716, DOI:

10.1016/S0006-3495(66)86689-4.

[Uylings, 1977] Uylings, H. “Optimization of Diameters and Bifurcation Angles in Lung and Vascular Tree Structures”, Bulletin of Mathematical

Biology, Vol. 39, pp. 509–520, 1977.

[Williams et al., 2012] S. Williams, D. D. Kalamkar, A. Singh, A. M. Deshpande, B. Van Straalen, M. Smelyanskiy, A. Almgren, P. Dubey, J. Shalf, and

L. Oliker, “Optimization of geometric multigrid for emerging multi-and manycore processors,” in Proceedings of the International Conference on

High Performance Computing, Networking, Storage and Analysis, p. 96, IEEE Computer Society Press, 2012.

[Yu et al., 2003] Dazhi Yu, Renwei Mei, Li-Shi Luo, and Wei Shyy, “Viscous flow computations with the method of lattice boltzmann equation,”

Progress in Aerospace Sciences, vol. 39, no. 5, pp. 329–367, 2003.

[***CFD challenge, 2014] ***, CFD challenge: predicting patient-specific hemodynamics at rest and stress through an aortic coarctation,

http://www.vascularmodel.org/miccai2013/, accessed 06/10/2014.

[***openMP, 2015] ***, http://openmp.org/wp/, accessed 06/10/2015.

[***vmtk, 2014] ***, The Vasular Modeling Toolkit. http://www.vmtk.org/, accessed 06/10/2014.