Universitatea din Craiova - old.unitbv.roold.unitbv.ro/Portals/31/Sustineri de...

Investeşte în oameni!

FONDUL SOCIAL EUROPEAN Programul Operaţional Sectorial Dezvoltarea Resurselor Umane 2007 – 2013 Axa prioritară 1 „Educaţie şi formare profesională în sprijinul creşterii economice şi dezvoltării societăţii bazate pe cunoaştere” Domeniul major de intervenţie 1.5. „Programe doctorale şi post-doctorale în sprijinul cercetării” Titlul proiectului: „Burse doctorale pentru dezvoltare durabila” BD-DD Numărul de identificare al contractului: POSDRU/107/1.5/S/76945 Beneficiar: Universitatea Transilvania din Braşov

Universitatea Transilvania din Braşov

Şcoala Doctorală Interdisciplinară

Departament: Automatică şi Tehnologia Informaţiei

Ing. Lucian Mihai ITU

Utilizarea procesării paralele în modelarea multiscalară a hemodinamicii coronariene

Parallel Processing in the Multiscale Modeling of Coronary Hemodynamics

Conducător ştiinţific

Prof. dr. ing. Florin MOLDOVEANU

BRAŞOV, 2013

MINISTERUL EDUCAŢIEI NAŢIONALE

UNIVERSITATEA “TRANSILVANIA” DIN BRAŞOV BRAŞOV, B-DUL EROILOR NR. 29, 500036, TEL. 0040-268-413000, FAX 0040-268-

410525 RECTORAT

D-lui (D-nei) ..............................................................................................................

COMPONENŢA Comisiei de doctorat

Numită prin ordinul Rectorului Universităţii „Transilvania” din Braşov Nr. 5971 din 25.07.2013

PREŞEDINTE: Conf. univ. dr. ing. Carmen GERIGAN

Universitatea “Transilvania” din Braşov CONDUCĂTOR ŞTIINŢIFIC: Prof. dr. ing. Florin MOLDOVEANU

Universitatea “Transilvania” din Braşov REFERENŢI: Prof. dr. ing. Dumitru POPESCU

Universitatea “Politehnica” din Bucureşti Prof. dr. ing. Dan ŞTEFĂNOIU

Universitatea “Politehnica” din Bucureşti Prof. dr. ing. DORU TALABĂ

Universitatea “Transilvania” din Braşov

Data, ora şi locul susţinerii publice a tezei de doctorat: 11.10.2013, ora 12, sala V III 9.

Eventualele aprecieri sau observaţii asupra conţinutului lucrării vă rugăm să

le transmiteţi în timp util, pe adresa [email protected] sau la numărul de fax: Departamentul de Automatică şi Tehnologia Informaţiei: 0268-418836

Totodată vă invităm să luaţi parte la şedinţa publică de susţinere a tezei de

doctorat.

Vă mulţumim.

mailto:[email protected]

III

CUPRINS

Pg. teză

Pg. rez.

1. INTRODUCERE ........................................................................................................1.1 Sistemul cardiovascular .........................................................................................1.2 Circulaţia coronariană............................................................................................1.3 Peretele arterial ......................................................................................................1.4 Sângele...................................................................................................................1.5 Simularea circuitului cardiovascular .....................................................................1.6 Patologii cardiovasculare.......................................................................................1.7 Obiectivele cercetării .............................................................................................1.8 Structura şi conţinutul tezei ...................................................................................

2. STADIUL ACTUAL ÎN SIMULAREA NUMERICĂ MULTISCALARĂ A HEMODINAMICII ARTERIALE ...........................................................................2.1 Modelarea multiscalară a circulaţiei arteriale........................................................

2.1.1 Modele tridimensionale..................................................................................2.1.2 Modele cvasi unidimensionale .......................................................................2.1.3 Modele cu parametrii distribuiţi .....................................................................2.1.4 Modele ale inimii ...........................................................................................2.1.5 Condiţii de frontieră utilizate în simulările hemodinamice............................

2.2 Simulări hemodinamice multiscalare cu geometrii arteriale specifice pacienţilor ..............................................................................................................

2.2.1 Imagistică prin rezonanţă magnetică ..............................................................2.2.2 Angiografia prin rezonanţă magnetică (MRA) ..............................................2.2.3 Angiografia prin tomografie computerizată ...................................................2.2.4 Angiografia prin raze X..................................................................................2.2.5 Ultrasunete .....................................................................................................

2.3 Procesare paralelă în modelarea hemodinamică multiscalară ...............................2.3.1 Introducere în procesarea paralelă..................................................................2.3.2 Utilizarea procesoarelor grafice pentru procesare paralelă ............................2.3.3 Utilizarea procesoarelor grafice pentru accelerarea simulărilor din

domeniul cardiovascular ................................................................................2.4 Concluzii................................................................................................................

3. MODELAREA MULTISCALARĂ DE ORDIN REDUS A HEMODINAMICII ARTERIALE ...........................................................................3.1 Modelul cvasi unidimensional...............................................................................

3.1.1 Joncţiuni arteriale ...........................................................................................3.1.2 Adimensionalizarea ecuaţiilor........................................................................

3.2 Condiţii de frontieră...............................................................................................3.2.1 Modelul Windkessel cu trei elemente ............................................................3.2.2 Arbore structurat ............................................................................................

3.3 Rezolvarea numerică a modelului hemodinamic unidimensional .........................3.3.1 Rezolvarea numerică a modelului unidimensional cu pereţi arteriali

elastici.............................................................................................................3.3.1.1 Metoda Lax-Wendroff ............................................................................3.3.1.2 Metoda descompunerii în serie Taylor....................................................

3.3.2 Rezolvarea numerică a modelului unidimensional cu pereţi arteriali vâscoelastici ...................................................................................................

3.4 Algoritmul de cuplare ............................................................................................

1 1 4 5 6 7 8 9 10

13 14 17 21 22 24 25 27 31 32 33 33 33 35 35 36 38 39

42 42 45 46 47 47 48 51 52 52 55 56 57

1 1 1 3 4

5 6 7

8 8 10 10 10 11 11 11 12

IV

3.5 Verificarea modelului unidimensional ..................................................................3.5.1 Curgerea Poisseuille .......................................................................................3.5.2 Comparaţie cu biblioteca LifeV .....................................................................

3.6 Rezultate ................................................................................................................3.6.1 Comparaţie între rezultatele obţinute cu diferite tipuri de condiţii de

frontieră de ieşire............................................................................................3.6.2 Influenţa algoritmului de cuplare asupra rezultatelor obţinute cu modele

vâscoelastice ale pereţilor arteriali .................................................................3.6.2.1 O singură arteră .......................................................................................3.6.2.2 Arbore arterial integral............................................................................

3.7 Concluzii................................................................................................................

4. PROCEDURĂ DE OPTIMIZARE EVOLUATĂ PENTRU PERSONALIZAREA SIMULĂRILOR HEMODINAMICE................................4.1 Introducere.............................................................................................................4.2 Metodologia procedurii de optimizare evoluate ....................................................

4.2.1 Estimarea şi compensarea rezistenţei .............................................................4.2.2 Estimarea şi compensarea complianţei ..........................................................

4.3 Rezultate ................................................................................................................4.3.1 Arteră carotidă comună idealizată cu raport fix al rezistenţei proxime .........4.3.2 Arteră carotidă comună idealizată cu valoare fixă a complianţei ..................4.3.3 Arteră carotidă comună idealizată cu complianţă adaptată indirect...............4.3.4 Bifurcaţie arterială iliacă idealizată ................................................................4.3.5 Geometrie aortică specifică unui pacient cu coarctaţie ..................................

4.4 Analiza rezultatelor................................................................................................4.5 Concluzii................................................................................................................

5. MODELAREA MULTISCALARĂ DE ORDIN REDUS A CIRCULAŢIEI CORONARIENE ŞI DIAGNOSTICAREA NON-INVAZIVĂ A STENOZELOR CORONARIENE ...........................................................................5.1 Introducere.............................................................................................................5.2 Model multiscalar al circulaţiei coronariene .........................................................

5.2.1 Model microvascular cu parametrii distribuiţi al circulaţiei coronariane ......5.2.2 Model de cădere de presiune pentru stenozele coronariene ...........................5.2.3 Model simplificat al inimii şi cuplarea cu circulaţia arterială ........................

5.3 Modelarea autoreglării coronariene în starea de repaus ........................................5.4 Personalizarea parametrilor modelului multiscalar ...............................................

5.4.1 Personalizarea parametrilor modelelor coronariene microvasculare .............5.4.1.1 Modelarea stării de repaus ......................................................................5.4.1.2 Modelarea stării de hiperemie.................................................................5.4.1.3 Estimarea rezistenţelor microvasculare de repaus ..................................5.4.1.4 Estimarea rezistenţelor microvasculare hiperemice................................

5.4.2 Personalizarea modelului inimii.....................................................................5.5 Procedură de optimizare pentru circulaţia coronariană .........................................

5.5.1 Optimizare cu modelul 0D detaliat ................................................................5.5.2 Optimizare cu modelul 1D cu geometrie redusă ............................................5.5.3 Optimizare cu modelul 1D cu geometrie completă........................................

5.6 Rezultate ................................................................................................................5.6.1 Testarea modelului simplificat al inimii.........................................................5.6.2 Validarea formei de undă a debitului .............................................................5.6.3 Validarea rezultatelor prin comparaţie cu ANSYS Fluent.............................5.6.4 Aplicarea procedurii de optimizare pentru geometria specifică unui

pacient ............................................................................................................5.6.5 Utilizarea modelului multiscalar pentru diagnosticarea pacienţilor cu

stenoze coronariene ........................................................................................

60 60 61 63 63 67 68 68 70

73 73 75 79 82 84 86 88 90 91 93 97 100

102 102 106 108 110 113 116 118 118 118 121 123 124 124 126 128 133 137 137 138 139 139 140 143

13 13 15 15 17

17 17 20 22 22 25 26

26 26 28 29 30 30 30 31 32 33 33 34 36 38 39 40 40

V

5.7 Concluzii................................................................................................................

6. DIAGNOSTICARE HEMODINAMICĂ NON-INVAZIVĂ A COARCTAŢIEI AORTICE......................................................................................6.1 Introducere.............................................................................................................6.2 Metodă non-invazivă de evaluare a coarctaţiilor aortice.......................................

6.2.1 Estimarea parametrilor pentru personalizarea modelului...............................6.2.2 Modelul căderii de presiune şi estimarea condiţiilor de frontieră ..................6.2.3 Estimare proprietăţilor pereţilor arteriali........................................................

6.3 Rezultate ................................................................................................................6.4 Analiza rezultatelor................................................................................................6.5 Concluzii................................................................................................................

7. UTILIZAREA UNUI PROCESOR GRAFIC PENTRU ACCELERAREA SIMULĂRILOR MULTISCALARE HEMODINAMICE.....................................7.1 Introducere.............................................................................................................7.2 Paralelizarea algoritmului numeric de rezolvare a modelului unidimensional .....

7.2.1 Algoritmul paralel hibrid CPU-GPU..............................................................7.2.2 Algoritmul paralel exclusiv GPU...................................................................

7.3 Rezultate ................................................................................................................7.3.1 Compararea execuţiei paralele şi secvenţiale cu diferite scheme numerice

de calcul..........................................................................................................7.3.2 Compararea strategiilor de copiere a datelor pentru algoritmul PHCG .........7.3.3 Compararea performanţelor obţinute cu algoritmii SCO, PHCG şi PGO......

7.4 Concluzii................................................................................................................

8. CONCLUZII FINALE...............................................................................................8.1 Concluzii................................................................................................................8.2 Contribuţii originale...............................................................................................

8.2.1 Modele cardiovasculare..................................................................................8.2.2 Diagnosticarea non-invazivă a patologiilor cardiovasculare..........................8.2.3 Proceduri de optimizare pentru personalizarea simulărilor hemodinamice ...8.2.4 Accelerarea simulărilor hemodinamice prin utilizarea unui procesor

grafic...............................................................................................................8.3 Diseminarea rezultatelor........................................................................................8.4 Direcţii viitoare de cercetare..................................................................................

A1. ADAPTAREA AUTOMATĂ A CONDIŢIEI DE FRONTIERĂ DE TIP ARBORE STRUCTURAT ........................................................................................

A2. ESTIMAREA PRESIUNII ARTERIALE VARIABILE ÎN TIMP FOLOSIND DATE ACHIZIŢIONATE PRIN REZONANŢĂ MAGNETICĂ ...

A3. IMPLEMENTAREA OPTIMIZATĂ PE UN PROCESOR GRAFIC A ALGORITMILOR BAZAŢI PE ŞABLOANE........................................................

A4. IMPLEMENTAREA OPTIMIZATĂ PE UN PROCESOR GRAFIC A ECUAŢIILOR NAVIER-STOKES ..........................................................................

BIBLIOGRAFIE ............................................................................................................

LISTĂ DE FIGURI........................................................................................................

LISTĂ DE TABELE ŞI ALGORITMI ........................................................................

LISTĂ DE ABREVIERI................................................................................................

LISTĂ DE NOTAŢII .....................................................................................................

REZUMAT .....................................................................................................................

ABSTRACT ....................................................................................................................

146

148 148 151 152 153 156 157 161 166

168 168 170 171 177 178 178 179 179 183

187 187 190 190 191 191 192 192 194

196

204

208

223

235

251

257

260

263

270

271

42

42 42 42 44 44 46 47 49 50

50 50 51 55 56 56 57 59

59 59 61 61 62 62 63 63

65

68

68

VI

CONTENTS....................................................................................................................

CURRICULUM VITAE................................................................................................

272

276

69

VII

CONTENTS

Pg. teză

Pg. rez.

1. INTRODUCTION......................................................................................................1.1 The cardiovascular system .....................................................................................1.2 The coronary circulation ........................................................................................1.3 The arterial wall .....................................................................................................1.4 Blood ......................................................................................................................1.5 Simulation of the cardiovascular system................................................................1.6 Cardiovascular pathologies ....................................................................................1.7 Research objectives ................................................................................................1.8 Structure and content of the thesis .........................................................................

2. STATE-OF-THE-ART MULTISCALE NUMERICAL COMPUTATIONS OF ARTERIAL HEMODYNAMICS.......................................................................2.1 Multiscale modeling of the arterial circulation ......................................................

2.1.1 Tridimensional models....................................................................................2.1.2 Quasi one-dimensional models .......................................................................2.1.3 Lumped parameter models ..............................................................................2.1.4 Heart models ...................................................................................................2.1.5 Boundary conditions in hemodynamic computations .....................................

2.2 Multiscale hemodynamic computations for patient-specific arterial geometries...2.2.1 Magnetic resonance imaging...........................................................................2.2.2 Magnetic resonance angiography....................................................................2.2.3 Computer tomography angiography................................................................2.2.4 X-ray angiography...........................................................................................2.2.5 Ultrasound .......................................................................................................

2.3 Parallel Processing for multiscale hemodynamic modeling ..................................2.3.1 Introduction in parallel processing..................................................................2.3.2 Graphics Processing Units for parallel processing..........................................2.3.3 Graphics Processing Units for the acceleration of cardiovascular

computations ..................................................................................................2.4 Conclusions ............................................................................................................

3. REDUCED-ORDER MULTISCALE MODELING OF ARTERIAL HEMODYNAMICS ...................................................................................................3.1 Quasi one-dimensional model ................................................................................

3.1.1 Arterial junctions.............................................................................................3.1.2 Non-dimensional equations.............................................................................

3.2 Boundary conditions ..............................................................................................3.2.1 Three-element Windkessel model...................................................................3.2.2 Structured tree .................................................................................................

3.3 Numerical solution of the hemodynamic one-dimensional model.........................3.3.1 Numerical solution of the one-dimensional model with elastic arterial

walls ...............................................................................................................3.3.1.1 Lax-Wendroff method..............................................................................3.3.1.2 Taylor series expansion method...............................................................

3.3.2 Numerical solution of the one-dimensional model with viscoelastic arterial walls ...................................................................................................

3.4 Coupling algorithm ................................................................................................3.5 Verification of the numerical solution of the one-dimensional model...................

1 1 4 5 6 7 8 9 10

13 14 17 21 22 24 25 27 31 32 33 33 33 35 35 36 38 39

42 42 45 46 47 47 48 51 52 52 55 56 57 60

1 1 1 3 4

5 6 7

8 8 10 10 10 11 11 11 12

VIII

3.5.1 Poisseuille flow ...............................................................................................3.5.2 Comparison with the LifeV library .................................................................

3.6 Results ....................................................................................................................3.6.1 Comparison of results obtained with different types of boundary

conditions .......................................................................................................3.6.2 Effect of coupling algorithm on the results obtained with viscoelastic

arterial wall models ........................................................................................3.6.2.1 Single artery .............................................................................................3.6.2.2 Full body arterial tree ...............................................................................

3.7 Conclusions ............................................................................................................

4. ADVANCED TUNING PROCEDURE FOR THE PERSONALIZATION OF HEMODYNAMIC COMPUTATIONS....................................................................4.1 Introduction ............................................................................................................4.2 Advanced Tuning Procedure..................................................................................

4.2.1 Resistance estimation and compensation ........................................................4.2.2 Compliance estimation and compensation......................................................

4.3 Results ....................................................................................................................4.3.1 Idealized common carotid artery model with chosen proximal resistance

fraction ...........................................................................................................4.3.2 Idealized common carotid artery model with chosen compliance ..................4.3.3 Idealized common carotid artery model with indirectly adapted

compliance .....................................................................................................4.3.4 Idealized iliac bifurcation................................................................................4.3.5 Patient-specific aortic geometry with coarctation...........................................

4.4 Discussion ..............................................................................................................4.5 Conclusion..............................................................................................................

5. REDUCED-ORDER MULTISCALE MODELING OF THE CORONARY CIRCULATION AND NON-INVASIVE ASSESSMENT OF CORONARY STENOSES .................................................................................................................5.1 Introduction ............................................................................................................5.2 Multiscale model of the coronary circulation.........................................................

5.2.1 Lumped parameter microvascular model of the coronary circulation ............5.2.2 Pressure-drop model for coronary stenoses ....................................................5.2.3 Lumped parameter heart model and coupling with arterial circulation ..........

5.3 Coronary rest state autoregulation modeling..........................................................5.4 Personalization of the parameters of the multiscale model....................................

5.4.1 Personalization of the parameters of the coronary microvascular models......5.4.1.1 Rest state modeling ..................................................................................5.4.1.2 Hyperemia modeling ................................................................................5.4.1.3 Estimation of rest state microvascular resistances ...................................5.4.1.4 Estimation of hyperemic microvascular resistances ................................

5.4.2 Personalization of the heart model ..................................................................5.5 Tuning procedure for the coronary circulation ......................................................

5.5.1 Tuning with the detailed lumped parameter model.........................................5.5.2 Tuning with the reduced-geometry one-dimensional model...........................5.5.3 Tuning with the full geometry one-dimensional model ..................................

5.6 Results ....................................................................................................................5.6.1 Lumped parameter heart model.......................................................................5.6.2 Validation of the coronary flow rate waveform..............................................5.6.3 Validation through comparison with ANSYS Fluent .....................................5.6.4 Tuning procedure results for a patient-specific geometry...............................5.6.5 Applying the multiscale model for the diagnosis of patients with coronary

stenoses...........................................................................................................

60 61 63 63 67 68 68 70

73 73 75 79 82 84 86 88 90 91 93 97 100

102 102 106 108 110 113 116 118 118 118 121 123 124 124 126 128 133 137 137 138 139 139 140 143

13 13 15 15 17

17 17 20 22 22 25 26

26 26 28 29 30 30 30 31 32 33 33 34 36 38 39 40 40

IX

5.7 Conclusions ............................................................................................................

6. NON-INVASIVE HEMODYNAMIC ASSESSMENT OF AORTIC COARCTATIONS .....................................................................................................6.1 Introduction ............................................................................................................6.2 Method for non-invasive assessment of aortic coarctations...................................

6.2.1 Parameter estimation for model personalization.............................................6.2.2 Pressure drop model and estimation of boundary conditions .........................6.2.3 Estimation of arterial wall properties ..............................................................

6.3 Results ....................................................................................................................6.4 Discussion ..............................................................................................................6.5 Conclusions ............................................................................................................

7 GRAPHICS PROCESSING UNIT ACCELERATED MULTISCALE HEMODYNAMIC COMPUTATIONS....................................................................7.1 Introduction ............................................................................................................7.2 Parallelization of the numerical solution scheme of the one-dimensional model..

7.2.1 Parallel hybrid CPU-GPU algorithm...............................................................7.2.2 Parallel GPU only algorithm...........................................................................

7.3 Results ....................................................................................................................7.3.1 Comparison of parallel and sequential execution for different numerical

schemes ..........................................................................................................7.3.2 Comparison of the memory copy strategies for the PHCG algorithm ............7.3.3 Comparison of the performance obtained with the SCO, MCO, PHCG and

PGO algorithms..............................................................................................7.4 Conclusions ............................................................................................................

8. FINAL CONCLUSIONS ...........................................................................................8.1 Conclusions ............................................................................................................8.2 Original contributions ............................................................................................

8.2.1 Cardiovascular models ....................................................................................8.2.2 Non-invasive assessment of cardiovascular pathologies ................................8.2.3 Tuning procedures for the personalization of hemodynamic computations ...8.2.4 Graphics processing unit based acceleration of hemodynamic

computations ..................................................................................................8.3 Dissemination of results .........................................................................................8.4 Future research activities........................................................................................

A1. TUNING PROCEDURE FOR STRUCTURED TREE BOUNDARY CONDITION...............................................................................................................

A2. TIME-VARYING ARTERIAL PRESSURE ESTIMATION USING DATA FROM MAGNETIC RESONANCE IMAGING ....................................................

A3. GRAPHICS PROCESSING UNIT OPTIMIZED COMPUTATION OF STENCIL BASED ALGORITHMS .........................................................................

A4. GRAPHICS PROCESSING UNIT OPTIMIZED COMPUTATION OF THE NAVIER-STOKES EQUATIONS...................................................................

REFERENCES...............................................................................................................

LIST OF FIGURES .......................................................................................................

LIST OF TABLES AND ALGORITHMS...................................................................

LIST OF ABBREVIATIONS........................................................................................

LIST OF NOTATIONS .................................................................................................

ABSTRACT ....................................................................................................................

146

148 148 151 152 153 156 157 161 166

168 168 170 171 177 178 178 179 179 183

187 187 190 190 191 191 192 192 194

196

204

208

223

235

251

257

260

263

270

42

42 42 42 44 44 46 47 49 50

50 50 51 55 56 56 57 59

59 59 61 61 62 62 63 63

65

68

X

CONTENTS....................................................................................................................

CURRICULUM VITAE................................................................................................

272

276

69


1

1. INTRODUCERE

1.1 SISTEMUL CARDIOVASCULAR

Sistemul cardiovascular transportă oxigen şi nutrimente către toate ţesuturile corpului, de unde preia dioxidul de carbon şi alte substanţe rezultate ca urmare a metabolismului celulelor. Din punct de vedere fizic sistemul constă din două pompe paralele sincronizate (partea stângă şi partea dreaptă a inimii), care propulsează un fluid vâscos (sângele) printr-o reţea de tuburi flexibile. Partea dreaptă a inimii pompează sânge deoxigenat prin circulaţia pulmonară, iar partea stângă a inimii pompează sânge oxigenat prin circulaţia sistemică. Figura 1.1 prezintă arterele şi venele principale ale sistemului cardiovascular.

Inima generează energia necesară transportării sângelui prin sistemul cardiovascular. Aceasta constă din patru cavităţi: două ventricule şi două atriumuri (v. fig. 1.2), ale căror dimensiuni variază de-a lungul ciclului cardiac datorită activităţii muşchiului inimii (miocardul). Atriumul drept primeşte sânge de la venele sistemice, acesta curgând mai departe prin ventriculul drept către vasele pulmonare. Din plămâni, sângele este direcţionat prin intermediul venelor pulmonare către atriumul stâng şi apoi către ventriculul stâng, de unde este pompat prin aortă şi transportat în restul corpului prin intermediul circulaţiei sistemice. Există patru valve în interiorul inimii, câte una la ieşirea fiecărei cavităţi, care controlează curgerea sângelui prin inimă şi asigură curgerea preponderent unidirecţională a sângelui prin sistem.

1.6 PATOLOGII CARDIOVASCULARE

Patologiile cardiovasculare reprezintă o clasă de afecţiuni care implică inima, arterele sau venele. În continuare, sunt descrise câteva dintre problemele cardiovasculare principale împreună cu posibilele tratamente, cu scopul de a arăta modul în care simulările numerice ale acestor stări patologice pot ajuta medicii în procesul de diagnosticare şi de luare a deciziilor referitoare la tratament [Nichols et al., 2011]: ateroscleroza: este o condiţie patologică în care peretele arterial se îngroaşă datorită

acumulării de grăsime (colesterol). Aceasta apare prin formarea unor plăci aterosclerotice multiple, care pot creşte şi care într-un stadiu avansat pot genera o stenoză severă sau ocluzia totală a vasului de sânge. Suplimentar, ruptura depunerilor conduce la expunerea circulaţiei la un material trombogenic, precum colagenul, care poate conduce la apariţia trombozei. Astfel, se pot bloca direct arterele, precum arterele coronariene iar, de cele mai multe ori se ajunge la o blocare a microcirculaţiei şi, ca urmare, la apariţia tromboembolismului. Procedurile medicale principale pentru tratarea arterelor stenozate sunt crearea unui bypass arterial sau plasarea unui stent (un dispozitiv care sprijină peretele arterial, menţinând aria originală a lumenului);

anevrismele: sunt dilatări locale graduale ale arterelor. Acestea apar, în principal, în arterele de la baza creierului, e.g. în cercul lui Willis, sau în aortă. Când dimensiunea unui anevrism creşte, există un risc sporit de ruptură, rezultând o hemoragie severă, alte complicaţii sau chiar moartea pacientului. Un anevrism poate fi fuziform sau sacular. În primul caz, de obicei, se plasează stenturi, iar în cel de al doilea caz anevrismul este tăiat sau umplut cu un fir metalic subţire. În ambele situaţii, simulările numerice pot ajuta medicul în vederea scăderii riscurilor post-operatorii, precum recanalizarea;

defecte congenitale: sunt condiţii patologice care sunt prezente de la naştere, sau care se dezvoltă în primele luni ale vieţii. Un defect congenital poate fi rezultatul unei anormalităţi genetice sau cromozomale, al mediului intrauterin, al unei erori de morfogeneză sau al unei infecţii. Un exemplu de astfel de defect este coarctaţia aortică. În general, defectele congenitale necesită o intervenţie chirurgicală, precum procedura Fontan în cazul atreziei


tricuspidiene [Fontan et al., 1971]. În acest caz, simulările numerice pot prezice aspectele hemodinamice principale, precum presiunea şi debitul, înainte şi după intervenţie;

insuficienţa cardiacă: reprezintă incapacitatea inimii de a furniza un debit cardiac suficient alimentării corpului. Poate conduce la numeroase simptome precum creşterea frecvenţei respiratorii, umflarea picioarelor şi intoleranţă la exerciţiu fizic. În cadrul tezei se vor studia două patologii diferite, şi anume ateroscleroza coronariană şi

coarctaţia aortică.

Fig. 1.1 Sistemul cardiovascular uman – arterele şi venele principale [***Wiki, 2012(a)].

2


3

Fig. 1.2 Componentele principale ale inmii [***Wiki, 2012(a)].

1.7 OBIECTIVELE CERCETĂRII

Afecţiunile cardiovasculare reprezintă principala cauză de deces la nivel mondial. Dintre multitudinea de boli cardiovasculare, afecţiunile circulaţiei coronariene reprezintă cauza a aproape 50% dintre decese [***WHO, 2011]. În România peste 60.000 de persoane decedează anual datorită problemelor cardiace. În ciuda îmbunătăţirilor semnificative ale tehnicilor de imagistică medicală şi a altor modalităţi de diagnosticare, rata de mortalitate rămâne ridicată, cauza principală fiind lipsa unor indicatori de diagnosticare in-vivo şi in-vitro personalizaţi. De exemplu, în cazul stenozelor coronariene, evaluarea pur anatomică a acestora poate conduce la o diagnosticare eronată a severităţii îngustărilor. De aceea, este necesară o evaluare funcţională a severităţii stenozelor. Evaluarea funcţională este însă o procedură invazivă care necesită utilizarea unui senzor de presiune şi inducerea stării de hiperemie coronariană. Această procedură este, pe de o parte, foarte costisitoare iar, pe de altă parte, prezintă şi o serie de riscuri pentru pacient.

Prin urmare, activitatea de cercetare ce a făcut obiectul prezentei teze de doctorat s-a centrat în jurul următorului obiectiv principal: studiul, dezvoltarea, implementarea, testarea şi validarea unui model multiscalar al circulaţiei coronariene pentru diagnosticarea non-invazivă a stenozelor coronariene. Având în vedere caracterul general al modelului multiscalar dezvoltat, obiectivul secundar al tezei a fost acela de a se dezvolta, implementa, testa şi valida un model multiscalar pentru diagnosticarea non-invazivă a coarctaţiilor aortice.

Având în vedere aceste două obiective, au fost formulate o serie de obiective specifice: studiul bibliografic al cercetărilor teoretice şi experimentale privind modalităţile de

modelare hemodinamică a circulaţiei arteriale în general şi a circulaţiei coronariene în particular;

studiul bibliografic al cercetărilor teoretice şi experimentale privind implementarea pe procesoare grafice a algoritmilor numerici utilizaţi în modelarea hemodinamică;

dezvoltarea, implementarea, personalizarea şi testarea unui model simplificat al inimii; dezvoltarea şi implementarea unui model multiscalar de ordin redus pentru simularea

hemodinamicii arteriale; dezvoltarea şi implementarea unui model multiscalar de ordin redus pentru simularea

hemodinamicii coronariene;


4

dezvoltarea, implementarea şi testarea unui algoritm de cuplare a segmentelor vasculare pentru cazul în care se foloseşte un model vâscoelastic al peretelui arterial;

dezvoltarea, implementarea şi testarea unui model de autoreglare a circulaţiei coronariene pentru starea de repaus a pacientului;

dezvoltarea, implementarea şi testarea unui algoritm de personalizare a simulărilor hemodinamice ale circulaţiei coronariene pentru stările de repaus şi de hiperemie ale pacienţilor;

dezvoltarea, implementarea şi testarea unui flux de lucru automatizat pentru efectuarea simulărilor hemodinamice coronariene personalizate corespunzătoare pacienţilor;

validarea simulărilor hemodinamice corespunzătoare circulaţiei coronariene prin intermediul indicatorului diagnostic FFR (Fractional Flow Reserve);

dezvoltarea, implementarea şi testarea unei proceduri de optimizare evoluate pentru personalizarea simulărilor hemodinamice;

dezvoltarea şi implementarea unui model multiscalar de ordin redus pentru simularea hemodinamicii arteriale a pacienţilor cu coarctaţie;

validarea simulărilor hemodinamice ale pacienţilor cu coarctaţie; implementarea şi testarea pe un procesor grafic a modelului multiscalar al circulaţiei

arteriale în vederea accelerării timpului de execuţie.

1.8 STRUCTURA ŞI CONŢINUTUL TEZEI

Teza de doctorat este formată din nouă capitole şi patru anexe. În capitolul 1, Introducere, sunt prezentate noţiuni introductive în ceea ce priveşte sistemul

cardiovascular, respectiv principalele componente ale acestuia (inima, circulaţia coronariană, pereţii arteriali, sângele), simularea numerică a circuitului cardiovascular şi principalele tipuri de patologii ale acestuia. De asemenea, sunt definite principalele obiective ale tezei de doctorat.

În capitolul 2, Stadiul actual în simularea numerică multiscalară a hemodinamicii arteriale, se realizează iniţial o introducere în domeniul modelelor hemodinamice multiscalare şi se prezintă componentele principale ale unui astfel de model. Deoarece, scopul principal al modelării hemodinamice este acela de a se realiza simulări specifice pacienţilor, este necesar ca geometria arterială a pacienţilor să fie achiziţionată prin tehnici de imagistică medicală. Se introduc principalele tehnologii utilizate în acest sens precum şi lucrările de referinţă din literatura de specialitate în care se prezintă simulări specifice pacienţilor. De asemenea, se introduc aspecte de procesare paralelă, cu accent pe procesarea paralelă pe procesoare grafice.

Având în vedere faptul că scopul final al dezvoltării simulatoarelor hemodinamicii arteriale este ca acestea să fie folosite în clinici de specialitate pentru diagnosticarea pacienţilor, un timp de execuţie redus este crucial. În acest sens, în capitolul 3, Modelarea multiscalară de ordin redus a hemodinamicii arteriale, se introduce un model de ordin redus care poate fi utilizat pentru a se analiza principalele caracteristici ale circulaţiei arteriale (debite şi presiuni). Se introduce modelul cvasi unidimensional precum şi metodele numerice de rezolvare a acestuia, dar şi o serie de condiţii de frontieră fiziologice care pot fi aplicate la frontierele de ieşire ale arborelui arterial.

În capitolul 4, Procedură de optimizare evoluată pentru personalizarea simulărilor hemodinamice, se introduce un nou algoritm de optimizare care poate fi utilizat în cadrul simulărilor hemodinamice cu scopul ca anumite mărimi din cadrul simulării (presiuni, debite) să aibă aceleaşi valori ca şi mărimile măsurate la pacienţi. Prin intermediul procedurii de optimizare propuse numărul de rulări repetitive ale unei simulări hemodinamice pentru un singur pacient este redus considerabil şi, de asemenea, riscul ca procedura de optimizare să nu găsească o soluţie pentru parametrii adaptaţi este diminuat.

În capitolul 5, Modelarea multiscalară de ordin redus a circulaţiei coronariene şi diagnosticarea non-invazivă a stenozelor coronariene, atenţia este îndreptată asupra circulaţiei coronariene. Deoarece scopul modelului multiscalar dezvoltat este să se analizeze circulaţia


coronariană în starea patologică (în prezenţa stenozelor), acesta este completat cu un model de cădere de presiune. Pentru a utiliza acest model multiscalar al circulaţiei coronariene pentru geometrii coronariene achiziţionate de la pacienţi, este necesară o personalizare a modelului. Astfel, se prezintă personalizarea modelului inimii, dar şi personalizarea condiţiilor de frontieră de ieşire aplicate în stările de repaus şi de hiperemie ale pacientului. În finalul capitolului se prezintă rezultatele simulărilor hemodinamice realizate pentru o serie de pacienţi. Aceste rezultate sunt validate prin intermediul indicatorului diagnostic FFR.

5

În capitolul 6, Diagnosticare hemodinamică non-invazivă a coarctaţiei aortice, se introduce un model multiscalar pentru evaluarea non-invazivă a pacienţilor cu coarctaţie aortică, atât pentru cazul pre-operatoriu, cât şi pentru cazul post-operatoriu. Componentele principale ale metodei propuse sunt modelul de ordin redus introdus în capitolul 3 şi adaptat la geometrii aortice specifice pacienţilor, o procedură de estimare a parametrilor dezvoltată pentru determinarea condiţiilor de frontieră personalizate şi a proprietăţilor pereţilor arteriali pe baza măsurătorilor non-invazive, precum şi un model comprehensiv de cădere de presiune cuplat cu modelul de ordin redus.

În capitolul 7, Utilizarea unui procesor grafic pentru accelerarea simulărilor multiscalare hemodinamice, se prezintă o implementare numerică bazată pe procesor grafic (GPU – Graphical Processing Unit) a modelului unidimensional, realizată în vederea reducerii timpului de execuţie. Au fost dezvoltaţi doi algoritmi, un algoritm hibrid CPU-GPU şi un algoritm exclusiv GPU, care sunt comparaţi cu algoritmul secvenţial bazat pe CPU. S-au evaluat diverse strategii de copiere a datelor între CPU şi GPU, conducând la trei variante ale algoritmului hibrid CPU-GPU.

În capitolul 8, Concluzii finale, sunt prezentate concluziile finale, contribuţiile personale din cadrul tezei de doctorat, diseminarea rezultatelor, precum şi direcţiile viitoare de cercetare.

2. STADIUL ACTUAL ÎN SIMULAREA NUMERICĂ MULTISCALARĂ A HEMODINAMICII ARTERIALE

Modelarea hemodinamicii arteriale are la bază teoria mecanicii fluidelor. Dinamica fluidelor computaţională (CFD) este, aşa cum s-a mai precizat o ramură a mecanicii fluidelor care foloseşte metode şi algoritmi numerici pentru rezolvarea şi analiza problemelor care implică curgeri de fluide [Ferziger et al., 2002]. Ecuaţiile fundamentale, care guvernează orice curgere de fluid sunt date de sistemul generalizat al ecuaţiilor Navier-Stokes (ecuaţia de conservare a masei, ecuaţiile de conservare a impulsului şi ecuaţia de conservare a energiei). Rezolvarea analitică a acestui sistem de ecuaţii este posibilă numai pentru cazuri idealizate, în care o serie de ipoteze simplificatoare conduc la o formă redusă a ecuaţiilor. De aceea, în general, rezolvarea sistemului de ecuaţii Navier-Stokes este realizat exclusiv sub formă numerică.

De-a lungul timpului au fost dezvoltate trei tipuri de metode numerice pentru rezolvarea acestor ecuaţii, şi anume: metoda diferenţelor finite, metoda elementelor finite şi metoda volumelor finite [Wendt, 2009].

În cazul aplicaţiilor hemodinamice, în care s-a avut în vedere simularea curgerii sângelui prin sistemul cardiovascular, sistemul de ecuaţii Navier-Stokes se rezolvă în interiorul unui domeniu închis, delimitat de pereţii vaselor de sânge. Incompresibilitatea sângelui, conduce la o simplificare a ecuaţiilor Navier-Stokes care, de obicei, se rezolvă fără a se considera ecuaţia de conservare a energiei (deoarece, în general, se consideră cazul izoterm) [Wendt, 2009]. Introducerea condiţiei de incompresibilitate duce însă şi la o dificultate suplimentară, şi anume la dispariţia termenului dependent de timp din ecuaţia de conservare a masei. Pentru rezolvarea acestei probleme au fost dezvoltate o serie de metode specializate, dintre care cea mai utilizată, la momentul actual, este metoda proiecţiei [Chung, 2002].

Simulările hemodinamice se concentrează, în general, pe un segment particular al circuitului cardiovascular, care este modelat tridimensional folosindu-se datele achiziţionate din imaginile preluate prin rezonanţă magnetică sau tomografie computerizată [Migliavacca et al.,


1999(b)]. Datorită timpului de execuţie foarte mare al simulărilor, acestea pot fi realizate doar pe segmente restrânse ale circuitului cardiovascular. Deoarece circuitul cardiovascular este un circuit închis, simulările realizate doar pentru un segment particular, fără a se ţine cont de influenţa restului circuitului asupra acestui segment, conduc la rezultate eronate, mai ales în ceea ce priveşte valorile de presiune [Quarteroni et al., 2001].

6

Pentru a rezolva această problemă au fost introduse modelele multiscalare, care se bazează pe o modelare tridimensională detaliată a segmentului de interes şi pe o modelare unidimensională sau cu parametrii distribuiţi pentru restul circuitului. În acest fel se poate ţine cont de influenţa restului sistemului circulatoriu asupra segmentului local [Formaggia et al., 2001].

Deoarece simularea numerică a hemodinamicii arteriale necesită un timp de execuţie foarte mare, există o cerere continuă de viteze de calcul tot mai mari. Procesoarele paralele conduc la o îmbunătăţire a vitezei de execuţie, deoarece majoritatea algoritmilor din dinamica fluidelor pot fi paralelizaţi relativ uşor dar, de multe ori, costul unor astfel de centre de calcul (clustere) este foarte mare [Zou et al., 2009]. Odată cu dezvoltarea procesoarelor grafice GPU, care conţin un număr foarte mare de core-uri, calculele pot fi efectuate cu viteze semnificativ mai mari şi, în special, de la introducerea limbajului CUDA (Compute Unified Device Architecture) de către NVIDIA în 2006, efortul de programare a unui GPU, folosit pentru calcule matematice comune, a fost redus semnificativ [Jesperson, 2009].

2.1 MODELAREA MULTISCALARĂ A CIRCULAŢIEI ARTERIALE

În ultimii ani, progresele realizate în domeniul simulărilor numerice ale dinamicii fluidelor, împreună cu dezvoltarea unor noi tehnologii în ceea ce priveşte achiziţia datelor experimentale, au permis obţinerea unor informaţii cantitative asupra curgerii sângelui. Cercetările s-au concentrat mai ales asupra determinării valorilor de viteză şi presiune, pornind de la mărimi globale medii, precum debitul sau presiunea, specificate la frontierele compartimentului vascular supus analizei [Vignon-Clementel et al., 2010].

În artere mari şi medii, sângele se comportă ca un fluid newtonian incompresibil. În aceste condiţii, ecuaţiile Navier-Stokes pentru cazul incompresibil sunt capabile să furnizeze viteza şi presiunea sângelui [Li, 2004].

Una dintre dificultăţile majore ale modelării precise a sistemului cardiovascular uman este faptul că acesta este un circuit închis cu un grad de interdependenţă mare între componentele individuale. Caracteristicile de curgere a sângelui într-un anumit segment al sistemului (hemodinamica locală) sunt strâns corelate cu dinamica globală a sistemului [Vignon-Clementel et al., 2006]. Distribuţia curgerii sângelui în diverse segmente vasculare este o caracteristică a întregului sistem, dar ea influenţează dinamica fiecărui segment. Studiul curgerii locale a sângelui este important deoarece anumite patologii precum îngroşarea locală a vasului de sânge sau formarea de stenoze sunt puternic influenţate de hemodinamica locală. Pe de altă parte, anumite modificări locale, precum schimbarea lumenului vascular, pot duce la redistribuirea globală a curgerii sângelui, declanşând un mecanism compensatoriu care asigură un debit suficient în zona de aval a stenozei.

Această influenţă reciprocă între hemodinamica globală sau sistemică şi cea locală a condus la conceptul de modelare geometrică multiscalară a circulaţiei.

Modelul local se bazează, de obicei, pe soluţia sistemului de ecuaţii Navier-Stokes (tridimensionale) pentru cazul fluidelor incompresibile. Pentru modelul simplificat al întregului circuit vascular există două posibilităţi: modele cvasi unidimensionale (modele 1D); modele cu parametri distribuiţi (modele 0D).

Modelele unidimensionale folosesc o serie de ipoteze simplificatoare pentru a transforma un domeniu tridimensional într-unul unidimensional de-a lungul axei longitudinale a vasului de sânge, rezultând, de obicei, un sistem format din trei ecuaţii cu derivate parţiale. Un aspect fundamental al acestor modele este faptul că se ţine cont de complianţa vaselor, ceea ce permite


descrierea fenomenelor de undă care apar în sistemul cardiovascular. Modelele 0D se bazează pe rezolvarea în domeniul timpului a unui sistem de ecuaţii diferenţiale ordinare pentru a determina valorile medii ale debitului şi presiunii în diverse compartimente ale sistemului cardiovascular.

7

Ambele modele au o importanţă deosebită pentru abordarea multiscalară deoarece permit reprezentarea fenomenelor principale sau fundamentale ale circulaţiei (precum mecanismul compensatoriu descris mai sus) într-un interval de timp rezonabil.

O problemă de mare interes este cuplarea modelelor simplificate (1D sau 0D) cu modelul local detaliat (3D) [Formaggia et al., 2001], [Quarteroni et al., 2001], un exemplu fiind prezentat în figura 2.1.

Fig. 2.1 Cuplarea unor modele tridimensionale, unidimensionale şi cu parametrii distribuiţi

(0D) pentru realizarea unei simulări multiscalare cu interacţiune fluid-structură [van de Vosse et al., 2011].

2.4 CONCLUZII

Principalul scop al simulărilor hemodinamice este acela de a realiza simulări specifice pacienţilor pentru a se putea înţelege caracteristicile hemodinamice critice în diverse situaţii patologice. Simulările hemodinamice pot astfel furniza doctorilor informaţii suplimentare care îi pot ajuta pe aceştia în stabilirea diagnosticului sau în alegerea tipului tratamentului, etc. Pentru a putea realiza astfel de simulări specifice pacienţilor este necesar să se achiziţioneze geometria arterială printr-o tehnică de imagistică medicală. După o scurtă trecere în revistă a lucrărilor de referinţă din literatura de specialitate, se prezintă principalele tehnologii de imagistică medicală: rezonanţa magnetică nucleară, angiografia prin rezonanţă magnetică, angiografia prin tomografie computerizată, angiografia prin raze X şi ultrasunete.

Chiar şi în cazul utilizării modelelor multiscalare, unde doar zona de interes este reprezentată în detaliu, timpii de execuţie necesari pentru rularea acestor modele se întinde de la câteva ore până la câteva zile. Evident, într-un cadru clinic, un astfel de timp de execuţie este inacceptabil în condiţiile în care rezultatele trebuie obţinute cât mai repede, atât pentru a facilita diagnosticarea rapidă a pacientului, cât şi pentru a se putea rula simulări pentru mai mulţi pacienţi. Prin urmare, pentru a se reduce timpul de execuţie calea cea mai utilizată este procesarea paralelă a modelelor multiscalare. În literatura de specialitate s-au prezentat implementări ale unor astfel de modele atât pe clustere cât şi pe procesoare grafice.

În cadrul tezei de doctorat atenţia este îndreptată spre o clasă specială a modelelor multiscalare, şi anume modelele multiscalare de ordin redus, care sunt compuse din modele unidimensionale şi modele cu parametrii distribuiţi. Acestea permit estimarea corectă a debitelor şi presiunilor arteriale pentru geometrii specifice pacienţilor la o fracţiune din timpul necesar modelelor multiscalare care înglobează modele tridimensionale (diferenţa este de cel puţin două ordine de mărime). Ele nu pot estima mărimi locale precum tensiunea intraparietală, dar majoritatea indicatorilor de diagnoză, utilizaţi în spitale, sunt formulaţi doar pe baza presiunilor şi debitelor.

Modelele multiscalare de ordin redus dezvoltate în cadrul tezei au ca obiectiv principal, aplicarea lor într-un cadru clinic pentru a sprijini medicii în luarea deciziilor referitoare la modul de tratare a pacienţilor. În acest context, activităţile de cercetare desfăşurate în cadrul studiilor doctorate au fost îndreptate spre următoarele direcţii:


1. Personalizarea simulărilor hemodinamice: în cadrul simulărilor hemodinamice este important ca nu doar geometria, ci şi parametrii modelului multiscalar (condiţiile de frontieră, proprietăţile peretelui arterial) să fie personalizaţi. De multe ori nu sunt disponibile măsurători directe ale mărimilor de interes, caz în care acestea trebuie estimate pe baza altor mărimi măsurate şi a caracteristicilor fiziologice ale circuitului cardiovascular.

2. Tehnici de optimizare a simulărilor hemodinamice specifice pacienţilor: pentru ca rezultatele simulărilor hemodinamice să fie cât mai apropiate de caracteristicile hemodinamice reale ale pacientului este necesar ca mărimile simulate să corespundă celor ale pacientului (presiunea arterială medie, sistolică, diastolică, debit mediu, maxim, etc.). Pentru a se obţine această corespondenţă în literatura de specialitate s-au introdus tehnici de optimizare, care însă conduc la o creştere a timpului total de execuţie datorită necesităţii rulării repetitive a unei simulări pentru acelaşi pacient. În cadrul tezei se introduc tehnici de optimizare avansate care permit reducerea timpului de execuţie

3. Procesarea paralelă pe procesoare grafice pentru reducerea timpului de execuţie: o altă direcţie care trebuie urmărită pentru reducerea timpului de execuţie este procesarea paralelă a algoritmilor numerici utilizaţi pentru rezolvarea modelelor multiscalare. Deşi modelele multiscalare de ordin redus sunt deja mult mai rapide decât modelele multiscalare de ordin complet, timpul lor de execuţie poate fi în continuare redus printr-o astfel de abordare.

3. MODELAREA MULTISCALARĂ DE ORDIN REDUS A HEMODINAMICII ARTERIALE

3.1 MODELUL CVASI UNIDIMENSIONAL

Modelul cvasi unidimensional se deduce din ecuaţiile Navier-Stokes tridimensionale pe baza ipotezelor simplificatoare prezentate în capitolul anterior. Ecuaţiile fundamentale ale acestui model reprezintă legile de conservare a masei şi de conservare a impulsului:

0),(),(

x

txq

t

txA, (3.1)

),(

),(),(

ρ

),(

),(

),(α

),( 2

txA

txqK

x

txptxA

txA

txq

xt

txqR

, (3.2)

unde, x reprezintă poziţia iar t reprezintă timpul. A(x,t) este aria secţiunii transversale, p(x,t) este presiunea şi q(x,t) este debitul. Coeficienţii α şi KR reprezintă corecţia impuls-flux şi respectiv pierderile vâscoase datorită frecării cu peretele arterial. Parametrul KR este calculat cu ajutorul expresiei:

dy

d2

sK R , (3.3)

şi depinde de tipul profilului ales. Pentru un profil de tipul γ1 12γγ yys se obţine

. În cazul unui profil parabolic, 2γπν2 RK 8RK şi 3/4 , fiind vâscozitatea cinematică a fluidului.

Ecuaţiile (3.1) şi (3.2) nu pot fi rezolvate analitic, ci doar numeric. Majoritatea schemelor numerice utilizate pentru rezolvarea modelului unidimensional necesită scrierea ecuaţiilor sub formă conservativă. În vederea reformulării ecuaţiilor se introduce mărimea B:

txp

ppxrAtxpxrB

,

000

dp,1

,,

, (3.4)

dxrxρx0

0

drBpAB

. (3.5)

8


Dependinţele spaţiale şi temporale ale mărimilor au fost omise pentru a o obţine o reprezentare mai concisă a modelului. Deoarece termenul

9

dxr 00 drB nu conţine derivate

parţiale ale mărimii p şi prin urmare ale mărimilor A şi q, acesta poate fi calculat direct şi este însumat în ambii membrii ai ecuaţiei (3.2):

dxr),(

),(α 0

0

2 drB

txA

txqKB

A

q

xt

qR

. (3.6)

Prin urmare, sistemul de ecuaţii al modelului unidimensional poate fi scris sub formă matriceală, pe baza ecuaţiilor (3.1) şi (3.6), astfel:

dxr),(

),(0

α0

0

2 drB

txA

txqKB

A

qq

xq

A

t R. (3.7)

Sistemul de ecuaţii (3.7) conţine trei necunoscute (A, q şi p) şi doar două ecuaţii şi, prin urmare, este necesară o a treia ecuaţie pentru a se putea determina mărimile necunoscute. Această ecuaţie este ecuaţia de stare, care introduce o relaţie între presiunea din interiorul vasului şi secţiunea transversală (ecuaţia de stare descrie complianţa vaselor). Atunci când peretele arterial este modelat ca material pur elastic se foloseşte următoarea expresie:

00

00 ),(

)(1)(

3

4)(Ψ),( p

txA

xAx

r

EhpAtxp el

, (3.8)

unde, E este modulul lui Young, h este grosimea peretelui arterial, r0 este raza iniţială corespunzătoare presiunii iniţiale p0 şi A0 este aria iniţială a secţiunii transversale. Proprietăţile pereţilor sunt, în continuare, estimate pe baza unei legi empirice determinate prin intermediul rezultatelor experimentale [Olufsen et al., 2000]:

30210

kxrkexpkxr

Eh , (3.9)

unde, k1, k2 şi k3 sunt trei constante determinate experimental. Alternativ, se poate folosi şi o lege vâscoelastică. Pentru a modela vâscoelasticitatea,

peretele arterial este considerat a fi un material de tip Voigt [Fung, 1993], pentru care tensiunea depinde atât de deformare cât şi de derivata în timp a deformării. În literatura de specialitate s-au propus următoarele două ecuaţii de stare: modelul V1 [Malossi et al., 2012]:

00

00 1

3

4)()( p

t

A

AAA

A

r

EhpAAp S

vel

; (3.10)

modelul V2 [Alastruey et al., 2011]:

0

0

0

00 1

3

4)()( p

t

A

AAA

A

r

EhpAAp S

vel

. (3.11)

Pentru a se realiza o comparaţie adecvată între cele două modele vâscoelastice, coeficientul vâscos γS este considerat egal pentru ambele modele:

214

tg

hET SS

S , (3.12)

unde, TS este timpul caracteristic de undă, ΦS este unghiul viscoelastic, iar σ raportul Poisson. La fel ca şi în cazul modelului elastic, coeficientul vâscoelastic este considerat a fi variabil în spaţiu


( ). Prezenţa termenului vâscoelastic în ecuaţiile (3.11) / (3.12) introduce un termen

suplimentar în ecuaţia de conservare a debitului:

10

)(xSS

dxr),(

),(Ψ

ρα 0

0

2 drB

txA

txqK

x

AB

A

q

xt

qR

v

. (3.13)

3.2 CONDIŢII DE FRONTIERĂ

În funcţie de disponibilitatea măsurătorilor in-vivo şi de ipotezele utilizate în modelare, cercetătorii folosesc în mod tipic una din următoarele abordări la frontiera de intrare: (i) debit variabil în timp, (ii) model simplificat al inimii cuplat [Formaggia et al., 2006] sau, (iii) condiţie de frontieră fără reflexii precum o undă de presiune sau de debit [Mynard et al., 2012]. Un profil de viteze (sau de debit) variabil în timp poate fi determinat relativ simplu într-un cadru clinic, fiind, de obicei, determinat în mod uzual ca parte a fluxului de lucru de diagnosticare a pacientului (2D/3D PC-MRI, ultrasunete Doppler).

Condiţiile de frontieră de ieşire sunt în mod tipic reprezentate de modele cu parametrii distribuiţi în cadrul simulărilor hemodinamice.

3.2.1 Modelul Windkessel cu trei elemente

Acest model, introdus în capitolul anterior, va fi aplicat ca şi condiţie de frontieră aperiodică [Westerhof et al., 1971] şi este reprezentat de următoarea ecuaţie:

CR

RRq

CR

p

t

qR

t

p

d

dp

dp

)(

, (3.29)

unde, Rp este rezistenţa proximă, Rd este rezistenţa distantă şi C este complianţa, în timp ce p şi q se referă la presiunea şi respectiv debitul de la intrarea în modelul Windkessel.

3.3 REZOLVAREA NUMERICĂ A MODELULUI HEMODINAMIC UNIDIMENSIONAL

În cazul utilizării unui model elastic pentru peretele arterial, se obţine sistemul de ecuaţii (3.18), care reprezintă un sistem de ecuaţii hiperbolice. În cazul utilizării unui model vâscoelastic pentru peretele arterial, natura hiperbolică a ecuaţiilor este pierdută datorită termenului suplimentar din relaţia arie-presiune. Abordările pentru rezolvarea numerică a ecuaţiilor modelului unidimensional pot fi împărţite în două mari categorii:

1. Metode care nu exploatează natura hiperbolică (originală) a ecuaţiilor: metoda Galerkin discontinuă a elementelor finite cu termeni de stabilizare [Raghu et al., 2011], metode implicite cu diferenţe finite sau cu elemente spectrale în interiorul cărora termenii neliniare sunt evaluaţi iterativ la fiecare moment de timp discret prin intermediul metodei Newton [Reymond et al., 2011], [Bessems et al., 2008], etc.; 2. Metode care exploatează natura hiperbolică a ecuaţiilor în cazul în care se foloseşte un model elastic pentru peretele arterial [Olufsen et al., 2000], [Mynard et al., 2008], sau care, pentru a reveni la o formă hiperbolică, reformulează ecuaţiile în cazul utilizării unui model vâscoelastic pentru peretele arterial, prin aplicarea unei metode de împărţire a operatorilor pentru ecuaţia de conservare a impulsului. Această abordare a fost iniţial propusă în lucrarea [Formaggia et al., 2003] pentru un singur vas şi apoi reutilizată în lucrările [Passerini, 2009(a)] şi [Alastruey et al., 2011]. Metodele explicite se bazează fie pe metode de ordinul unu (metoda caracteristicilor), fie

pe metode de ordinul doi (metoda Lax-Wendroff în doi paşi [Olufsen et al., 2000], sau metoda descompunerii în serie Taylor [Mynard et al., 2008], [Mynard et al., 2012]). Datorită preciziei mai mari, metodele de ordinul doi sunt preferate în general şi au fost folosite şi în cadrul acestui


studiu. Metoda caracteristicilor este, în schimb, folosită la frontierele de intrare, de ieşire şi la bifurcaţiile arteriale.

3.3.1 Rezolvarea numerică a modelului unidimensional cu pereţi arteriali elastici

Pentru rezolvarea numerică a modelului unidimensional se foloseşte sistemul de ecuaţii (3.27), în care mărimile fizice sunt adimensionalizate. Pentru simplitatea notaţiilor, semnul tilda va fi în continuare omis. În vederea prezentării algoritmilor numerici, ecuaţia (3.27) este scrisă în forma matriceală:

SRU

xt, (3.44)

,dxdr

π2drdr

π2π8

0

,α

0

000

0002

1

00

2

2

1

drdfAfr

dfA

dfAfA

A

q

ReS

S

AAAfA

qq

R

R

q

A

S

RU

(3.45)

unde, U este vectorul mărimilor necunoscute, R este termenul de flux şi S este partea dreaptă. Pentru implementarea numerică a sistemului de ecuaţii (3.44) se consideră două metode de ordinul doi: metoda Lax-Wendroff şi metoda descompunerii în serie Taylor.

3.3.1.1 Metoda Lax-Wendroff

Prin

se notează soluţia numerică la momentul de timp discret , în

nodul m al grid-ului segmentului vascular curent. Metoda Lax-Wendroff (LW) constă din doi paşi. La pasul întâi se calculează mărimile la nodurile intermediare ale grid-ului; aceste mărimi sunt calculate la locaţii situate între nodurile grid-ului, şi prin urmare există doar valori interioare segmentului vascular:

tnxmnm Δ,ΔUU tnΔ

2222/12/12/12/12/12/12/1

nj

nj

nj

nj

nj

njn

j x

t SSRRUUU , (3.46)

unde, este pasul de timp discret iar t 21 /mj şi m fac referire la locaţiile pe grid.

La cel de-al doilea pas se calculează mărimile în nodurile grid-ului. Acest pas utilizează atât valorile de la pasul anterior cât şi de la nodurile intermediare:

2/12/1

2/12/1

2/12/1

2/12/1

2/12/1

2/12/1

1

2

Δ

Δ

Δ

nm

nm

nm

nm

nm

nm

nm

nm

t

x

tSSSSRRUU . (3.47)

unde, Δx este distanţa dintre două puncte ale grid-ului.

3.3.2 Rezolvarea numerică a modelului unidimensional cu pereţi arteriali vâscoelastici

În cazul utilizării unui model vâscoelastic pentru peretele arterial, se aplică o schemă de împărţire a operatorilor pentru ecuaţia de conservare a impulsului, pentru a reveni la o formă hiperbolică a ecuaţiilor. Prin urmare ecuaţia (3.6), este reformulată astfel sub formă adimensională (folosind şi notaţiile din secţiunea anterioară):

22 γ S

x

q

xA

x

R

t

q

. (3.59)

11


Această ecuaţie nu este hiperbolică şi nu poate fi reformulată sub formă conservatoare. Schema de împărţire a operatorilor are la bază ipoteza conform căreia influenţa termenului vâscoelastic este mică în comparaţie cu influenţa termenului elastic al modelului peretelui arterial. Debitul este considerat a fi compus dintr-o componentă elastică şi una vâscoelastică ( ) iar ecuaţia (3.59) este împărţită, la rândul ei, în două ecuaţii:

12

ve qqq

22 S

x

R

t

qe

. (3.60)

0

x

q

xA

t

qv . (3.61)

Prin urmare, determinarea soluţiei numerice, la fiecare moment de timp discret, este compusă din doi paşi secvenţiali. La primul pas sistemul de ecuaţii compus din ecuaţia (3.18a) şi (3.60) este rezolvat ca şi în cazul utilizării unui model elastic pentru peretele arterial, determinându-se astfel mărimile A(x,t) şi qe(x,t). La cel de-al doilea pas, ecuaţia (3.61) este rezolvată pentru a se determina mărimea qv(x,t) şi astfel şi debitul total q(x,t).

3.4 ALGORITMUL DE CUPLARE

Rezultate contradictorii au fost raportate în literatura de specialitate în ceea ce priveşte efectul termenului vâscoelastic al modelului peretelui arterial în comparaţie cu utilizarea unei legi pur elastice: undă de presiune similară dar o variaţie mai mică a ariei transversale [Raghu et al., 2011], respectiv o presiune de puls mai mare dar variaţii similare ale ariei transversale [Passerini, 2009(a)].

Lucrările [Alastruey et al., 2011] şi [Passerini, 2009(a)] au neglijat componenta vâscoasă a modelului peretelui arterial la punctele de cuplare a segmentelor vasculare (frontiere de intrare şi ieşire, joncţiuni), şi acest aspect a fost identificat ca posibilă cauză pentru concluziile contradictorii menţionate mai sus. Pentru a se verifica dacă această ipoteză simplificatoare nu influenţează rezultatele, s-a dezvoltat un algoritm de cuplare nou (v. alg. 3.1), care poate fi aplicat în egală măsură la frontierele de intrare, de ieşire cât şi la bifurcaţii [Itu et al., 2013(b)]. Algoritmul se bazează pe metodologia introdusă anterior în lucrarea [Passerini et al., 2009(b)] pentru cuplarea modelelor tridimensionale şi unidimensionale.

Pentru domeniul din amonte se consideră o condiţie de frontieră de presiune, în timp ce, pentru domeniile din aval se aplică o condiţie de frontieră de debit. În continuare, exponentul n face referire la soluţia la momentul tn, în timp ce exponentul n+1 face referire la soluţia la momentul tn+1. Indicele k face referire la iteraţiile realizate la fiecare moment de timp pentru determinarea mărimilor necunoscute ale celor m segmente.

În cadrul acestor ecuaţii, debitul vasului din amonte a fost considerat pozitiv în direcţia de ieşire din vas, în timp ce, debitul vaselor din aval a fost considerat pozitiv în direcţia de intrare în vas (χ este un parametru de relaxare).

Algoritmul 0.1 Cuplarea a m domenii.

(1) Iniţializare: , . 0k miAAppqq ni

ni

ni

ni

ni

ni ,...,1,,, 1

0,1

0,1

0,

(2) Iterare (k)

(2.1) Calculul debitului total, care trebuie distribuit vaselor din

aval:

, (3.63)

m

i

nki

nk

nkm qqq

2

1,

1,1

11,2 χ1χ

(2.2) Rezolvarea modelelor din aval având la bază următorul set de

condiţii:


13

, (3.64) 11,2

2

11,

n

km

m

i

nki qq

. (3.65) 1,...,2,11,1

11,

mipp n

kin

ki

Ca urmare, se determină şi , pentru 11,

nkip 1

1,

nkiq mi ,...,2 .

(2.3) Rezolvarea modelului din amonte folosind condiţia de frontieră:

→ . (3.66) 11,2

11,1

n

kn

k pp 11,1

nkq

(2.4) 1 kk

(3) Testarea condiţiei de convergenţă: dacă

m

i

nki

nk qq

2

11,

11,1 , schema

de calcul trece la următorul moment de timp (pasul 1), altfel

algoritmul iterativ de cuplare trece la următoarea iteraţie

(pasul 2).

3.6 REZULTATE

Rezultatele prezentate în această secţiune se concentrează asupra a două aspecte. În primul rând va fi prezentată o comparaţie a rezultatelor obţinute cu cele două tipuri de condiţii de frontieră prezentate în secţiunea 3.2 (modelul Windkessel cu 3 elemente şi arborele structurat). În al doilea rând se va analiza influenţa utilizării algoritmului de cuplare introdus în secţiunea anterioară asupra rezultatelor obţinute cu modele vâscoelastice ale peretelui arterial.

Sângele a fost modelat ca fluid Newtonian cu o densitate de 1.055 g/cm3 şi o vâscozitate dinamică de 0.045 dynes/cm2s pentru toate simulările prezentate în această secţiune. Rezultatele sunt raportate pentru arborele arterial detaliat prezentat în lucrarea [Stergiopulos et al., 1992] şi reprezentat în figura 3.7. Acesta este compus din 51 artere. La frontiera de intrare s-a impus un debit variabil în timp [Olufsen et al., 2000].

3.6.1 Comparaţie între rezultatele obţinute cu diferite tipuri de condiţii de frontieră de ieşire

Figura 3.8 prezintă presiunea, debitul şi aria transversală variabile în timp la cele cinci locaţii marcate cu un cerc albastru în figura 3.7. Fiecare figură conţine două grafice, care au fost obţinute cu condiţiile de frontieră de tip Windkessel (WK), respectiv de tip arbore structurat (ST). În ambele cazuri s-a folosit un model elastic al peretelui arterial. Deoarece rezistenţa totală introdusă de oricare dintre cele două condiţii de frontieră este similară, valorile medii ale celor trei mărimi reprezentate sunt aproximativ egale la orice locaţie din interiorul arborelui arterial.

Astfel, presiunile obţinute cu condiţia de frontieră ST scad mai târziu în cadrul unui ciclu cardiac, indicând faptul că undele reflectate sosesc la locaţia respectivă mai târziu. Acest aspect este mai pronunţat pentru arterele proxime ale arborelui arterial şi poate fi explicat astfel: condiţia de frontieră de tip ST simulează propagarea undelor până la nivelul arteriolelor unde au loc principalele reflexii, în timp ce condiţia de frontieră WK, fiind un model cu parametrii distribuiţi, nu este capabilă să modeleze propagarea undelor dincolo de frontierele de ieşire, introducând astfel reflexii doar la frontierele de ieşire. Ca urmare a sosirii întârziate a undelor reflectate, şi presiunea maximă aste atinsă mai târziu. Aceste aspecte conduc, de asemenea, la oscilaţii mai pronunţate în cadrul formelor de undă ale debitelor afişate în cea de-a doua coloană a figurii 3.8. În ceea ce priveşte ariile transversale, în general, variaţia în cadrul unui ciclu cardiac este mai mare în cazul utilizării unei condiţii de frontieră de tip ST. Deoarece s-a folosit un model elastic, presiunea şi aria transversală sunt în fază, iar variaţia mai pronunţată a ariei transversale este reflectată de o variaţie mai mare a presiunilor. Presiunea de puls mai mare obţinută pentru condiţia de frontieră de tip ST indică o complianţă mai scăzută decât cea impusă


prin condiţia de frontieră de tip WK. Se subliniază faptul că, pentru arterele modelate cvasi unidimensional complianţa este independentă de condiţia de frontieră utilizată.

Fig. 3.8 Variaţia în timp a presiunii, debitului şi ariei transversale la rădăcina aortei (a),

aorta descendentă (b), aorta abdominală (c), artera femurală (d) şi artera subclavie (e) (corespunzător locaţiilor A-E din figura 3.7), pentru condiţii de frontieră de tip Windkessel şi arbore structurat.

14


Suplimentar, în figura 3.9 sunt prezentate variaţiile în timp şi spaţiu de-a lungul segmentului 27 al arborelui arterial (aorta descendentă). Se poate observa propagarea undelor atât în timp cât şi în spaţiu.

Fig. 3.9 Variaţiile în timp şi spaţiu ale presiunii (a) şi debitului în aorta descendentă (b).

3.6.2 Influenţa algoritmului de cuplare asupra rezultatelor obţinute cu modele vâscoelastice ale pereţilor arteriali

Parametrii coeficientului vâscoelastic γS au fost setaţi după cum urmează: TS = 0.30 s, ΦS = 10°, σ = 0.5. Toate punctele de cuplare au fost rezolvate prin intermediul algoritmului de cuplare introdus în secţiunea anterioară. Parametrul de relaxare χ a fost setat la valoarea 0.5.

S-au folosit în total cinci configuraţii de calcul: (a) model elastic; (b) model V1 aplicat la punctele interioare şi de cuplare; (c) model V1 aplicat doar la punctele interioare (V1-int); (d) model V2 aplicat la punctele interioare şi de cuplare; (e) model V2 aplicat doar la punctele interioare (V2-int). Acestea au fost testate prin intermediul a două configuraţii arteriale: (i) o singură arteră; (ii) arbore arterial complet.

Suplimentar, s-au comparat rezultatele obţinute prin intermediul algoritmului clasic de cuplare, bazat pe metoda caracteristicilor [Olufsen et al., 2000] cu rezultatele obţinute prin intermediul noului algoritm de cuplare, prezentat în secţiunea anterioară (cu un model elastic al peretelui arterial). S-a calculat norma L2 a diferenţelor absolute dintre cele două soluţii numerice, rezultatele fiind de ordinul 10-6 cm2 pentru aria transversală, 10-5 ml/s pentru debit şi 10-5 mmHg pentru presiune, demonstrându-se astfel că ambii algoritmi de cuplare conduc practic la aceleaşi rezultate.

3.6.2.2 Arbore arterial integral

În continuare sunt prezentate rezultatele pentru arborele arterial complet din figura 3.7. În figura 3.11 sunt prezentate variaţiile în timp ale presiunilor şi debitelor la locaţiile A, B, D şi respectiv E ale arborelui din figura 3.7. Se poate astfel observa că modelarea efectului vâscoelastic conduce la reducerea oscilaţiilor de înaltă frecvenţă ale formelor de undă de debit şi presiune, fenomenul acesta fiind mai pronunţat la locaţiile distante ale arborelui arterial. Aceste observaţii corespund rezultatelor raportate în literatură [Reymond et al., 2011]. Suplimentar, aria transversală şi presiunea nu mai sunt în fază, valoarea maximă a ariei transversale fiind, în general, atinsă la un moment de timp ulterior (acest aspect poate fi observat în figura 3.12(a)). O variaţie mai mică a ariei transversale a fost obţinută pentru modelele V1 şi V2, în timp ce, pentru modelele V1-int şi V2-int, variaţia este similară celei obţinute pentru o lege elastică. Suplimentar, în figura 3.12(b) se prezintă bucla de histerezis obţinută pentru relaţia presiune-arie, tot în aorta descendentă (în cazul unui model elastic se obţine o variaţie liniară arie-presiune). Rezultatele sunt similare pentru modelele V1 şi V2. Aria buclei de histerezis este proporţională cu energia disipată.

15


16

Fig. 3.11 Variaţiile în timp ale presiunii şi debitului la rădăcina aortei, în aorta

descendentă, în artera femurală şi respectiv în artera subclavie stângă.

Fig. 3.12 Aria transversală variabilă în timp în aorta descendentă (a); variaţia

presiune-arie în aorta descendentă (b).


Pentru a cuantifica diferenţa dintre rezultatele obţinute cu modelele V1 şi V1-int, respectiv cu modelele V2 şi V2-int, s-au calculat diferenţele relative maxime şi medii la locaţiile A, B, D şi E din figura 3.7. Pentru a evalua corect aceste rezultate, s-au calculat, de asemenea, diferenţele relative maxime şi medii dintre rezultatele obţinute cu modelele V1 şi V2 şi respectiv cu modelul elastic. Rezultatele sunt prezentate în tabelul 3.6. Deşi diferenţele dintre modelele V1 şi V1-int, respectiv dintre modelele V2 şi V2-int par mici, ele sunt comparabile cu diferenţele dintre modelele vâscoelastice şi modelul elastic. Pentru aorta descendentă de exemplu, diferenţele sunt chiar de acelaşi ordin de mărime. Chiar şi la locaţiile distante, unde influenţa componentei de vâscozitate este mai mare, neglijarea acesteia la punctele de cuplare introduce erori care reprezintă 1/5 sau mai mult din diferenţa dintre modelele vâscoelastice şi modelele elastice.

Tab. 3.6 Diferenţele relative medii şi maxime obţinute între diferite implementări ale modelelor vâscoelastice şi între modelele vâscoelastice şi elastice la locaţiile A, B, D, E identificate în figura 3.7.

V1 – V1-int V1 – elastic V2 – V2-int V2 – elastic Arteră Diferenţă

P [%] Q [%] P [%] Q [%] P [%] Q [%] P [%] Q [%]

Rel. Med. 0.545 0.0 0.954 0.0 0.620 0.0 1.088 0.0 Aorta asc. Rel. Max. 1.227 0.0 3.896 0.0 1.397 0.0 4.353 0.0

Rel. Med. 0.499 0.851 0.897 1.343 0.568 0.956 1.019 1.472 Aorta desc. Rel. Max. 0.961 2.970 3.955 3.447 1.100 3.289 4.410 3.887

Rel. Med.. 0.310 0.984 1.528 3.866 0.342 1.079 1.697 4.278 Artera femur. Rel. Max. 1.104 3.030 7.812 11.66 1.241 3.365 8.484 12.96

Rel. Med.. 0.522 1.019 1.128 3.982 0.594 1.090 1.287 4.387 Artera subcl. Rel. Max. 1.231 3.417 5.774 17.74 1.390 3.652 6.365 19.12

3.7 CONCLUZII

În acest capitol s-a introdus un model multiscalar de ordin redus pentru simularea hemodinamicii arteriale. Iniţial s-a introdus modelul cvasi unidimensional, care reprezintă componenta de bază a modelului multiscalar. Mărimile fundamentale ale acestui model sunt aria transversală, debitul şi presiunea. Dat fiind faptul că aria transversală poate varia în timp, modelul unidimensional permite surprinderea principalelor aspecte ale interacţiunii dintre sânge şi peretele arterial. Astfel, pentru modelarea peretelui arterial s-au considerat o lege elastică şi două legi vâscoelastice.

4. PROCEDURĂ DE OPTIMIZARE EVOLUATĂ PENTRU PERSONALIZAREA SIMULĂRILOR HEMODINAMICE

4.2 METODOLOGIA PROCEDURII DE OPTIMIZARE EVOLUATE

Se propune o procedură de optimizare evoluată pentru obţinerea unor caracteristici hemodinamice dorite în cadrul modelelor multiscalare, prezentată în figura 4.1. În continuare, se utilizează un model multiscalar de ordin redus, care are la bază modelul unidimensional prezentat pe larg în capitolul 3.

La frontiera de intrare s-au folosit profiluri de debit variabil în timp iar la frontierele de ieşire s-au impus modele Windkessel cu trei elemente (se reia ec. (3.29)):

CR

)RR(q

CR

pqR

p

d

dp

dp

tt, (4.1)

17


unde, p reprezintă presiunea la frontiera de ieşire, q este debitul şi Rp şi Rd sunt rezistenţele proxime respectiv distante, iar C este complianţa. Modelele Windkessel au fost impuse ca şi condiţii de frontieră non-periodice.

Caracteristicile undelor de presiune şi de debit au fost folosite ca obiective ale procedurii de optimizare. Parametrii adaptaţi în vederea obţinerii anumitor caracteristici hemodinamice au fost, fie parametrii modelelor Windkessel, fie mărimi formulate pe baza acestora. Astfel, s-au adaptat rezistenţa totală a modelului Windkessel, Rt, raportul dintre rezistenţa proximă şi rezistenţa distantă, ρ, şi constanta de timp a scăderii exponenţiale a presiunii atunci când debitul de intrare în modelul Windkessel este zero, τ:

dpt RRR , (4.2)

tp R/R , (4.3)

CRd . (4.4)

Determinarea valorilor acestor parametri a fost formulată ca soluţie a unui sistem de ecuaţii neliniare cu o rădăcină în punctul în care mărimile simulate şi cele de referinţă au aceeaşi valoare. Pentru a se determina valorile reziduurilor (ale funcţiilor obiectiv: f(xi)) pentru un set de valori ale parametrilor xi, este necesară rularea unei simulări hemodinamice cu valorile xi.

După cum este prezentat în figura 4.1(a), sistemul de ecuaţii neliniare este iniţial rezolvat pentru un model 0D, compus din modelele Windkessel impuse ca şi condiţii de frontieră de ieşire în cadrul modelului multiscalar. Pentru a se determina o soluţie iniţială pentru modelul 0D, s-a considerat un grid de valori fiziologice ale parametrilor. Setul de valori care conduce la cea mai mică normă L2 pentru reziduurile obiectivelor (normalizată prin intermediul valorilor tipice alese), a fost ales ca soluţie iniţială pentru metoda dogleg trust region [Nocedal et al., 2006], prin intermediul căreia s-a determinat soluţia x0 utilizată în cadrul următorilor paşi. Figura 4.1(b) prezintă modelul 0D corespunzător modelului multiscalar al unei geometrii arteriale a aortei proxime (în cadrul modelului 0D s-a folosit aceeaşi condiţie de frontieră ca şi în cadrul modelului multiscalar).

În continuare, s-a folosit o abordare de tip punct fix pentru a se determina o matrice Jacobiană calculată cu diferenţe finite pe baza valorilor variaţiilor tipice ale parametrilor şi care să fie consistentă cu valorile tipice ale reziduurilor. Componentele Jacobian-ului aproximat au fost calculate cu ajutorul ecuaţiei:

ijtipjj

tipjtip

jij ss

sJ eexfexf

2

1

2

1100 , (4.5)

unde, ei şi ej reprezintă vectorii unitate în direcţiile i şi j iar este variaţia tipică a parametrului

j calculată prin intermediul ecuaţiei:

tipjs

eqn

i

typiij

typj fJs

1

1 , (4.6)

unde, este valoarea tipică a reziduului obiectivului fi. typ

if

După cum este prezentat în figura 4.1(a), procedura compusă din ecuaţiile (4.5) şi (4.6) este o procedură iterativă care este finalizată când norma euclidiană a diferenţei dintre două matrice Jacobiane consecutive, normalizate cu valorile corespunzătoare ale şi , este mai

mică decât 10-6. Valorile variaţiilor tipice ale parametrilor, împreună cu valorile tipice ale reziduurilor obiectivelor sunt folosite pentru normalizarea mărimilor.

typif typ

js

În continuare, se creează modelul multiscalar şi sunt introduşi doi paşi suplimentari (evidenţiaţi prin îngroşarea textului), înainte de a rula simulările cu modelul multiscalar. Modelul multiscalar corespunzător geometriei aortei proxime cu coarctaţie, specifice pacientului, este

18


prezentat în figura 4.1(c) (obiectivele acestei configuraţii au fost presiunile sistolice şi diastolice şi raporturile de divizare a debitului mediu pentru fiecare arteră supraortică şi pentru aorta descendentă).

Fig. 4.1 Fluxul de lucru al procedurii de optimizare evoluate (a); model 0D utilizat în cadrul

primilor paşi ai procedurii de optimizare (b); model multiscalar utilizat pentru determinarea valorilor finale ale parametrilor în cadrul simulării hemodinamice realizate pentru geometria aortei proxime cu coarctaţie (c).

19


Ideea de bază a procedurii de optimizare evoluate este de a compensa proprietăţile hemodinamice ale modelului cu ordinul cel mai mare din interiorul modelului multiscalar, atunci când se trece de la modelul 0D la modelul multiscalar. În particular, sunt considerate două proprietăţi ale modelului multiscalar: complianţa şi rezistenţa. Pentru cazul din figurile 4.1(b) şi 4.1(c), este crucial să se compenseze complianţa (aceasta are valorile cele mai mari în aorta proximă) şi rezistenţa coarctaţiei. Pe de altă parte, compensarea rezistenţei segmentelor sănătoase nu este crucială deoarece arterele mari au cea mai mică rezistenţă. Compensarea rezistenţei segmentelor sănătoase devine importantă atunci când zona de interes este reprezentată de artere mici şi/sau zona de interes conţine un număr mare de segmente arteriale.

La primul pas suplimentar se estimează complianţa şi rezistenţa modelului cu ordinul cel mai mari (acestea nu pot fi determinate exact apriori). Aceste două mărimi nu au fost luate în considerare la rezolvarea sistemului de ecuaţii neliniare pentru modelul 0D. Prin urmare, la cel de al doilea pas suplimentar, valorile parametrilor modelelor Windkessel utilizate la frontierele de ieşire sunt recalculate pentru a compensa proprietăţile hemodinamice ale modelului cu ordinul cel mai mare. Se subliniază faptul că valorile parametrilor sistemului neliniar nu sunt modificate la acest pas ci se modifică doar parametrii modelelor Windkessel.

Odată ce valorile parametrilor modelelor Windkessel au fost adaptate corespunzător, se rulează o simulare cu modelul multiscalar şi se evaluează reziduurile obiectivelor. La prima rulare cu modelul multiscalar, variabilele modelului multiscalar sunt iniţializate pe baza rezultatelor obţinute cu modelul 0D. Fiecare simulare, realizată pentru un anumit set de valori ale parametrilor, este rulată până când normele L2 ale diferenţelor normalizate dintre profilurile de presiune şi de debit obţinute la ciclul cardiac curent şi la cel anterior devin mai mici de 10-5.

Dacă toate reziduurile obiectivelor sunt mai mici decât limita de toleranţă (considerată în acest caz egală cu ), procedura de optimizare este finalizată. Dacă, pe de altă parte,

obiectivele nu sunt îndeplinite, se utilizează o metodă cvasi-Newton pentru actualizarea valorilor parametrilor:

10typ /f i

iiii xfJxx 11

. (4.7)

Deoarece determinarea matricei Jacobiene prin intermediul unor formule cu diferenţe finite este prea costisitoare din punct de vedere al timpului de calcul, Jacobianul, determinat cu ajutorul abordării de tip punct fix, este actualizat la fiecare iteraţie astfel:

is

ti

t

isiiiiii sDs

sDsJxfxfJJ

2

21

1

, (4.8)

unde, este pasul curent iar Ds este o matrice diagonală de scalare: iii xxs 1

.,0

,,/1

ji

jisD

typj

ijs (4.9)

Obiectivele utilizate pentru diversele modele multiscalare din această lucrare sunt: presiunea maximă, Pmax, presiunea minimă, Pmin, şi raportul de divizare a debitului pentru vasele terminale, Φ. Valoarea tipică a tuturor obiectivelor formulate pe baza presiunilor a fost setată la 1 mmHg şi la 0.005 pentru raporturile de divizare a debitelor.

4.2.2 Estimarea şi compensarea complianţei

Pentru estimarea complianţei unui segment vascular i al modelului cu ordinul cel mai mare din cadrul modelului multiscalar, , se folosesc proprietăţile mecanice ale peretelui arterial: iHOC

iL

ii

iiiHO dl

hE

lrlrC

0

2

2

3, (4.14)

20


unde, Ei este modulului lui Young şi hi este grosimea peretelui. Alternativ, se poate folosi şi metoda introdusă în lucrarea [Ismail et al., 2013]. Majoritatea complianţei din interiorul unui arbore arterial/venos se află în vasele mari şi, prin urmare, complianţa modelului cu ordinul cel mai mare poate deveni semnificativă în comparaţie cu complianţa elementelor Windkessel, în special, atunci când se folosesc domenii cu artere mari. În continuare, se determină complianţa totală a modelului cu ordinul cel mai mare

21

i

iHOHO CC care este apoi distribuită la

frontierele de ieşire ale modelului geometric, pe baza raportului ariilor transversale:

n

jjjHOjaMS rrCC

1

22 , (4.15)

unde, reprezintă valoarea presupusă a complianţei suplimentare introdusă la fiecare

frontieră de ieşire a modelului cu ordinul cel mai mare din cadrul modelului multiscalar. Scopul este acela să se determine complianţa modelelor Windkessel folosite la frontierele de ieşire ale modelului multiscalar, , astfel încât să se obţină pentru modelul multiscalar aceeaşi

complianţă totală ca şi în cazul modelului 0D,

jaMSC

jMSC

jDC0 . Figura 4.2(a) prezintă reprezentarea

clasică a modelului multiscalar, compusă din modelul cu ordinul cel mai mare şi elementul Windkessel impus la frontiera de ieşire. În figura 4.2(b), complianţa modelului cu ordinul cel mai mare este concentrată într-o complianţă 0D la frontiera de ieşire j, jaMSC , plasată înaintea

modelului Windkessel. Pentru determinarea complianţelor elementelor Windkessel ale modelului multiscalar, , trebuie transferată în interiorul modelului Windkessel, după cum

este prezentat în figura 4.2(c). Pentru ca modelele din figurile 4.2(c) şi 4.2(d) să fie echivalente, trebuie să fie valabilă următoare expresie:

jMSC jaMSC

jWKaMSjDjMS CCC 0 . (4.16)

Prin urmare, pentru a determina jMSC , trebuie calculată jWKaMSC .

reprezintă valoarea presupusă a complianţei suplimentare introdusă de modelul multiscalar la frontiera de ieşire j şi transferată în interiorul modelului Windkessel.

jWKaMSC

Au fost testate trei abordări diferite pentru calculul mărimii jWKaMSC : o abordare

directă, o abordare analitică şi o abordare numerică.

Fig. 4.2 Reprezentarea clasică a modelului multiscalar (a); complianţa

modelului cu ordinul cel mai mare concentrată într-o complianţă 0D, plasată în faţa modelului Windkessel (b); complianţa modelului cu ordinul cel mai mare transferată în interiorul modelului Windkessel (c); reprezentarea clasică a modelului Windkessel (d).


Abordarea 1 – Compensarea directă a complianţei. Influenţa rezistenţei proxime a modelului Windkessel este neglijată:

22

jaMSjWKaMS CC . (4.17)

Abordarea 2 – Compensarea analitică a complianţei. Abordarea 1 nu ia în considerare faptul că rezistenţa proximă a elementului Windkessel diminuează influenţa complianţei. Prin urmare, folosind o metodă introdusă anterior [Grinberg et al., 2008], este

determinată după cum urmează:

jWKaMSC

jMStjMSpjaMSjWKaMS RRCC 1 , (4.18)

unde, jMSpR este rezistenţa proximă a modelului Windkessel de la frontiera j.

Abordarea 3 – Compensarea numerică a complianţei. Se propune o abordare numerică bazată pe metoda presiunii de puls (PPM) [Stergiopulos et al., 1994], care este descrisă în algoritmul 4.2. Metoda PPM a fost introdusă pentru estimarea complianţei din aval de o locaţie anume din interiorul arborelui arterial (debitul variabil în timp şi presiunea medie şi de puls la locaţia respectivă sunt cunoscute). În continuare, se foloseşte o metodă PPM modificată bazată pe modelul Windkessel cu 3 elemente (PPM-WK3), care înlocuieşte modelul Windkessel cu două elemente al metodei PPM originale.

Algoritmul 4.2 Estimarea bazată pe PPM a complianţei.

(1) for fiecare vas terminal j

(2) Rulare simulare cu model WK cu două elemente

jaMSjMStj C,R,tQ D0 → jPPref

(3) Iniţializare PPM-WK3 → jtQ D0 , jMSpR , jMSdR , jaMSC ,

jPPref

(4) Execuţie PPM-WK3 → jPPMaMSC

(5) ← jWKaMSC jPPMaMSC

(6) final (for)

4.3 REZULTATE

Procedura de optimizare evoluată a fost testată cu ajutorul a cinci configuraţii diferite. Cel de al cincilea proces de optimizare (v. fig. 4.3) foloseşte o geometrie arterială a unui pacient, formată din rădăcina aortei, trei artere supraortice (artera brahiocefalică care se bifurcă în artera subclavie dreaptă şi artera carotidă comună dreaptă, artera carotidă comună stângă şi artera subclavie stângă) şi aorta descendentă care prezintă o coarctaţie (v. fig. 4.3c).

Rezultate sunt prezentate în continuare pentru patru proceduri de optimizare diferite: procedura de optimizare iniţială (OTP) propusă în lucrarea [Spilker et al., 2010], procedura de optimizare evoluată cu compensare a rezistenţelor şi cu compensare directă a complianţei (ATPDCC), procedura de optimizare evoluată cu compensare a rezistenţelor şi cu compensare analitică a complianţei (ATPACC) şi procedura de optimizare evoluată cu compensare a rezistenţelor şi cu compensare numerică a complianţei (ATPNCC).

4.3.5 Geometrie aortică specifică unui pacient cu coarctaţie

Această geometrie, prezentată în figura 4.3(c), a fost folosită anterior în literatura de specialitate [***MICCAI, 2012], debitul variabil în timp impus la frontiera de intrare fiind cunoscut. Suplimentar faţă de presiunea sistolică şi diastolică, raporturile de divizare ale debitului mediu către fiecare arteră supraortică şi către aorta descendentă au fost folosite ca obiective ale procedurii de optimizare.


23

Fig. 4.3 Model idealizat al arterei carotide comune şi profilul debitului de intrare corespunzător (a); Model idealizat al bifurcaţiei iliace şi profilul debitului de intrare corespunzător (b); Model al aortei proxime cu coarctaţie (c).

Rezistenţa totală a fiecărei artere supraortice, rezistenţa totală a aortei descendente şi complianţa totală au fost parametrii adaptaţi. Prin urmare, sistemul de ecuaţii neliniare rezolvat prin intermediul procedurilor de optimizare a fost formulat astfel:

0

0

0

0

0

refDAosimDAo

refLCCsimLCC

refBCsimBC

refminsimmin

refmaxsimmax

DAot

LSt

LCCt

BCt

PP

PP

C

R

R

R

R

f , (4.25)

unde, indicele BC reprezintă artera brahiocefalică, indicele LCC face referire la artera carotidă comună stângă iar indicele DAo face referire la aorta descendentă. Doar trei dintre cele patru rapoarte de divizare ale debitelor au fost folosite ca obiective în ecuaţia (4.25), deoarece cel de al


patrulea raport este obţinut ca diferenţă. Valorile scalate ale reziduurilor obiectivelor precum şi ale parametrilor, obţinute cu cele patru proceduri de optimizare sunt prezentate în figura 4.9.

Fig. 4.9 Evoluţia procedurii de optimizare pentru modelul aortei proxime cu coarctaţie.

24


În timp ce, sunt necesare trei iteraţii cu procedura OTP, două iteraţii sunt necesare cu procedurile ATPDCC şi ATPACC şi doar o singură iteraţie este necesară cu procedura ATPNCC.

La fel ca şi în cazurile anterioare, când se foloseşte procedura OTP, presiunea de puls la iteraţia zero este semnificativ mai mică decât presiunea de puls de referinţă, datorită complianţei suplimentare a modelului cu ordinul cel mai mare. În special, presiunea diastolică iniţială este mult apropiată de valoarea finală a acesteia atunci când se folosesc versiunile ATP. Acest aspect este, de asemenea, reflectat de faptul că valoarea complianţei la iteraţia zero este considerabil mai aproape de valoarea finală atunci când se foloseşte una din procedurile ATP propuse.

Suplimentar, toate valorile iniţiale ale raporturilor de divizare ale debitului sunt mai apropiate de valorile lor finale atunci când se folosesc versiunile ATP. Acest aspect este cauzat, în principal, de estimarea iniţială superioară a rezistenţei totale a aortei descendente Rt-DAo. Estimarea iniţială a Rt-DAo în cazul versiunilor ATP este superioară, deoarece ecuaţia (4.12) este utilizată pentru calculul rezistenţei segmentului de coarctaţie. Pe de altă parte, rezistenţele segmentelor sănătoase din aorta proximă sunt cele mai mici din întreaga circulaţie arterială şi compensarea lor are o importanţă secundară, în comparaţie cu compensarea complianţei şi a rezistenţei segmentului patologic.

Din cadrul procedurilor ATP, versiunea ATPNCC conduce la cele mai bune performanţe. Figura 4.10 prezintă o comparaţie a variaţiilor în timp ale presiunii la frontiera de ieşire a arterei subclavie stângi, corespunzătoare soluţiei iniţiale obţinute cu procedura OTP, soluţiei iniţiale obţinute cu procedura ATPNCC precum şi soluţiei finale (identică pentru ambele proceduri de optimizare). Soluţia iniţială furnizată de procedura ATPNCC este evident superioară soluţiei iniţiale obţinute cu procedura OTP şi aproape identică cu soluţia finală.

Fig. 4.10 Comparaţie între profilurile de presiune obţinute la frontiera de intrare a modelului aortei proxime, corespunzătoare soluţiei iniţiale cu OTP, soluţiei iniţiale cu ATPNCC şi soluţiei finale.

4.4 ANALIZA REZULTATELOR

Rezultatele prezentate în secţiunea anterioară demonstrează avantajele procedurii de optimizare evoluate, aceasta conducând la rezultate superioare pentru toate cele cinci configuraţii. Soluţiile iniţiale obţinute pentru modelele multiscalare sunt îmbunătăţite pentru toate cele cinci configuraţii iar numărul iteraţiilor este redus de la două sau trei iteraţii la una singură. După cum a fost analizat în lucrarea [Spilker et al., 2010], o componentă crucială pentru adaptarea cu succes a parametrilor este utilizarea unui model 0D la începutul procedurii de optimizare. Dezavantajul modelului 0D, compus doar din modelele Windkessel aplicate la frontierele de ieşire ale modelului multiscalar, este acela că nu ia în considerare proprietăţile hemodinamice ale domeniului proxim. Cele mai importante două proprietăţi ale domeniului

25


proxim sunt complianţa şi rezistenţa. Pornind de la ideea că aceste două proprietăţi trebuie compensate, au fost dezvoltate trei versiuni diferite ale procedurii de optimizare evoluate. Procedura de estimare şi compensare a rezistenţelor modelului cu ordinul cel mai mare este identică pentru toate cele trei versiuni, care diferă doar din punctul de vedere al compensării complianţiei. Diferenţa dintre procedurile ATP şi procedura OTP este reprezentată doar de cei doi paşi suplimentari prezentaţi în figura 4.1(a). Prin urmare, parametrii sistemului neliniar nu mai reprezintă direct proprietăţi ale modelelor Windkessel utilizate la frontierele de ieşire ci proprietăţi generale ale întregului model multiscalar.

26

Din punct de vedere computaţional, procedurile ATPDCC şi ATPACC sunt mai uşor de implementat şi au un timp de execuţie mai mic deoarece sunt metode analitice. Procedura ATPNCC este o abordare numerică, bazată pe o variantă modificată a metodei presiunii de puls. În comparaţie cu timpii de calcul mari necesari pentru rularea simulărilor cu modelul multiscalar, timpul de execuţie suplimentar necesar pentru procedura ATPNCC, în comparaţie cu procedurile ATPDCC şi ATPACC, este neglijabil. Deoarece, pentru cele cinci configuraţii diferite considerate în secţiunea 4.3, procedura ATPNCC conduce la cele mai bune performanţe, această variantă este considerată a fi cea mai bună dintre cele trei propuse. Un alt avantaj al procedurii de optimizare evoluate este faptul că, atât configuraţia modelului 0D cât şi configuraţia modelului multiscalar, nu sunt afectate de paşii suplimentari, arhitectura software dezvoltată pentru procedura OTP putând fi complet reutilizată.

Motivaţia dezvoltării procedurii de optimizare avansate nu a fost doar să se reducă numărul de iteraţii şi, prin urmare, şi timpul total de execuţie, dar şi să se obţină o soluţie iniţială mai apropiată de soluţia finală.

4.5 CONCLUZII

În acest capitol s-a introdus o procedură de optimizare evoluată în vederea obţinerii unor caracteristici hemodinamice dorite în cadrul simulărilor multiscalare realizate pentru geometrii arteriale sănătoase şi patologice.

5. MODELAREA MULTISCALARĂ DE ORDIN REDUS A CIRCULAŢIEI CORONARIENE ŞI DIAGNOSTICAREA NON-INVAZIVĂ A STENOZELOR CORONARIENE

5.1 INTRODUCERE

Circulaţia coronariană are o serie de caracteristici specifice care o deosebesc de restul componentelor sistemului cardiovascular. Aceasta este, pe de o parte, responsabilă de generarea presiunii arteriale necesară pentru alimentarea cu sânge a arborelui arterial sistemic dar, în acelaşi timp, presiunea generată împiedică curgerea sângelui prin arborele arterial coronarian. Deoarece contracţia miocardului este strâns legată de debitul coronarian şi de oxigenul livrat, echilibrul dintre necesarul de oxigen şi oxigenul livrat reprezintă un factor critic pentru funcţionarea normală a inimii.

De-a lungul unui ciclu cardiac debitul coronarian arterial este în anti-fază cu debitul coronarian venos. În timpul sistolei, contracţiile miocardului împiedică dezvoltarea unui debit arterial normal, dar amplifică debitul venos. În timpul diastolei, ca urmare a relaxării miocardului, debitul arterial creşte iar debitul venos scade.

Spre deosebire de celelalte componente ale sistemului cardiovascular, gradul de extracţie a oxigenului din sânge, în starea de repaus, este aproape de valoarea maximă (~75%). Datorită acestui aspect, o creştere a necesarului de oxigen nu poate fi satisfăcută decât printr-o creştere a debitului alimentat. Mărimile principale care determină necesarul de oxigen sunt frecvenţa cardiacă, presiunea sistolică şi contractilitatea ventriculului stâng. Dublarea valorii oricăruia dintre aceşti parametri conduce la o creştere cu aproximativ 50% a necesarului de oxigen. În vederea diagnosticării stenozelor (a determinării stenozelor semnificative şi nesemnificative), în


prezent se utilizează mai multe metode. La nivel mondial cea mai utilizată metodă este angiografia coronariană, prin intermediul căreia se estimează gradul de îngustare al stenozei (v. fig. 5.3):

[%]100N

S

D

DQCA , (5.1)

unde, DS reprezintă diametrul minim în zona stenozei, iar DN reprezintă diametrul normal (sănătos) al arterei coronariene. Dacă QCA este mai mare de 50%, se consideră că stenoza este semnificativă, iar dacă este mai mic de 50%, se consideră că stenoza este nesemnificativă. Această metodă are însă o serie de dezavantaje, printre care faptul că nu se iau în considerare aspectele fiziologice → metoda realizează doar o evaluare morfologică, nu şi funcţională a stenozelor;

Incapacitatea indicatorului QCA de a diagnostica corect stenozele coronariene a fost demonstrată recent într-un studiu [De Bruyne et al., 2012]. Deoarece, dezavantajele menţionate mai sus erau cunoscute şi anterior acestui studiu, de-a lungul ultimelor două decenii s-a introdus şi dezvoltat un nou indicator diagnostic al stenozelor: rezervă fracţionară de debit (FFR – Fractional Flow Reserve) [Pijls et al., 2000].

Fig. 5.3 Estimarea gradului de îngustare al unei stenoze coronariene prin intermediul angiografiei coronariene.

Determinarea indicatorului FFR presupune realizarea unei proceduri invazive, de-a lungul căreia se introduce un cateter (prin artera femorală sau radială) împreună cu un senzor de presiune, care este apoi avansat în arterele coronariene ale epicardului până la locaţia stenozei. Astfel, se determină presiunile medii în amonte şi în aval de stenoză. Aceste măsurători se realizează după administrarea unei substanţe care induce hiperemie (debit maxim), de obicei, folosindu-se adenozina. Indicatorul FFR reprezintă raportul dintre debitul maxim printr-o arteră stenozată ( ) şi debitul maxim ipotetic obţinut în absenţa stenozei ( ): SQmax

NQmax

N

S

Q

QFFR

max

max . (5.2)

Pe baza echivalentului legii lui Ohm în hidraulică, se cunoaşte faptul că debitul este raportul dintre presiune şi rezistenţă. Astfel:

N

maxva

Smaxvd

Nmax

Smax

R/PP

R/PP

Q

QFFR

, (5.3)

unde, Pd este presiunea medie măsurată în aval de stenoză, Pa este presiune aortică medie, Pv este presiunea venoasă, este rezistenţa microvasculară în prezenţa stenozei, în starea de S

maxR

27


hiperemie, iar este rezistenţa microvasculară ipotetică în absenţa stenozei, de asemenea, în

starea de hiperemie. Presupunându-se că cele două rezistenţe sunt egale şi că presiunea venoasă este neglijabilă în raport cu cele două presiuni arteriale, indicatorul FFR va fi determinat astfel [Pijls et al., 1996]:

28

NmaxR

a

dNmax

Smax

Q

Q

P

PFFR . (5.4)

Valoarea normală pentru indicatorul FFR este 1.0 iar valoarea critică este de aproximativ 0.8. Dacă indicatorul FFR se află sub această valoare, atunci se va implanta un stent.

Considerând, pe de o parte, importanţa clinică a indicatorului FFR în diagnosticarea stenozelor coronariene pe de o parte şi costul ridicat al unei măsurători a indicatoruluui FFR, dar şi riscurile asociate cu această intervenţie invazivă, pe de altă parte, modelul multiscalar al circulaţiei coronariene prezentat în acest capitol are rolul de a evalua pe cale non-invazivă valoarea indicatorului FFR.

5.2 MODEL MULTISCALAR AL CIRCULAŢIEI CORONARIENE

În continuare, se introduce un model multiscalar de ordin redus al circulaţiei coronariene. În cadrul acestui model, prezentat în figura 5.4 [Itu et al., 2012(b)], aorta precum şi arterele mari ale circulaţiei coronariene sunt modelate prin intermediul modelului cvasi unidimensional introdus în capitolul anterior. Circulaţia microvasculară este modelată prin intermediul unui model cu parametri distribuiţi specializat, care reprezintă condiţia de frontieră de ieşire pentru modelul unidimensional. Peretele arterial este modelat prin intermediul unei legi pur elastice, valorile celor trei constante din ecuaţia (3.9) fiind preluate din lucrarea [Mynard, 2011]. După cum a fost prezentat în secţiunea anterioară, contracţiile inimii au o influenţă importantă asupra caracteristicilor hemodinamice ale circulaţiei coronariene. Prin urmare, modelul multiscalar propus conţine şi un model simplificat al inimii, bazat pe modelul elastanţei variabile în timp. Acest model este util din două puncte de vedere: pe de o parte, furnizează condiţia de frontieră de intrare pentru aortă iar, pe de altă parte, permite estimarea presiunii intraventriculare care conduce la apariţia unui debit coronarian redus de-a lungul sistolei. Modelul inimii este cuplat cu frontiera de intrare a aortei prin intermediul unui model cu parametrii distribuiţi al valvei aortice.

Fig. 5.4 Model multiscalar al circulaţiei coronariene.


Contracţiile ventriculului stâng acţionează, în special, asupra circulaţiei microvasculare şi prin urmare, pentru modelarea microcirculaţiei coronariene nu se poate folosi modelul Windkessel clasic cu trei elemente. În literatura de specialitate s-a introdus un model cu parametrii distribuiţi specializat pentru circulaţia coronariană [Mantero et al., 1992], care este prezentat în figura 5.4 şi care este folosit în cadrul modelului multiscalar propus. În vederea modelării contracţiilor ventriculului, acesta înglobează presiunea intraventriculară. Modelul multiscalar modelează atât arborele coronarian stâng (LCA – Left Coronary Artery) cât şi arborele arterial drept (RCA – Right Coronary Artery). Deoarece arterele arborelui coronarian stâng alimentează, în principal, miocardul ventriculului stâng, presiunea intraventriculară folosită în cadrul modelelor microvasculare cuplate la frontierele de ieşire ale acestui arbore este aceea a ventriculului stâng. Similar, în cadrul modelelor cuplate la frontierele de ieşire ale arborelui coronarian drept se foloseşte presiunea intraventriculară a ventriculului drept.

Scopul principal al modelului multiscalar propus nu este să se simuleze hemodinamica circulaţiei coronariene doar în cazul unui arbore sănătos, ci mai ales în cazul unui arbore care prezintă stenoze coronariene în arterele epicardului. Una dintre supoziţiile care stau la baza deducerii modelului unidimensional, atunci când se porneşte de la modelul tridimensional, este că viteza axială este dominantă în raport cu viteza radială. Această supoziţie este valabilă însă numai în cazul unui arbore arterial sănătos, fără stenoze. În cazul unor modificări semnificative ale razei, aşa cum se întâmplă în cazul stenozelor, componenta radială a vitezei nu mai poate fi neglijată. Cel mai important aspect este căderea de presiune de-a lungul stenozei. În acest sens, pentru segmentele arteriale care conţin stenoze, modelul unidimensional este înlocuit cu un model analitic al căderii de presiune. Detaliile acestuia precum şi cuplarea cu modelul unidimensional sunt descrise în secţiunile următoare.

În vederea evaluării indicatorului diagnostic FFR trebuie modelată atât starea de repaus cât şi starea de hiperemie. Prin urmare, atunci când modelul multiscalar este construit pentru o geometrie specifică unui pacient, condiţiile de frontieră ale modelului multiscalar trebuie personalizate atât pentru starea de repaus cât şi pentru starea de hiperemie.

5.2.3 Model simplificat al inimii şi cuplarea cu circulaţia arterială

Modelul simplificat al inimii, bazat pe elastanţa variabilă în timp a fost menţionat pe scurt în secţiunea 2.1.4. Ideea de bază a acestui model este aceea că presiunea şi volumul ventriculelor sunt cuplate printr-un coeficient variabil în timp, numit elastanţă. Ecuaţia fundamentală a modelului este reluată în această secţiune:

0VtV

tPtE

LV

LV

, (5.26)

unde, E(t) este elastanţa, PLV(t) şi VLV(t) sunt presiunea şi respectiv volumul ventriculului stâng, iar V0 este volumul mort.

Elastanţa are valori ridicate în timpul sistolei şi valori scăzute pe durata diastolei. În lucrări din literatura de specialitate s-a arătat că elastanţa este practic independentă de sarcina aplicată ventriculului [Senzaki et al., 1996]. Acest model prezintă o serie de avantaje: permite modelarea conceptelor fundamentale ale celor două ventricule ale inimii, poate fi personalizat pe baza valorilor achiziţionate de la pacient şi poate fi cuplat cu modele ale circulaţiei arteriale de orice ordin (modele cu parametrii distribuiţi, modele 1D sau modele 3D) [Formaggia et al., 2006]. Deoarece, în realitate, inima este cuplată cu circulaţia arterială prin intermediul valvelor aortice, pentru cuplarea modelului inimii cu un model al circulaţiei arteriale, se foloseşte şi un model al valvei aortice, compus dintr-o rezistenţă şi o inertanţă hidraulică:

dt

dQLtQRtVtEtP artVartVLVa , (5.27)

unde, Pa este presiunea arterială, RV-art este rezistenţa valvei, iar LV-art este inertanţa valvei.

29


30

5.4 PERSONALIZAREA PARAMETRILOR MODELULUI MULTISCALAR

În vederea utilizării modelului multiscalar dezvoltat pentru simulări hemodinamice specifice pacienţilor, modelul trebuie personalizat. Trei aspecte principale trebuie personalizate pentru modelul multiscalar din figura 5.4: geometria coronariană; parametrii modelelor coronariene microvasculare; parametrii modelului inimii.

Pentru studiul prezentat în acest capitol geometria coronariană este achiziţionată prin tomografie computerizată.

5.4.1 Personalizarea parametrilor modelelor coronariene microvasculare

Deoarece indicatorul diagnostic FFR este evaluat în practică în starea de hiperemie a pacientului, în vederea evaluării non-invazive a acestui indicator şi simularea hemodinamică a circulaţiei coronariene trebuie să corespundă tot stării de hiperemie a pacientului. În laboratorul de cateterizare, unde se efectuează măsurătorile FFR, pacientul este adus în starea de hiperemie prin administrarea unui medicament vasodilator, înainte de a se evalua valoarea FFR-ului prin măsurarea presiunilor. Hiperemia poate fi indusă prin administrarea substanţei vasodilatoare pe cale intravenoasă (efect de durată lungă) sau intracoronariană (efect de durată scurtă). Administrarea intravenoasă a vasodilatorului conduce la o uşoară creştere a frecvenţei cardiace şi la o scădere a presiunii sângelui [Pijls et al., 2000]. Deoarece, pentru o simulare, efectul administrării intracoronariene a vasodilatorului poate fi extins la infinit şi această alternativă de a obţine hiperemia nu influenţează deloc, sau doar foarte puţin frecvenţa cardiacă şi presiunea sanguină şi este astfel mai uşor de modelat, descrierea care urmează se concentrează asupra acestei abordări. Deoarece FFR foloseşte valori medii ale presiunii (măsurate în amonte şi în aval faţă de stenoză), complianţele joacă un rol secundar, nefiind necesară estimarea precisă a acestora (deoarece influenţează doar forma variaţiei a presiunii şi debitului, nu însă şi valorile medii, care sunt determinate doar de rezistenţe).

5.4.1.1 Modelarea stării de repaus

În vederea modelării stării de repaus se poate folosi fie produsul frecvenţă-presiune fie produsul stres-masă-frecvenţă.

Utilizarea produsului frecvenţă-presiune

Presiunea aortică medie (PAM) poate fi estimată simplu folosind presiunea sistolică (PS) şi presiunea diastolică (PD), care pot fi măsurate non-invaziv:

)(3

1PDPSPDPAM . (5.32)

O expresie mai precisă a fost introdusă şi validată în lucrarea [Razminia et al., 2004], care ia în considerare şi frecvenţa cardiacă (HR), pornind de la ideea că durata relativă a perioadelor sistolice şi diastolice se modifică odată cu frecvenţa bătăilor:

PDPSHRPDPAM

0012.03

1 . (5.33)

Pentru a se estima rezistenţa coronariană, se poate folosi următoarea formulă (echivalentul legii lui Ohm în hidraulică):

correpaus RQPAM , (5.34)

unde, Qrepaus este valoarea medie a debitului coronarian, iar Rcor este rezistenţa coronariană totală, ambele în starea de repaus. Deoarece PAM este cunoscută, doar debitul coronarian mai trebuie determinat pentru a se calcula rezistenţa totală. După cum a fost precizat anterior, o


cerere de oxigen crescută poate fi satisfăcută doar printr-un debit crescut. Prin urmare, se poate considera că debitul coronarian este proporţional cu cererea de oxigen. Cererea şi consumul de oxigen în circulaţia coronariană sunt dificil de determinat prin măsurători şi pot fi obţinute doar invaziv. În literatura de specialitate au fost propuse o serie de metode pentru estimarea consumului de oxigen pe baza unor variabile mecanice. S-a arătat astfel că frecvenţa cardiacă este determinantul principal al consumului de oxigen. Al doilea determinant major este presiunea (generarea presiunii consumă mai mult oxigen decât contracţia muşchilor, i.e. decât generarea debitului). Prin urmare, cel mai utilizat index pentru estimarea consumului miocardic de oxigen este produsul frecvenţă-presiune. Utilizând acest produs, debitul relativ în starea de repaus poate fi estimat astfel [Anderson et al., 2000]:

}4.0]10)(7.0{[8 3 PSHRqrepaus . (5.35)

Această expresie este valabilă numai dacă debitul coronarian reuşeşte să satisfacă cererea de oxigen a pacientului. Teoretic, cererea miocardică de oxigen ar trebui iniţial calculată folosind produsul frecvenţă-presiune şi apoi ar trebui stabilită o corespondenţă între cererea de oxigen şi debit. Deoarece interesul este însă îndreptat către debit, ecuaţia (5.35) este de preferat, aceasta stabilind o relaţie directă între produsul frecvenţă-presiune şi debitul relativ.

Pentru a se determina valoarea absolută a debitului în starea de repaus, qrepaus, trebuie înmulţit cu masa miocardului. Pentru o inimă normală, se presupune, în general, că ventriculul stâng reprezintă aproximativ două treimi din masa totală, în timp ce, ventriculul drept şi atriile reprezintă cealaltă treime [Jauhiainen et al., 2002]. În concluzie, trebuie determinată masa ventriculului stâng MVS, şi apoi se poate determina valoarea absolută a debitului coronarian în starea de repaus, după cum urmează:

VSrestrest MqQ 5.1 . (5.36)

Prin urmare, rezistenţa coronariană totală poate fi determintă astfel:

repauscor Q

PAMR . (5.37)

Pentru geometria coronariană specifică unui anumit pacient, scopul este de a distribui această rezistenţă între diverse modele cu parametrii distribuiţi impuse la frontierele de ieşire ale geometriei. În acest sens, se poate folosi legea lui Murray [Murray, 1926(a)], [Murray, 1926(b)], care presupune că energia necesară pentru curgerea sângelui şi energia necesară pentru a menţine structura vaselor de sânge este minimă şi, prin urmare, se obţine:

3~ ii rkQ , (5.38)

unde, k este o constantă iar r este raza vasului de sânge.

5.4.1.2 Modelarea stării de hiperemie

Există diverse substanţe vasodilatoare care pot fi utilizate în laboratorul de cateterizare pentru evaluarea valorii FFR: dipyridamol, papaverină, adenozină, etc. Substanţa cea mai utilizată la momentul actual este adenozina, care este administrată, în mod preferenţial, pe cale intracoronariană. Prin urmare, următoarele paragrafe se concentrează asupra efectului intracoronarian al administrării adenozinei.

Efectul administrării adenozinei

După cum este descris în lucrarea [Wilson et al., 1990], adenozina conduce la o creştere a vitezei fluxului sanguin coronarian de aproximativ 4.5 ori pentru persoane normale, sănătoase (fără afecţiuni coronariene). Raportul dintre viteza sanguină la starea de repaus şi viteza la starea de hiperemie se notează în mod curent cu CFVR şi a fost determinat după administrarea intravenoasă sau intracoronariană a adenozinei. Valoarea de 4.5 a fost confirmată în toate

31


experimentele efectuate. Creşterea vitezei coronariene este egală cu creşterea debitului, deoarece se pot presupune profiluri de viteză asemănătoare. Deoarece, în starea de hiperemie, presiunea arterială scade uşor, o creştere de 4.5 ori a debitului nu înseamnă şi o scădere de 4.5 a rezistenţei. După cum este descris în lucrarea [Wilson et al., 1990], se poate calcula un index al rezistenţei coronariene totale (TCRI) după urmează:

32

repauscor

hipercor

repaus

repaus

hiper

hiper

R

R

Q

PAM

Q

PAMTCRI . (5.49)

Valorea medie a indexului TCRI determinată în studiul menţionat mai sus este 0.22. Unul din motivele pentru care CFVR nu are o valoare critică clară atunci când este folosit

pentru diagnosticarea funcţională a stenozelor este faptul că viteza în starea de repaus variază în mod semnificativ la diverse momente de timp (nu există o stare de repaus clar definită). Prin urmare, nu este posibil să se folosească valoarea 0.22 la calculul rezistenţei hiperemice pentru fiecare pacient. După cum a fost precizat anterior, frecvenţa cardiacă este determinantul principal al consumului de oxigen (împreună cu presiunea sistolică). De aceea, în studiul descris în lucrarea [Wilson et al., 1990], una dintre condiţii era ca frecvenţa cardiacă să fie mai mică de 80 bpm. Prin urmare, trebuie determinat dacă diferiţi parametrii precum frecvenţa cardiacă sau presiunea arterială influenţează sau nu valoarea CFR şi prin urmare valoarea TCRI.

5.4.1.3 Estimarea rezistenţelor microvasculare de repaus

Pentru estimarea rezistenţelor microvasculare de repaus se propune o metodologie formată din cinci paşi [Sharma et al., 2012]:

1. Presiunea arterială medie este calculată cu ecuaţia (5.33). 2. Masa ventriculului stâng este calculată cu ecuaţia (5.45). 3. Debitul coronarian total în starea de repaus este determinat cu ecuaţiile (5.35) şi (5.36)

dacă doar HR şi PS sunt cunoscute. Dacă este disponibil şi stresul (calculat conform ecuaţiei (5.46)), atunci se va folosi ecuaţia (5.47).

4. Deoarece debitul este proporţional cu raza la puterea a treia, debitul total de repaus, care este suma debitelor tuturor vaselor coronariene terminale, poate fi calculat astfel:

n

ii

n

iirest QrkQ

11

3 . (5.53)

5. Pasul final este de a calcula rezistenţa terminală totală a fiecărui vas de sânge, folosindu-se ecuaţia:

ii Q

PAMR . (5.54)

Pentru a determina Qi, ecuaţia (5.38) este împărţită la ecuaţia (5.53):

n

jj

in

jj

i

rest

i

r

r

rk

rk

Q

Q

1

3

3

1

3

3

, (5.55)

şi, prin urmare:

n

jj

iresti

r

rQQ

1

3

3

. (5.56)

Astfel, rezistenţa terminală totală a fiecărui vas de sânge este calculată, după cum urmează:


33

3

1

3

irest

n

jj

ii rQ

rPAM

Q

PAMR

. (5.57)

5.4.1.4 Estimarea rezistenţelor microvasculare hiperemice

Informaţia de intrare pentru cea de-a doua parte a metodei propuse este reprezentată de rezistenţele microvasculare de repaus, calculate cu ecuaţia (5.57). Se propun următorii paşi pentru calculul rezistenţelor microvasculare hiperemice:

1. Se calculează valoarea medie a vitezei maxime pe o secţiune transversală în starea de repaus (rAPV), folosind ecuaţia (5.48).

2. Se calculează valorile CFVR pentru fiecare arteră coronariană principală folosind ecuaţiile (5.50) – (5.52).

3. Se calculează indexul TCRI folosind o valoare de 5 mmHg pentru [Wilson et al., 1990]:

PAM

.1

1

CFVRPAM

PAMPAM

CFVRPAM

PAM

Q

Q

PAM

PAM

Q

PAM

Q

PAMTCRI

repaus

repaus

repaus

hiper

hiper

repaus

repaus

hiper

repaus

repaus

hiper

hiper

(5.58)

Alternativ, în locul paşilor 1-3, se poate folosi o expresie pentru determinarea valorii TCRI, bazată pe rezultatele experimentale descrise în lucrarea [McGinn et al., 1990], rezultate care au o deviaţie standard foarte mică (0.01 – 0.02):

bpm.10016.0001.0

bpm;1001.00016.0

HRpentruHR

HRpentruHRTCRI (5.59)

4. Se calculează rezistenţele microvasculare hiperemice folosind următoarea expresie:

TCRIRR repausihiperi . (5.60)

unde, este valoarea calculată cu ecuaţia (5.57). repausiR

5.4.2 Personalizarea modelului inimii

Pe baza informaţiilor prezentate în secţiunea 5.2.3 se poate defini un model generalizat al inimii. În vederea utilizării unui model personalizat al inimii, denormalizarea elastanţei prin algoritmul 5.1 trebuie realizată cu valori maxime şi minime ale elastanţei specifice pacientului. Pentru determinarea acestor valori personalizate se foloseşte procedura descrisă în continuare.

Pentru aplicarea acestei metode trebuie determinat momentul de timp la care se atinge valoarea maximă a elastanţei, tmax. Acesta este calculat pe baza unei expresii liniare, pornind de la frecvenţa cardiacă în starea de repaus a pacientului [Kim et al., 2010]:

,120,/15.0

;120,17.0/116.0

HRdacaHR

HRdacaHRtmax (5.61)

unde, HR este frecvenţa cardiacă. În figura 5.8 se prezintă variaţia tipică a elastanţei normalizate. Pentru a se deduce valoarea maximă a elastanţei denormalizate se porneşte de la următoarele două supoziţii: elastanţa maximă are o variaţie liniară de la un ciclu cardiac la altul; punctul de intersecţie al dreptei care uneşte punctele corespunzătoare elastanţei maxime

este fix (V0-SB). Folosind ecuaţia (5.26), se poate scrie pentru două puncte de pe graficul P-V:


34

SBNNN VtVtPtE 0/

, (5.62)

SBmaxmaxmaxmax VtVtPEtE 0/ , (5.63)

unde, ecuaţia (5.62) reprezintă un punct oarecare de pe grafic, iar ecuaţia (5.63) reprezintă punctul corespunzător elastanţei maxime. Efectuând raportul celor două ecuaţii, se poate scrie:

SBN

SBmaxNNNN

SBN

SBmax

max

N

max

N

VtV

VtVtPtE

VtV

VtV

tP

tP

tE

tE

0

0

0

0 . (5.64)

Prin urmare, se poate determina V0-SB (volumul mort – volumul de sânge rămas în ventricul la finalul sistolei):

NNNN

NNNmaxNNSB tEtP

tVtEtVtPV

0 . (5.65)

Fig. 5.8 Elastanţa normalizată (a); Variaţia presiune-volum de-a lungul unui ciclu cardiac în interiorul ventriculului [Senzaki et al., 1996] (b).

Punctul corespunzător ecuaţiei (5.62) a fost ales ca fiind punctul corespunzător momentului de timp la care începe sistola. Prin urmare VFDtV N , VFD fiind volumul

intraventricular la finalul diastolei, iar PDtP N este presiunea diastolică aortică. Presiunea

corespunzătoare punctului de elastanţă maximă a fost estimată pe baza unei expresii validate anterior [Kelly et al., 1992]:

PS.PtP vSmax 90 , (5.66)

unde, PvS este presiunea maximă din interiorul ventriculului, iar PS este presiunea sistolică aortică. Prin urmare, folosind expresia maxNNN tP/tPtP ecuaţia (5.65) poate fi rescrisă

astfel:

NN

NNSB tEPSPD

tEVFDPSVFSPDV

9.0/

9.0/0 , (5.67)

unde, VFS este volumul intraventricular la finalul sistolei, care poate fi calculat după cum urmează:

60/HRVBVFDVFS , (5.68)

unde, VB este volumul de bătaie al unui ciclu. În final, se calculează valoarea maximă a elastanţei:


35

SBmax VVFSPSE 0/9.0 . (5.69)

5.5 PROCEDURĂ DE OPTIMIZARE PENTRU CIRCULAŢIA CORONARIANĂ

În secţiunea precedentă au fost prezentate o serie de activităţi de personalizare a parametrilor modelului multiscalar. Pentru ca o simulare hemodinamică să furnizeze rezultate care să corespundă îndeaproape sistemului cardiovascular al unui pacient, este nevoie ca, suplimentar, anumite mărimi simulate să aibă aceleaşi valori ca şi mărimile achiziţionate de la pacient. În acest sens, în figura 5.9 se prezintă o vedere de ansamblu asupra procesului necesar a fi realizat pentru estimarea valorilor FFR ale stenozelor prezente în geometria coronariană: Procesarea datelor achiziţionate prin CT: geometria achiziţionată prin tomografie

computerizată şi descrisă prin axele centrale ale vaselor precum şi prin contururile secţiunilor transversale este preprocesată pentru a se construi modelul arterial 1D. Se identifică segmentele stenozelor;

Optimizarea cu modelul 0D: pe baza geometriei obţinute în cursul paşilor precedenţi, se construieşte un model detaliat cu parametrii distribuiţi (0D) compus dintr-un model al inimii, un model al valvei precum şi din circuite RLC. Metoda dogleg trust region, descrisă ulterior, este folosită pentru a adapta presiunea arterială medie şi procentajul reprezentat de debitul coronarian din debitul cardiac total. De-a lungul procesului de optimizare se adaptează volumul diastolic al ventriculului stâng şi rezistenţa sistemică totală (rezistenţa impusă la frontiera de ieşire a aortei). Rezistenţa sistemică totală determinată în această etapă este finală şi nu va fi modificată în paşii următori. Valoarea estimată a volumului ventriculului stâng trebuie adaptată, în continuare, deoarece modelul cu parametrii distribuiţi nu este capabil să surprindă fenomenele de propagare a undelor;

Fig. 5.9 Procesul de estimare a presiunilor coronariene şi a valorilor FFR.


36

Utilizarea modelului cu geometrie 1D redusă: este construită o geometrie vasculară 1D redusă, compusă doar din segmentul aortic şi din două vase lungi de 4.0 cm corespunzătoare arborilor coronarieni stâng şi respectiv drept. Metoda dogleg trust region şi o metodă cvasi-Newton sunt utilizate în această etapă pentru optimizarea efectuată în vederea respectării valorilor de referinţă ale presiunilor medii, maxime şi minime ale pacientului. Astfel, volumul diastolic al modelului inimii este adaptat, în continuare, împreună cu rezistenţa proximă şi rezistenţa distantă, dar şi cu complianţa totală a condiţiei de frontieră sistemice (modelul Windkessel cu 3 elemente). Toţi parametri modificaţi în această etapă vor fi adaptaţi şi la pasul cinci folosind modelul 1D cu geometrie completă;

Iniţializarea variabilelor modelului 1D cu geometrie completă pe baza rezultatelor obţinute cu modelul 1D cu geometrie redusă: variabilele aplicaţiei de calcul a valorilor FFR (debitul, ariile transversale şi presiunea din cadrul segmentelor 1D precum şi debitul din cadrul modelelor microvasculare coronariene) sunt iniţializate pe baza soluţiilor obţinute cu geometria 1D redusă. Astfel, nu mai este necesară o perioada de iniţializare (care necesită simularea a 5-7 cicluri cardiace suplimentare) atunci când se porneşte simularea cu modelul 1D cu geometrie completă;

Utilizarea modelului 1D cu geometrie completă pentru simularea stării de repaus: simularea cu geometria completă, corespunzătoare stării de repaus este pornită cu valorile parametrilor determinate în cadrul etapelor anterioare. În această etapă, cei trei parametri adaptaţi în cadrul procesului de optimizare realizat cu modelul 1D cu geometrie redusă, sunt modificaţi în continuare pentru a respecta cele trei presiuni de referinţă. Când sunt îndeplinite condiţiile de convergenţă, rezistenţele microvasculare coronariene, corespunzătoare stării de repaus, sunt înlocuite cu valorile corespunzătoare stării de hiperemie;

Utilizarea modelului 1D cu geometrie completă pentru simularea stării de hiperemie: în cadrul acestei etape toate valorile parametrilor determinate în cadrul etapelor anterioare sunt menţinute constante şi nu se mai utilizează procedura de optimizare. Simularea este rulată până la apariţia convergenţei, sunt înregistrate valorile presiunii din cadrul arborilor coronarieni şi sunt calculate valorile FFR corespunzătoare segmentelor de stenoză;

Includerea valorilor de presiune şi de FFR în fişierul XML: rezultatele obţinute în etapele anterioare sunt incluse în fişierul XML iniţial care conţine geometria arborilor coronarieni şi care este salvat cu denumirea ‘results.xml’.

5.5.1 Optimizare cu modelul 0D detaliat

În continuare, atenţia este îndreptată spre pasul al doilea din figura 5.9. Iniţial, se construieşte un model 0D detaliat compus dintr-un model al inimii şi al valvelor, circuite RLC pentru arterele coronariene proxime, şi modele microvasculare pentru arteriole şi capilare (v. fig. 5.10) [Milisic et al., 2004].

Primul obiectiv al acestui model este să se obţină o anumită valoare pentru procentajul reprezentat de debitul coronarian din debitul cardiac total. Doar valorile rezistenţelor prezintă importanţă pentru distribuţia debitului mediu între artere iar valorile parametrilor modelelor microvasculare sunt constante şi sunt setate conform procedurii descrise în secţiunea 5.4.1. Prin urmare, pentru îndeplinirea primului obiectiv, se va adapta rezistenţa sistemică totală.

Al doilea obiectiv al acestui model este de a se obţine o anumită valoare medie a presiunii arteriale. Pentru a determina parametrul care va fi adaptat în vederea îndeplinirii acestui obiectiv, s-a realizat o analiza de sensibilitate, care a demonstrat faptul că există diverşi parametri care permit modificarea debitului cardiac (şi, prin urmare, a presiunii arteriale) al modelului redus al inimii cuplat la frontiera de intrare a modelului arterial: momentul de timp la care se atinge valoarea maximă a elastanţei; elastanţa maximă; volumul mort al inimii; volumul ventriculului la finalul diastolei; rezistenţa sistemică; presiunea atriumului stâng (frecvenţa cardiacă este fixată şi nu poate fi modificată). Cea mai mare sensibilitate s-a obţinut pentru volumul ventriculului la finalul diastolei, parametru care a fost, în continuare, adaptat pentru a îndeplini cel de al doilea


obiectiv. În studii anterioare [Coogan et al., 2011], presiunea atriumului stâng a fost adaptată dar, deoarece rolul modelului inimii este doar de a genera un anumit debit şi o anumită presiune şi de a furniza o estimare a presiunii intramiocardice (care este dată de ventricul), atriumul nu trebuie modelat.

În timp ce valoarea rezistenţei sistemice, determinate la acest pas, este finală, valoarea volumului diastolic al ventriculului stâng reprezintă doar o estimare iniţială, care va fi adaptată în continuare.

Procedura de determinare a valorilor celor doi parametri a fost formulată ca o problemă de determinare a soluţiei unui sistem de ecuaţii neliniare cu o rădăcină în punctul în care valorile simulate sunt egale cu valorile de referinţă. Valoarea de referinţă a procentajului reprezentat de debitul coronarian din debitul total a fost setată la 4.5% [Pijls et al., 2000], în timp ce valoarea de referinţă a presiunii medii este calculată pe baza presiunii diastolice şi sistolice şi a frecvenţei cardiace (toate aceste mărimi sunt măsurate direct la pacient), conform ecuaţiei (5.33):

Fig. 5.10 Model 0D detaliat folosit în timpul procedurii de optimizare de la pasul doi din figura 5.9 (folosind metoda dogleg trust region).

Pentru a se determina valorile reziduurilor (ale funcţiilor obiectiv: f(xi)) pentru un anumit set de valori ale parametrilor xi, trebuie realizată o simulare cu modelul 0D. Sistemul neliniar de ecuaţii a fost formulat după cum urmează:

0

0

refsim

refmedsimmed

t

D

ProcQProcQ

PP

R

Vf , (5.70)

unde, VD este volumul ventriculului stâng la finalul diastolei, Pmed este presiunea arterială medie, ProcQ este procentajul reprezentat de debitul coronarian, sim face referire la valorile

determinate prin intermediul simulării cu modelul 0D, iar ref face referire la valorile de

referinţă.

37


38

5.5.2 Optimizare cu modelul 1D cu geometrie redusă

Odată ce este îndeplinită condiţia de convergenţă pentru procesul de optimizare din figura 5.11, se trece la pasul trei al procesului prezentat în figura 5.9. Acest pas este format din metoda dogleg trust region şi dintr-o metodă cvasi-Newton, operaţiile individuale fiind prezentate în figura 5.12. Algoritmul de optimizare este bazat pe metodologia introdusă în lucrarea [Spilker et al., 2010].

Fig. 5.12 Schema - bloc a procedurii de optimizare bazată pe metoda dogleg trust region (în partea stângă) şi pe o metodă cvasi-Newton (în partea dreaptă).

De-a lungul acestui pas se foloseşte modelul cu geometrie redusă prezentat în figura 5.13(b). Acesta este compus dintr-un model redus al inimii (bazat pe modelul elastanţei variabile), un model al valvei aortice, partea proximă a aortei, precum şi două vase lungi de 4.0 cm ale arterei coronariene stângi, respectiv drepte. Pentru fiecare dintre aceste trei artere, la impunerea condiţiei de frontieră de ieşire, se folosesc modele cu parametrii distribuiţi. Procedura


de optimizare se concentrează asupra volumului de sânge din ventriculul stâng la finalul diastolei şi asupra modelului microvascular impus la frontiera de ieşire a aortei, adică asupra rezistenţelor proxime şi distante Rp şi Rd, precum şi asupra complianţei C. Scopul procedurii de optimizare este acela de a se realiza o simulare hemodinamică ale cărei presiuni medii, maxime şi minime să fie aceleaşi cu cele ale pacientului.

Rezistenţa aortică totală, Rt = Rp + Rd, a fost determinată anterior folosindu-se un model detaliat cu parametrii distribuiţi şi este menţinută constantă de-a lungul optimizării (pentru a se obţine un anumit procentaj al debitului coronarian din debitul total). Prin urmare, pentru a se respecta cele trei presiuni de referinţă, se vor adapta următorii trei parametri: VD – volumul ventriculului stâng la finalul diastolei; ρ = Rp / Rt – raportul dintre rezistenţa proximă şi rezistenţa totală; τ = CRd – constanta de timp corespunzătoare scăderii exponenţiale a presiunii

elementului Windkessel în cazul în care debitul este zero. La fel ca şi pentru pasul doi din figura 5.9, procedura de determinare a valorilor celor trei

parametri a fost formulată ca o problemă de determinare a soluţiei unui sistem de ecuaţii neliniare cu o rădăcină în punctul în care valorile simulate ale presiunii sunt egale cu valorile de referinţă (valorile măsurate la pacient). Pentru determinarea valorilor reziduurilor (ale funcţiilor obiectiv f(xi)) pentru un anumit set al valorilor parametrilor xi, este nevoie să se efectueze o simulare cu geometria redusă din figura 5.13(b).

Sistemul neliniar f(xi) = 0 a fost rezolvat cu o metodă cvasi-Newton, care actualizează valorile parametrilor la fiecare iteraţie. Astfel, la iteraţia i, pasul si este folosit pentru a se actualiza valorile parametrilor xi, folosind ultima estimare a Jacobian-ului Ji al sistemului neliniar:

iiiii xfJxxs 11

. (5.80)

Capacitatea acestei metode de a actualiza parametri este independentă de diferenţele dintre valorile reale ale parametrilor şi ale funcţiilor obiectiv. Sistemul ecuaţiilor neliniare este formulat astfel:

0

0

0

refminsimmin

refmaxsimmax

refmedsimmeds

PP

PP

PPV

f . (5.81)

Fig. 5.13 Model cu parametrii distribuiţi utilizat pentru iniţializarea procedurii de

optimizare (folosind metoda dogleg trust region) (a); Model 1D cu geometrie redusă (b).

5.5.3 Optimizare cu modelul 1D cu geometrie completă

Odată ce condiţia de convergenţă a procesului de optimizare aplicat asupra modelului 1D cu geometrie redusă este îndeplinită şi modelul 1D cu geometrie completă a fost iniţializat pe baza valorilor obţinute cu modelul 1D cu geometrie redusă, se începe procesul de optimizare cu

39


geometria 1D completă, corespunzător stării de repaus a pacientului. Metoda cvasi-Newton utilizată anterior pentru geometria 1D redusă este aplicată şi pentru geometria 1D completă (procesul de optimizare utilizat este prezentat în figura 5.14), sistemul neliniar rezolvat prin această abordare fiind cel introdus prin ecuaţia (5.81). Matricea Jacobiană şi valorile parametrilor sunt iniţial setate pe baza ultimelor valori obţinute la procesul de optimizare desfăşurat pentru geometria 1D redusă. Din nou, testul de convergenţă al procesului de optimizare foloseşte praguri egale cu 1/10 din valorile tipice ale obiectivelor, iar ecuaţia (5.83) este folosită pentru actualizarea matricei Jacobiene, în timp ce ecuaţia (5.80) este folosită pentru calculul variaţiei parametrilor şi pentru actualizarea valorilor acestora.

Se subliniază faptul că fiecare simulare realizată cu modelul 0D, modelul 1D redus, respectiv modelul 1D complet este rulată până când este îndeplinit criteriul de convergenţă pentru fiecare obiectiv, care este dat de o variaţie mai mică de 0.01% a mărimii respective de la un ciclu cardiac la altul. Fiecare simulare ulterioară, cu noi valori ale parametrilor foloseşte ca valori iniţiale rezultatele obţinute pentru vechiul set de parametri.

Fig. 5.14 Schema - bloc a procedurii de optimizare bazată pe metoda cvasi-Newton, aplicată pentru modelul 1D cu geometrie completă.

5.6 REZULTATE

5.6.5 Utilizarea modelului multiscalar pentru diagnosticarea pacienţilor cu stenoze coronariene

În continuare, sunt prezentate rezultatele simulărilor hemodinamice efectuate prin intermediul modelului multiscalar introdus în acest capitol pentru geometriile coronariene a patru pacienţi. Accentul este pus asupra indicatorului diagnostic FFR care a fost măsurat pentru aceşti patru pacienţi şi a cărui valoare este comparată cu valoarea obţinută în cadrul simulării hemodinamice. Toate rezultatele prezentate corespund stării de hiperemie intracoronariană.

40


Primul pacient prezintă o stenoză coronariană în artera interventriculară anterioară cu o valoare măsurată a FFR-ului de 0.74. Deoarece valoarea este mai mică decât 0.8, stenoza este semnificativă din punct de vedere funcţional. Înainte de a prezenta valorile FFR obţinute în urma simulării, în figura 5.21 se prezintă variaţia în timp a presiunii aortice în comparaţie cu variaţia în timp a presiunii arteriale din aval de stenoză, precum şi variaţia în timp a debitului coronarian prin artera stenozată. Se poate observa că presiunea coronariană din aval de stenoză este semnificativ mai mică decât presiunea aortică ca urmare a influenţei stenozei. Diferenţa cea mai mare dintre presiunea aortică şi presiunea coronariană se obţine la începutul diastolei, când debitul coronarian are valorile cele mai mari. Debitul coronarian în starea de hiperemie prezintă o formă de undă tipică, cu valori scăzute de-a lungul sistolei şi valori ridicate de-a lungul diastolei.

Fig. 5.21 Presiunile aortice şi coronariene şi debitul coronarian al primului pacient analizat.

În figura 5.22 sunt reprezentate valorile indicatorului FFR în geometria coronariană. Se poate observa că valoarea FFR obţinută în urma simulării hemodinamice cu modelul multiscalar este de 0.726. Pe de o parte, această valoare este apropiată de valoarea măsurată (0.74) iar, pe de altă parte, clasificarea stenozei realizată ca urmare a valorii FFR simulate, este aceeaşi ca şi pentru valoarea măsurată, şi anume stenoza este semnificativă.

Fig. 5.22 Valorile indicatorului FFR obţinute în urma simulării

hemodinamice în arborii arteriali ai pacientului 1.

41


42

5.7 CONCLUZII

În acest capitol s-a introdus un model multiscalar al circulaţiei coronariene. Modelul este format din trei componente fundamentale: un model simplificat al inimii, modelul unidimensional introdus în capitolul trei precum şi un model microvascular al circulaţiei coronariene.

Scopul principal al modelului introdus este acela de a se realiza o diagnosticare non-invazivă a stenozelor coronariene. În varianta clasică de diagnosticare se foloseşte indicatorul QCA, care ţine cont de aspecte exclusiv anatomice, neglijând aspectele funcţionale ale circulaţiei coronariene. În ultimul deceniu s-a impus tot mai mult un alt indicator diagnostic, şi anume FFR, care prezintă avantajul că evaluează funcţional stenozele coronariene. Dacă valoarea acestuia este mai mică de 0.8, atunci se recomandă plasarea unui stent, iar dacă valoarea acestuia este mai mare de 0.8 se prescrie un tratament medicamentos. Dezavantajul acestui indicator este că necesită efectuarea unor măsurători invazive de presiune, care, pe de o parte, supun pacientul unui anumit iar, pe de altă parte, sunt costisitoare. Prin urmare, obiectivul modelului propus în cadrul capitolului a fost acela de evalua pe cale non-invazivă valoarea indicatorului FFR.

6. DIAGNOSTICARE HEMODINAMICĂ NON-INVAZIVĂ A COARCTAŢIEI AORTICE

6.1 INTRODUCERE

Coarctaţia aortei este un defect congenital care constă într-o îngustare a aortei, care apare la 5 până la 8% dintre toţi pacienţii cu afecţiunea cardiacă congenitală (v. fig. 6.1) [Ringel et al., 2012]. Nu există certitudine referitoare la factorii care cauzează apariţia acestei afecţiuni. Totuşi, o explicaţie general acceptată este aceea că patologia este legată de anormalităţi în procesul dezvoltării arterelor şi în hemodinamica punctului de ramificaţie al vasului ductus arteriorus în timpul dezvoltării fătului. Se pare că o serie de factori au influenţă, dintre care cel mai important este cel genetic [Völkl et al., 2005].

Datorită îngustării induse de coarctaţia aortei, presiunea sângelui în partea superioară a corpului şi în ventriculul stâng este mai mare. Acest aspect conduce la rândul lui la apariţia hipertrofiei miocardice pentru a se menţine la un nivel normal stresul lângă pereţii arteriali, dar şi debitul cardiac. Dacă hipertensiunea nu este tratată, se poate ajunge la insuficienţă cardiacă congestivă.

Pacienţii născuţi cu coarcţie aortică necesită îngrijire medicală continuă, care constă în vizualizarea invazivă sau non-invazivă a aortei, terapie medicamentoasă şi, în cazul în care coarctaţia se reformează, cateterizare invazivă sau intervenţie chirurgicală pentru a reduce presiunea sanguină în aorta ascendentă. Pentru evaluarea pre-operatorie a severităţii coarctaţiei şi pentru evaluarea post-operatorie a îngustării reziduale există o serie de tehnici utilizate în clinici. Evaluarea anatomică se bazează, de obicei, pe rezonanţă magnetică (MRI) sau tomografie computerizată (CT), în timp ce, evaluarea funcţională constă în măsurarea căderii de presiune (Δp) de-a lungul coarctaţiei. Cea mai precisă evaluare a căderii de presiunii (determinată ca diferenţa de presiune dintre valoarea sistolică maximă din amonte de coarctaţie şi valoarea sistolică maximă din aval de coarctaţie) este realizată invaziv prin măsurarea presiunilor cu ajutorul unui cateter. Dacă valoarea este mai mare de 20 mmHg, atunci se consideră că îngustarea este semnificativă.

6.2 METODĂ NON-INVAZIVĂ DE EVALUARE A COARCTAŢIILOR AORTICE

După cum a fost precizat şi anterior, pentru a se obţine o metodă fezabilă de determinare non-invazivă a căderii de presiune, timpul de execuţie al algoritmului este un parametru critic. Având în vedere acest aspect, s-a ales abordarea bazată pe modelul cvasi unidimensional, care împreună cu elementele Windkessel impuse la frontierele de ieşire ale geometriei reprezintă un model de ordin de redus al fluxului sanguin în cadrul aortei. Abordarea de ordin redus este cel


puţin câteva ordine de mărime mai rapidă decât calculele bazate un model tridimensional, putând, pe de altă parte, estima corecta presiunile şi debitele variabile în timp în cadrul unor modele specifice pacienţilor [Reymond et al., 2011]. Următorul pas în efectuarea unor simulări specifice pacienţilor a fost acela de a se dezvolta o procedură de estimare a parametrilor necesari personalizării modelului prin utilizarea unor informaţii obţinute pe cale non-invazivă de la pacient (această strategie de personalizare a fost introdusă într-o formă preliminară în lucrările [Ralovich et al., 2012(a)], [Ralovich et al., 2012(b)]). Estimarea este realizată în conjuncţie cu un model comprehensiv al căderii de presiune care este cuplat cu modelul de ordin redus al fluxului sanguin în vederea estimării căderii de presiune pentru fiecare pacient.

În cadrul metodei propuse, se impune o formă de undă de debit variabil în timp la frontiera de intrare în cadrul modelului (determinată din date 2D PC-MRI) şi modele Windkessel cu trei elemente la fiecare frontieră de ieşire (artera brahiocefalică, artera carotidă comună stângă, artera subclavie stângă şi aorta descendentă). Abordarea propusă poate fi aplicată atât pentru un model pre-operatoriu cât şi pentru un model post-operatoriu. În continuare, ne vom îndrepta atenţia însă asupra modelului pre-operatoriu, deoarece acesta este mai semnificativ din punct de vedere clinic. Figura 6.2 prezintă abordarea propusă pentru cazul pre-operatoriu.

Fig. 6.2 Prezentarea generală a metodei de personalizare pentru cazul pre-operatoriu.

Simularea fluxului sanguin în segmentele sănătoase ale geometriei aortice se realizează cu ajutorul modelului unidimensional prezentat în capitolul 3. Una dintre supoziţiile utilizate la deducerea modelului de ordin redus este aceea că viteza axială este dominantă şi, prin urmare, că vitezele pe direcţia radială sunt neglijabile. Această supoziţie este valabilă pentru vasele normale, sănătoase, dar în cazul unei modificări bruşte a diametrului, precum în cazul unei coarctaţii, componentele radiale nu mai pot fi excluse. Astfel, s-a introdus un model al căderii de presiune (descris în secţiunile următoare) de-a lungul coarctaţiei pentru a modela corect

43


rezistenţa introdusă de coarctaţie. La implementare, acest segment a fost cuplat cu segmentele proxime şi distante prin impunerea continuităţii presiunii totale şi a debitului. La frontierele de ieşire s-a aplicat ecuaţia Windkessel pentru a completa sistemul de ecuaţii.

6.2.1 Estimarea parametrilor pentru personalizarea modelului

Simulări specifice pacienţilor ale fluxului sanguin necesită condiţii fiziologice la frontierele de intrare şi de ieşire ale domeniului de calcul. În funcţie de disponibilitatea măsurătorilor in-vivo, precum şi de supoziţiile modelului, activităţile de cercetare raportate în literatură folosesc în mod tipic una din următoarele două tipuri de condiţii la frontiera de intrare: viteză (sau debit) variabil în timp [LaDisa et al, 2011], [Olufsen et al., 2000]; model simplificat al inimii [Formaggia et al., 2006] [Coogan et al., 2011].

Prima dintre cele două tipuri de condiţii poate fi determinată simplu într-un cadru clinic, făcând, de obicei, parte din fluxul de lucru clinic de rutină (2D/3D PC-MRI, ultrasunete Doppler). Aceste măsurători sunt apoi impuse la frontiera de intrare pe baza unui profil de viteze plat, parabolic sau de tip Womersley. Abordarea alternativă este de a cupla un model simplificat al circulaţiei din amonte (e.g. un model cu parametrii distribuiţi al inimii) şi de a adapta parametrii modelului pentru a obţine debite şi presiuni fiziologice în interiorul domeniului de interes. Pentru condiţia de frontieră de ieşire, elementele de tip Windkessel surprind aspectele fiziologice ale circulaţiei şi au fost folosite pe scară largă în cadrul studiilor raportate [Stergiopulos et al., 1992], [Vignon-Clementel, 2010]. Pe baza informaţiilor non-invazive trebuie astfel estimaţi trei parametri (două rezistenţe: rezistenţa proximă – Rp, cea distantă – Rd şi o complianţă – C) pentru fiecare frontieră de ieşire. Componentele principale ale metodologiei de personalizare sunt estimarea condiţiilor de frontieră la intrarea şi la ieşirile din domeniu, cuplarea unui model al căderii de presiune precum şi estimarea proprietăţilor mecanice ale pereţilor arteriali pe baza datelor achiziţionate de la pacient.

6.2.2 Modelul căderii de presiune şi estimarea condiţiilor de frontieră

Presiunea arterială medie (MAP), definită ca presiune medie de-a lungul unui ciclu cardiac este responsabilă de transportul sângelui către vasele distante ale circulaţie şi, în final, către ţesuturi. MAP este legată de rezistenţa arterială prin următoarea lege: , unde Q este debitul mediu într-un anumit punct al circulaţiei arteriale, iar R este rezistenţa arterială totală distantă. În cazul aortei, se poate scrie pentru fiecare frontieră de ieşire i:

RQMAP

iti RQMAP , (6.3)

unde, Qi este debitul mediu prin frontiera de ieşire i, iar itR este rezistenţa totală distantă, care

reprezintă suma celor două rezistenţe ale modelului Windkessel (Rt = Rp+ Rd). În cadrul aortei ascendente, MAP este estimată pe baza valorilor non-invazive determinate prin tensiometru [Razminia et al., 2004]:

),(0012.03

1DBPSBPHRDBPMAP

(6.4)

unde, HR este frecvenţa cardiacă, iar SBP şi DBP reprezintă presiunile sistolice şi diastolice. Valorile medii ale debitului în aorta ascendentă (Qasc) şi aorta descendentă (Qdesc) sunt

determinate din secţiunile PC-MRI. Astfel, debitul total către arterele supra-aortice poate fi determinat prin expresia Qsupra-aortic = Qasc – Qdesc.

Pentru primele ramuri care se bifurcă din aortă, debitul este distribuit către arterele copil direct proporţional cu raza la pătrat a vaselor respective [Zamir et al., 1992]. Astfel:

3

1

22 /i

iiaorticsuprai rrQQ , (6.5)

44


unde, ri este raza arterei supra-aortice i. Deoarece diferenţa de presiune dintre aorta ascendentă şi cele trei ramuri supra-aortice este nesemnificativă (pierderile vâscoase sunt neglijabile), aceeaşi presiune medie este folosită pentru estimarea rezistenţei totale:

i

it Q

MAPR . (6.6)

În cazul pacienţilor cu coarctaţie, supoziţia aceasta nu mai este valabilă pentru aorta descendentă deoarece îngustarea introdusă de coarctaţie introduce, la rândul ei, o cădere de presiune de-a lungul aortei, care poate fi transformată într-o rezistenţă dependentă de debit Rc(Q). Astfel, rezistenţa totală a aortei descendente, care reprezintă suma rezistenţelor coarctaţiei şi a modelului Windkessel impus la frontiera de ieşire, este estimată după cum urmează:

desc

cdesct Q

MAPQRR )( . (6.7)

Rezistenţa dependentă de debit este estimată pe baza unui model al căderii de presiunii. Pe baza acestor model, este propus următorul model comprehensiv al căderii de presiune

de-a lungul coarctaţiei [Itu et al., 2013(a)]:

qRKt

qLKqq

A

A

A

KqRKP vccuu

c

tvcv )α(||1

2

ρ)α(Δ

2

020

, (6.8)

unde, primul termen modelează pierderile vâscoase, al doilea termen modelează pierderile turbulente, al treilea termen modelează efectul inerţial iar al patrulea termen este o componentă

continuă. este un coeficient de viscozitate iar 20/053.01 AAK cv

cL

vc dllr

R0

4 )(

18

este

rezistenţa vâscoasă; este un coeficient de turbulenţă; 52.1tK 2.1uK este un coeficient de

inertanţă iar

cL

u dllr

L0

2

1

este inertanţa; este un coeficient de continuitate, iar

este numărul Womersley. Secţiunile de start şi de final ale coarctaţiei au fost considerate ca fiind locaţia unde raza scade sub 95% din valoarea de referinţă a locaţiei respective, şi respectiv locaţia unde raza creşte la peste 95% din valoarea de referinţă a locaţiei respective. Formulările termenilor de turbulenţă şi de inerţie din ecuaţia (6.8) au fost alese datorită posibilităţii de a ţine cont de forma specifică a coarctaţiei. Astfel, modelul căderii de presiune a putut fi personalizat pentru o geometrie a coarctaţiei specifică pacientului. Termenul de turbulenţă a fost folosit cu succes în diferite studii independente realizate in-vitro [Seeley et al., 1976] sau in-vivo [Steele et al., 2003]. De asemenea, s-a folosit un termen continuu, care a fost introdus anterior, ca urmare a diferenţei de fază dintre debit şi căderea de presiune, identificată într-un studiu computaţional [Bessems, 2007]. Deoarece modelul căderii de presiune conţine atât un termen liniar, cât şi un termen pătratic în raport cu debitul, s-au investigat două abordări diferite pentru evaluarea rezistenţei introduse de coarctaţie:

20018.0 cK

rezistenţa este calculată pe baza debitului mediu al aortei descendente:

descdescc QQPQR /)(Δ)(1 ; (6.9)

rezistenţa este calculată pe baza valorii medii a rezistenţelor calculate la fiecare moment de timp:

ntqtqPQRn

descdescc /)(/))((Δ)(1

2

, (6.10)

unde, este calculat cu ecuaţia (6.8), iar n este numărul cadrelor achiziţionate prin PC-MRI. )(P

45


Rezistenţa proximă la fiecare frontieră de ieşire este setată egală cu rezistenţa caracteristică a vasului (pentru a minimiza reflexiile) şi este calculată astfel:

ii

ip r

hE

rR

3

ρ2

π

12

, (6.11)

estimarea valorii E·h/ri fiind prezentată mai jos. În continuare, rezistenţa distantă este calculată prin scăderea rezistenţei proxime din rezistenţa totală.

Pentru estimarea complianţelor, se determină iniţial o complianţă totală [Stergiopulos et al., 1992]. În continuare, complianţa vaselor modelate ca segmente unidimensionale (Cprox) este calculată prin însumarea complianţelor de volum ale segmentelor corespunzătoare:

dl

hE

rAC j

jprox 2

3, (6.12)

unde, A este secţiunea transversală. În final, complianţa totală impusă la frontierele de ieşire (Cout) este determinată prin scăderea Cprox din Ctot. Valoarea astfel obţinută este distribuită către cele patru frontiere de ieşire:

3

1

22

iiioutiout rrCC . (6.13)

6.2.3 Estimarea proprietăţilor pereţilor arteriali

Un aspect important al unui studiu computaţional bazat pe pereţi arteriali complianţi este estimarea proprietăţilor mecanice ale acestora. În acest sens, s-a folosit o metodă care are la bază determinarea vitezei de undă [Olufsen et al., 2000], precum şi relaţia dintre viteza de undă şi proprietăţile peretelui arterial:

0ρ3

2

r

hEc

, (6.14)

unde, c este viteza de undă. Pentru estimarea vitezei de undă, s-a folosit metoda timpului de tranziţie [Ibrahim et al., 2010], unde txc . Aici Δx reprezintă distanţa (măsurată de-a lungul liniei centrale a aortei) dintre frontiera de intrare la rădăcina aortică şi frontiera de ieşire la aorta descendentă, iar Δt reprezintă timpul necesar pentru ca unda de debit să parcurgă distanţa dintre cele două frontiere. Timpul Δt este determinat ca interval de timp între punctele de ascensiune ale celor două forme de undă. Locaţia punctului de ascensiune este determinată ca intersecţie dintre curba de ascendentă de debit şi valoarea minimă a debitului (v. fig. 6.3). Curba ascendentă este aproximată cu o linie dreaptă trasată între punctele situate la 20 şi 80% din valoarea maximă a debitului la locaţia respectivă. Odată ce viteza de undă a fost calculată, mărimea E·h/r0, din ecuaţiile (6.11), (6.12), şi (6.14) este calculată astfel:

2

ρ3 2

0

c

r

hE

. (6.15)

Proprietăţile tuturor segmentelor vasculare aparţinând aortei sunt determinate pe baza acestei ecuaţii. Pentru estimarea proprietăţilor vaselor supra-aortice, se foloseşte o abordare uşor modificată, proprietăţile fiind determinate separat pentru fiecare arteră supra-aortică. Se foloseşte această metodă pentru a limita reflexiile obţinute la bifurcaţiile dintre aortă şi aceste artere. Astfel, iniţial se calculează coeficientul de reflexie [Mynard et al., 2008]:

iidp

iidp

YY

YYΓ , (6.16)

46


unde, Yp (Yd) este admitanţa caracteristică a vaselor părinte, respectiv copil. Admitanţa caracteristică reprezintă inversul rezistenţei caracteristice (calculată conform ecuaţiei (6.11)). Sunt în total trei bifurcaţii, una pentru fiecare arteră supra-aortică, iar rezistenţa caracteristică a fiecărui segment supta-aortic este calculată prin setarea coeficientului Γ la valoarea 0:

paortadaortadaortapaortaaorticsupra RRRRR / . (6.17)

Fig. 6.3 Estimarea timpului de tranziţie a debitului între aorta ascendentă şi cea descendentă.

Cunoscându-se rezistenţa caracteristică, E·h/r0 este determinat, după cum urmează:

ρ2

π3 40

2

0

rZ

r

hE aorticsupra . (6.18)

Pentru a se evita viteze de undă cu valori care nu se află în intervalul fiziologic, s-a setat un prag minim de 200 cm/s precum şi un prag maxim de 1200 cm/s pentru vitezele de undă ale arterelor supra-aortice.

Figura 6.4 prezintă o vedere de ansamblu asupra metologiei descrise în cadrul ultimelor două secţiuni, în timp ce tabelul 6.2 prezintă toţi parametrii de intrare necesari, împreună cu modul de achiziţie al acestora. Pentru cazul post-operatoriu, dacă se poate identifica o îngustare reziduală, atunci se va aplica aceeaşi procedură.

6.3 REZULTATE

Metodologia prezentată în secţiunea anterioară a fost validată prin investigarea a cinci pacienţi din cadrul studiului COAST [Ringel et al., 2012] care prezentau coarctaţii înnăscute sau recurente în zona istmului aortei sau în primul segment al aortei descendente.

Suplimentar, studiul include şi următoarele informaţii imagistice: angiograme 3D obţinute prin rezonanţă magnetică după folosirea unei substanţe de contrast (MRA) şi imagini de rezonanţă magnetică 2D cu contrast de fază (PC-MR). Angiogramele vizualizează toracele, inclusiv TAA şi arterele supra-aortice şi permit segmentarea precisă a vaselor de sânge arteriale. Modelul geometric segmentat 3D al arborelui arterial a fost folosit pentru a se determina linia centrală a arterelor precum şi razele (secţiunile transversale) la diverse locaţii. Imaginile PC-MR intersectează, în mod tipic, aorta de două ori. Diversele imagini au fost înregistrate pe baza coordonatelor echipamentului MR, iar după înregistrare, imaginile PC-MR au fost segmentate şi integrate pentru a se determina profilurile personalizate de debit la frontiera de intrare a aortei ascendente şi la frontiera de ieşire a aortei descendente.

Pentru a construi grid-ul geometric discretizat pe baza liniilor centrale ale arterelor şi a secţiunilor transversale, s-a folosit o metodă similară celei introduse în lucrarea [Steele et al.,

47


2003]. Astfel, pentru fiecare segment al modelului arterial, s-au folosit o serie de segmente unidimensionale cu secţiuni transversale variabile în spaţiu, astfel încât să se obţină o geometrie similară geometriei 3D determinată prin MRI. Soluţia numerică la locaţiile interfeţelor dintre segmentele 1D separate a fost determinată prin aplicarea legilor de conservare a debitului şi a impulsului. După revizuirea imaginilor pacientului 5, s-a putut observa o localizare greşită a planului PC (care intersecta ventriculul stâng şi nu aorta ascendentă), care a condus la valori eronate ale debitului de intrare. Astfel, doar patru pacienţi au fost incluşi în procedura finală de evaluare, care este descrisă în continuare.

Tabelul 6.3 prezintă datele specifice pacienţilor. Trei seturi de măsurători bazate pe tensiometru au fost efectuate la braţe şi picioare, iar valorile medii ale acestor măsurători au fost folosite pentru a calcula diferenţa dintre presiunile din partea superioară (braţ) şi inferioară (picior) a corpului la sistolă şi diastolă. Măsurătorile prin ecocardiografie Doppler ale vitezei maxime din amonte şi din aval de coarctaţie au fost folosite pentru a calcula căderea de presiune de-a lungul coarctaţiei pe baza ecuaţiei Bernoulli modificate [Seifert et al., 1999].

Tab. 6.3 Datele achiziţionate de la pacienţi.

Gradientul de presiune [mmHg]

Nr. pac.

SBP [mmHg]

DBP [mmHg]

HR [bpm]

% Îngust. coarctaţie

ascQ [ml/s]

descQ

[ml/s] Cateter Doppler Braţ-

picior

1 124 65 71 59.1 163.9 59.4 55 80 42 2 112 71 61 46.3 97.8 50.2 8 23 -3 3 124 71 118 47.9 88.6 61.5 30 67 40 4 89 50 74 39.5 199.5 85.7 14 27 27

Sângele a fost modelat ca fluid Newtonian incompresibil cu o densitate de 1.055 g/cm3 şi o vâscozitate dinamică de 0.045 dynes/cm2·s în cadrul tuturor celor patru simulări. Spaţiul dintre două noduri consecutive ale grid-ului a fost de 0.05cm, ceea ce a condus la un model de calcul cu 1200-1600 de grade de libertate (valorile variabile în timp şi spaţiu ale ariilor transversale şi ale debitului), valoarea finală depinzând de geometria specifică pacientului. Deoarece s-a adoptat o schemă numerică explicită, pasul de timp este limitat de condiţia CFL, şi a fost fixat la valoarea 2.5e-5s.

Primul pas a fost acela de a evalua cele două abordări diferite pentru calculul rezistenţei introduse de coarctaţie: Rc1 şi Rc2. Tabelul 6.4 prezintă valoarea debitului mediu al aortei descendente determinat prin PC-MRI şi calculat pe baza modelului de ordin redus [Ralovich et al., 2012(c)]. Pentru evaluarea critică a celor două abordări, s-a determinat eroarea medie relativă precum şi eroarea medie absolută a debitului mediu în cadrul aortei descendente. Eroarea relativă a fost calculată pe baza expresiei |Qmasurat – Qcalculat| / | Qmasurat | x 100 iar eroarea absolută pe baza expresiei |Qmasurat – Qcalculat|. Rezultatele sunt prezentate pe ultimele două rânduri ale tabelului 6.4 şi arată faptul că, deşi diferenţele dintre valorile măsurate şi cele calculate sunt mici pentru ambele metode, rezultatele sunt mai precise atunci când se foloseşte Rc2. Prin urmare, în continuare, rezistenţa coarctaţiei va fi calculată pe baza debitului variabil în timp al aortei descendente.

În continuare, se compară, pentru cei patru pacienţi, căderile de presiune calculate pe cale non-invazivă cu ajutorul modelului de ordin redus cu: (a) valorile obţinute conform standardul clinic pe baza măsurătorilor efectuate cu cateter; (b) valorile obţinute prin ecocardiografie Doppler şi ecuaţia modificată a lui Bernoulii; (c) valorile obţinute pe baza diferenţei de presiune dintre partea superioară şi inferioară a corpului. Căderile de presiune determinate pe baza cateterelor şi pe baza modelului de ordin redus au fost calculate ca diferenţă între valoarea maximă în aorta ascendentă şi aorta descendentă. Figura 6.6. prezintă rezultatele acestei

48


49

comparaţii. Rezultatele arată că diferenţele dintre metoda propusă şi măsurătorile invazive sunt foarte

mici, având o eroare medie absolută de 1.45 mmHg şi o eroare medie relativă de 10% pentru ΔP AAo-DAo. Căderile de presiune determinate pe baza Doppler şi valorilor obţinute prin tensiometru au o eroare medie absolută de 23 mmHg, respectiv de 11.75 mmHg, în timp ce, erorile relative medii sunt de 112% şi respectiv de 72%.

Tab. 6.4 Compararea debitului mediu în cadrul aortei descendente şi evaluarea metodelor de calcul a rezistenţei coarctaţiei.

Debit calculat [ml/s] Nr. pacient PC-MRI [ml/s] Folosind Rc1

Folosind Rc2

1 59.39 57.62 58.24 2 50.22 49.79 49.93 3 61.51 61.23 61.40 4 85.69 84.00 84.27 Eroare relativă medie [%] 1.56 1.09 Eroare absolută medie [ml/s] 1.04 0.74

Fig. 6.6 Comparaţie între căderile de presiune dintre aorta ascendentă şi aorta descendentă

(AAo-DAo) pe baza valorilor maxime de presiune, folosind patru metode diferite.

6.4 ANALIZA REZULTATELOR

Rezultatele excelente obţinute în cadrul validării metodei non-invazive de determinare a căderii de presiune de-a lungul coarctaţiilor demonstrează fezabilitatea acesteia pentru o evaluarea clinică precisă. Timpul de calcul a variat între 6 şi 8 minute, ceea ce permite includerea metodei într-un flux de lucru clinic existent. Discrepanţa dintre valorile standard şi cele obţinute prin metodele alternative de evaluare non-invazivă a coarctaţiilor, subliniază suplimentar necesitatea unei metode de evaluare non-invazivă precisă.

Metoda de personalizare a simulărilor prezentată în secţiunea 6.2 nu este aplicabilă doar în cazul modelelor unidimensionale, ci şi în cazul abordărilor tridimensionale. Dacă modelul geometriei utilizat în cadrul simulării este mai detaliat (adică dacă acesta are mai multe ramuri),


metodologia utilizată pentru determinarea condiţiilor de frontieră de ieşire nu se modifică, singura diferenţă fiind aceea că rezistenţele/complianţele totale sunt distribuite către mai multe ramuri, respectiv către mai multe frontiere de ieşire. Suplimentar, dacă razele vaselor terminale scad, distribuirea rezistenţelor poate fi realizată pe baza supoziţiei de tensiune intraparietală constantă [Ibrahim et al., 2010].

50

Deoarece căderea de presiunea de-a lungul coarctaţiei este calculată ca diferenţă dintre valorile maxime obţinute în aorta ascendentă şi aorta descendentă, valoarea finală nu este determinată doar de debitul maxim şi de geometrie, ci de o interacţiune complexă dintre aceste două aspecte, întârzierea de fază introdusă de complianţă [Kadem et al., 2006], viteza de propagare a undei precum şi undele de debit şi de presiune care traversează aorta în sens invers. Deoarece viteza de undă este determinată individual pentru fiecare pacient, metoda propusă este capabilă să modeleze corect momentul sosirii undelor reflectate, care modifică unda de debit şi care pot augmenta presiunea maximă în aorta ascendentă şi descendentă, influenţând astfel rezultatul final al căderii de presiune. Pentru cazul post-operatoriu, în situaţia în care se plasează un stent, acesta va avea, în general, alte proprietăţi elastice decât peretele aortei. Aceasta conduce la o incompatibilitate a impedanţelor la cele două interfeţe cu stentul şi prin urmare la unde reflectate suplimentare la cele două interfeţe, influenţând astfel formele de undă ale debitului şi presiunii. Aceste aspecte motivează alegerea unei metode bazate pe un model unidimensional capabil să surprindă fenomenele de propagare a undelor, în favoarea unei abordări bazate pe un model cu parametrii distribuiţi.

6.5 CONCLUZII

În acest capitol a fost prezentată o abordare bazată pe CFD, cuplată cu o strategie nouă de personalizare pentru evaluarea non-invazivă pre- şi post-operatorie a pacienţilor cu coarctaţie. S-a realizat un studiu de validare pe baza măsurătorilor clinice in-vivo obţinute în cadrul cateterizării de rutină, rezultatele fiind promiţătoare. Abordarea propusă este complet automatizată, nu necesită proceduri iterative de adaptare a parametrilor, precum şi un timp de execuţie total de aproximativ 6-8 minute (obţinându-se astfel o metodologie fezabilă pentru un flux clinic de diagnosticare).

7. UTILIZAREA UNUI PROCESOR GRAFIC PENTRU ACCELERAREA SIMULĂRILOR MULTISCALARE HEMODINAMICE

7.2 PARALELIZAREA ALGORITMULUI NUMERIC DE REZOLVARE A MODELULUI UNIDIMENSIONAL

Anterior, în capitolul 3, au fost prezentate două scheme numerice de ordinul doi pentru rezolvarea modelului unidimensional: Lax-Wendroff (LW) şi dezvoltare în serie Taylor (TS – Taylor Series). De asemenea, au fost prezentaţi un model elastic şi unul vâscoelastic pentru deformarea pereţilor arteriali. În vederea implementării numerice a modelului unidimensional cu o lege vâscoelastică a peretelui arterial, s-a folosit metoda împărţirii operatorilor, care necesită determinarea unui debit de corecţie la fiecare pas. Determinarea debitului de corecţie presupune rezolvarea unui sistem tridiagonal de ecuaţii liniare, care necesită un timp de execuţie relativ mare.

Deoarece metoda LW este formată din doi paşi şi necesită determinarea debitului de corecţie de două ori la fiecare moment de timp discret, în cazul utilizării unei legi vâscoelastice, s-a folosit doar schema numerică TS.

De asemenea, s-au considerat cele două tipuri generale de condiţii de frontieră, aperiodică şi periodică, reprezentate prin două condiţii de frontieră fiziologice: Windkessel (WK) şi respectiv arbore structurat (ST – Structured Tree).


Diversele configuraţii de calcul, pentru care s-au evaluat diferite strategii de paralelizare, sunt prezentate în tabelul 7.1. În continuare, prin algoritm secvenţial exclusiv CPU (SCO – Sequential CPU Only) se va face referire la implementările clasice, bazate pe CPU, a schemelor numerice.

51

Tab. 7.1 Configuraţii de calcul pentru care se evaluează reducerea timpului de execuţie printr-o implementare paralelă.

Caz Schemă numerică Legea peretelui arterial

Condiţie de frontieră la ieşire

1 Lax-Wendroff Elastică Windkessel 2 Lax-Wendroff Elastică Arbore structurat 3 Serie Taylor Elastică Windkessel 4 Serie Taylor Elastică Arbore structurat 5 Serie Taylor Vâscoelastică Windkessel 6 Serie Taylor Vâscoelastică Arbore structurat

Deoarece cele două scheme numerice menţionate anterior sunt explicite, ele pot fi paralelizate în mod eficient. Astfel, au fost testate două abordări diferite: 1. Algoritm paralel hibrid CPU-GPU (PHCG – Parallel Hybrid CPU-GPU): valorile mărimilor

necunoscute în punctele interioare ale grid-ului sunt calculate pe GPU, iar valorile în punctele de la frontierele de intrare/ieşire şi de la bifurcaţii (numite în continuare puncte de joncţiune), sunt calculate pe CPU. Avantajul acestei abordări este că fiecare dispozitiv este folosit pentru tipul de calcule pentru care a fost dezvoltat (CPU – secvenţial, GPU – paralel), dar dezavantajul este că la fiecare moment de timp discret trebuie copiate date între zonele de memorie ale celor două dispozitive pentru a inter-schimba valorile de la sau de lângă punctele de joncţiune;

2. Algoritm paralel exclusiv GPU (PGO – Parallel GPU Only): valorile tuturor mărimilor necunoscute sunt calculate pe GPU, iar CPU-ul este folosit doar pentru a iniţializa zonele de memorie şi pentru a controla execuţia pe GPU. Avantajul este că nu trebuie realizate copieri de date între zonele de memorie ale celor două dispozitive, dar dezavantajul este că operaţii cu un paralelism redus trebuie executate pe GPU.

7.2.1 Algoritmul paralel hibrid CPU-GPU

Iniţial se va face referire la cazul în care se foloseşte un model elastic pentru peretele arterial. Pornind de la schema numerică descrisă în secţiunea 3.3, s-a dezvoltat fluxul de lucru prezentat în figura 7.1 [Itu et al., 2012(c)]. Primul pas constă în iniţializarea modelului arterial (se alocă memorie host pentru fiecare punct al grid-ului: raza iniţială, aria transversală iniţială, elasticitatea peretelui, se calculează derivatele spaţiale ale razei şi ale elasticităţii etc.) şi în alocarea şi iniţializarea memoriei GPU-ului. În continuare, se execută o buclă while, care parcurge momentele de timp discrete până când au fost simulate un anumit număr de cicluri cardiace. În interiorul buclei while, CPU-ul şi GPU-ul procesează datele în paralel până la atingerea unei bariere de sincronizare. În timpul execuţiei paralele, CPU-ul calculează noile valori la punctele de joncţiune iar GPU-ul efectuează calculele corespunzătoare punctelor interioare ale grid-ului (ecuaţiile (3.46) şi (3.47) pentru schema LW şi ecuaţia (3.56) pentru schema TS). Operaţiile realizate de GPU la o iteraţie sunt incluse într-un singur kernel GPU, deoarece nu este necesară o sincronizare la nivel de bloc de fire de execuţie. Deoarece operaţiile efectuate pe GPU sunt asincrone, nu este nevoie de o abordare specială pentru a genera paralelism-ul la nivel de task între CPU şi GPU. Pentru a se putea determina pe CPU valorile mărimilor necunoscute în punctele de joncţiune, trebuie cunoscute valorile de la punctele situate lângă punctele de joncţiune, de la momentul de timp anterior. În mod analog, pentru a se putea


determina pe GPU valorile mărimilor necunoscute în punctele situate lângă punctele de joncţiune, trebuie cunoscute valorile de la punctele de joncţiune, de la momentul de timp anterior. Prin urmare, la începutul şi finalul fiecărei iteraţii, se realizează operaţii de copiere a datelor între zonele de memorie ale celor două dispozitive. Suplimentar, la finalul iteraţiei se introduce o barieră de sincronizare pentru a asigura faptul că valorile de la, respectiv de lângă punctele de joncţiune, au fost anterior inter-schimbate între dispozitivele de calcul. La ultimul ciclu cardiac simulat, după ce s-a atins convergenţa, rezultatele sunt salvate în fişiere pentru a putea fi vizualizate sau post-procesate. Deoarece numărul de iteraţii la fiecare ciclu cardiac este foarte mare (datorită naturii explicite a schemei numerice), rezultate sunt salvate în fişiere numai după un anumit număr de iteraţii (la fiecare 20 până la 50 de iteraţii).

Fig. 7.1 Fluxul de lucru al algoritmului PHCG în cazul în care se foloseşte un model elastic pentru peretele arterial.

Pentru îmbunătăţirea timpului de execuţie necesar pe GPU, a fost optimizat kernel-ul. Scopul a fost acela de a spori intensitatea de calcul şi de a scădea numărul operaţiilor efectuate asupra memoriei globale. Această abordare este necesară pentru a se obţine o performanţă sporită şi atunci când se folosesc arbori arteriali mici, pentru care paralelismul nu este foarte pronunţat.

Configuraţia de execuţie a kernel-ului a fost organizată după cum urmează: fiecare fir de execuţie este responsabil de un singur punct al grid-ului; un bloc de fire de execuţie este creat pentru fiecare segment arterial; şi atât grid-ul de blocuri cât şi grid-ul de fire al unui bloc sunt unidimensionale. Rezolvarea numerică a modelului unidimensional pentru un arbore arterial este practic o metodă de descompunere a domeniului, cu câte un domeniu pentru fiecare segment arterial. Prin urmare, informaţiile trebuie inter-schimbate între segmentele arteriale doar la interfeţele domeniilor. Deoarece pentru algoritmul PHCG punctele de joncţiune sunt rezolvate pe

52


CPU şi, prin urmare, nu există necesitatea comunicării şi sincronizării între operaţiile efectuate pentru punctele interioare ale unor segmente arteriale diferite, asocierea dintre un bloc de fire de execuţie şi un segment arterial a fost o alegere logică. Suplimentar, deoarece paralelismul este redus şi intensitatea de calcul este ridicată, nu a fost luată în considerare o abordare, la care un singur fir de execuţie determină valorile mărimilor necunoscute în mai multe puncte interioare ale grid-ului. Un segment arterial este împărţit în mai multe domenii de calcul numai dacă resursele hardware ale unui multiprocesor sunt insuficiente pentru a rula blocul de fire de execuţie respectiv (în acest caz se introduc puncte de joncţiune suplimentare, iar mărimile necunoscute de la interfeţele acestor domenii sunt determinate prin asigurarea conservării debitului şi a egalităţii presiunii totale). Suplimentar, pentru a se evita accesări redundante ale memoriei globale, realizate de diferite fire de execuţie, valorile de la momentul de timp anterior şi rezultatele intermediare sunt stocate în memoria partajată a multiprocesorului.

Un aspect foarte important pentru algoritmul PHCG este transferul de date dintre CPU şi GPU. Deşi cantitatea de date care trebuie transferate este mică (numai valorile de la sau de lângă punctele de joncţiune trebuie inter-schimbate), timpul total de execuţie necesar pentru realizarea transferului de date este mare. Acest aspect este datorat faptului că numărul total de operaţii de copiere este mare şi faptului că datele care trebuie copiate se află în locaţii izolate ale memoriei globale. Pentru scăderea timpului de execuţie necesar operaţiilor de copiere s-au dezvoltat trei abordări diferite, prezentate în figura 7.2. Acestea au condus la trei variante de implementare a algoritmului PHCG: 1. PHCG cu copieri separate (PHCGCS – PHCG Copy Separately): fiecare locaţie este copiată

separat, obţinându-se în total 8 operaţii de copiere pentru fiecare segment arterial, la fiecare iteraţie. Figura 7.2(a) prezintă locaţiile pentru care valorile sunt interschimbate la fiecare iteraţie precum şi direcţia de copiere. Tablourile de date reprezentate sunt generice, corespunzând fie ariei secţiunilor transversale, fie debitelor unui singur segment arterial;

2. PHCG cu copiere integrală (PHCGCA – PHCG Copy All): tablourile de date utilizate pentru ariile secţiunilor transversale şi pentru debite sunt copiate integral la fiecare iteraţie: 4 operaţii de copiere la fiecare iteraţie (v. fig. 7.2(b));

3. PHCG cu copiere compactă (PHCGCC – PHCG Copy Compact): se alocă tablouri suplimentare pentru locaţiile care trebuie copiate: 4 operaţii de copiere la fiecare iteraţie. Figura 7.2(c) prezintă tablourile suplimentare care trebuie alocate precum şi locaţiile din care/în care se copiază datele. Pentru un arbore arterial, valorile tuturor mărimilor necunoscute de un anumit tip (aria secţiunii transversale sau debit), corespunzătoare locaţiilor care trebuie copiate, sunt stocate într-un singur tablou.

Fig. 7.2 Variante de copiere a datelor între CPU şi GPU: (a) operaţii de copiere separate

pentru fiecare locaţie; (b) copierea integrală a tablourilor; (c) copierea tablourilor suplimentare compacte.

53


Figura 7.4 prezintă fluxul de lucru în cazul utilizării unei legi vâscoelastice pentru peretele arterial [Itu et al., 2013(c)].

Fig. 7.4 Fluxul de lucru al algoritmului PHCG în cazul în care se foloseşte un model vâscoelastic pentru peretele arterial.

Datorită utilizării metodei de împărţire a operatorilor, sunt necesare o serie de modificări în raport cu cazul utilizării unei legi elastice. În primul rând se folosesc două kernel-uri diferite: unul pentru calculul ariilor secţiunilor transversale şi componentei elastice a debitului (ecuaţiile (3.18a) şi (3.60)), şi un al doilea kernel pentru calculul componentei vâscoelastice a debitului (ecuaţia (3.62)). Din nou, CPU-ul şi GPU-ul lucrează în paralel la începutul fiecărei iteraţii. După atingerea unei prime bariere de sincronizare, valorile de la punctele de joncţiune sunt copiate pe GPU pentru a se pregăti calculul componentei vâscoelastice a debitului (în ecuaţia (3.62) se folosesc valorile ariilor secţiunilor transversale şi ale componentei elastice a debitelor de la momentul de timp curent).

Se calculează apoi componenta vâscoelastică a debitului precum şi debitul total, după care valorile de debit de la punctele situate lângă punctele de joncţiune sunt copiate înapoi pe CPU pentru a pregăti efectuarea calculelor de la următoarea iteraţie (debitul total din punctele de joncţiune este dat doar de componenta elastică deoarece se folosesc condiţii de frontieră de tip Dirichlet omogene pentru componenta vâscoelastică a debitului). Pentru rezolvarea sistemului de

54


ecuaţii tridiagonal pe GPU, se foloseşte un algoritm CR-PCR (Cyclic Reduction – Parallel Cyclic Reduction) [Zhang et al., 2010].

7.2.2 Algoritmul paralel exclusiv GPU

Necesitatea execuţiei unor operaţii de copiere la fiecare iteraţie reduce semnificativ performanţa globală a algoritmului PHCG. Prin urmare, s-a realizat o implementare în cadrul căreia toate punctele grid-ului sunt procesate pe GPU. Această abordare elimină operaţiile de copiere între zonele de memorie (sunt necesare doar operaţiile de copiere din cadrul iteraţiilor la care se salvează rezultatele) dar, pe de altă parte, obligă GPU-ul să efectueze calcule cu paralelism redus necesare punctelor de joncţiune. Un alt dezavantaj al algoritmului PGO, comparativ cu algoritmul PHCG, este că, deoarece toate operaţiile sunt efectuate pe GPU, se pierde paralelismul la nivel de task dintre CPU şi GPU. Figura 7.6(a) prezintă fluxul de lucru pentru cazul cel mai complex, şi anume atunci când se foloseşte o lege vâscoelastică pentru peretele arterial împreună cu o condiţie la frontieră de tip arbore structurat. Până la trei kernel-uri trebuie executate la fiecare iteraţie: 1. Calculul integralei de convoluţie (o operaţie de tip scanare înmulţire-adunare [Sengupta et

al., 2007]); 2. Calculul noilor valori ale secţiunilor transversale şi ale componentei elastice a debitului:

ecuaţiile (3.18a) şi (3.60); 3. Calculul componentei vâscoelastice a debitului: ecuaţia (3.62).

Fig. 7.6 Flux de lucru generic al GPU-ului pentru cazul în care se foloseşte o condiţie de

frontieră de tip arbore structurat împreună cu un model vâscoelastic al peretelui arterial. Toate operaţiile sunt efectuate pe GPU, iar CPU-ul realizează doar coordonarea acestora (a). Operaţiile kernel-ului utilizate pentru calculul noilor valori în toate punctele grid-ului (b).

55


În cazul în care se foloseşte condiţia de frontieră WK, atunci primul kernel nu este apelat, iar dacă se foloseşte un model elastic pentru peretele arterial, atunci cel de-al treilea kernel nu este apelat.

Figura 7.6(b) prezintă operaţiile efectuate de cel de-al doilea kernel pentru calculul noilor valori ale secţiunilor transversale şi ale debitelor, accentul fiind pus asupra punctelor de joncţiune. Iniţial, se procesează punctele interioare ale grid-ului după modelul din figura 7.3 (operaţiile nu au mai fost detaliate). În continuare, primul fir de execuţie al primului bloc procesează punctul corespunzător frontierei de intrare, iar ultimul fir de execuţie al fiecărui bloc procesează punctele corespunzătoare frontierelor de ieşire sau ale bifurcaţiilor/conexiunilor (un punct de conexiune este introdus dacă un segment arterial este împărţit în mai multe domenii). Astfel, pentru punctele de joncţiune, paralelismul este prezent doar la nivel de bloc şi nu la nivel de fir de execuţie.

7.3 REZULTATE

Sângele a fost modelat ca fluid incompresibil Newtonian cu o densitate de 1.055 g/cm3 şi o vâscozitate dinamică de 0.045 dynes/cm2s.

Pentru compararea performanţelor celor trei algoritmi (SCO, PHCG, PGO) s-a folosit arborele arterial detaliat prezentat în lucrarea [Stergiopulos et al., 1992], compus din 51 de artere. La frontiera de intrare s-a impus un debit variabil în timp [Olufsen et al., 2000], iar pentru frontierele de ieşire s-au folosit condiţii de frontieră de tip WK şi ST (s-au utilizat valorile parametrilor prezentate în tabelul 3.4). Rezistenţa totală şi complianţa au fost preluate din lucrarea [Bessems, 2007], în timp ce raza minimă utilizată pentru generarea arborilor structuraţi a fost adaptată astfel încât să se obţină o rezistenţă totală similară celei utilizate în cazul condiţiei de frontieră de tip WK (rezistenţă totală: 310371 . sdynes /cm5). Acest aspect a permis compararea variaţiilor în timp ale debitelor şi presiunilor obţinute cu cele două tipuri de condiţii de frontieră.

Pentru frontiera de tip ST s-a folosit un factor exponenţial egal cu 2.7. Constantele care caracterizează asimetria arborelui binar au fost setate la valorile 0.908 şi 0.578, în timp ce raportul dintre lungime şi rază a avut valoarea egală cu 50. Proprietăţile elastice ale peretelui arterial au fost egale atât pentru domeniul unidimensional cât şi pentru domeniul reprezentat de arborele structurat. Împreună cu raza minimă la care arborele structurat este întrerupt, aceşti parametrii determină complianţa arborelui structurat.

Algoritmul secvenţial (SCO) a fost executat pe un procesor Intel i5 CPU cu 2.8GHz, în timp ce pentru algoritmii paraleli (PHCG, PGO) s-a folosit un GPU NVIDIA GTX 460 (336 de core-uri, 48KB memorie partajată şi 32K regiştrii). Toate calculele au fost efectuate cu structuri de date în virgulă mobilă cu precizie dublă, deoarece precizia simplă ar afecta rezultatele mai ales la punctele de joncţiune unde metoda caracteristicilor este aplicată pe baza metodei Newton.

7.3.2 Compararea strategiilor de copiere a datelor pentru algoritmul PHCG

În continuare, s-a evaluat performanţa celor trei strategii de copiere a datelor pentru algoritmul PHCG. Tabelul 7.2 prezintă timpii de execuţie ai operaţiilor efectuate pe GPU, corespunzători unui ciclu cardiac simulat cu pereţi elastici, cu schema LW şi cu condiţia de frontieră de tip WK (s-a ales acest caz ca şi caz reprezentativ). Pentru algoritmul PHCGCS, execuţia kernel-ului ocupă doar 2% din timpul total de execuţie al GPU-ului, performanţa fiind puternic limitată de transferul de date dintre CPU şi GPU. Deşi cantitatea totală de date transferate este mai mare în cazul algoritmului PHCGCA, acesta reprezintă o îmbunătăţire, deoarece numărul operaţiilor de copiere este redus drastic. Cele mai bune rezultate însă se obţin prin intermediul algoritmului PHCGCC, deoarece cantitatea de date de transferat este redusă aşa cum s-a observat în primul caz iar numărul de operaţii de copiere este redus ca şi în al doilea caz. Singurul dezavantaj este acela că anumite fire de execuţie trebuie să introducă date în tablourile suplimentare reprezentate în figura 7.2(c), respectiv să le citească din acestea. Aceasta conduce la o creştere cu 7% a timpului de execuţie total, dar această creştere este mai mult decât

56


compensată de timpul redus necesar pentru operaţiile de copiere. Prin urmare, este considerată numai varianta PHCGCC a algoritmului PHCG în cadrul analizei prezentate în continuare.

Tab. 7.2 Timpii de execuţie [s] ai operaţiilor efectuate pe GPU, obţinuţi pentru simularea unui ciclu cardiac cu cele trei variante ale algoritmului PHCG.

Operaţie PHCGCS PHCGCA PHCGCC

Copiere H→D 32.7 4.62 1.09 Kernel 1.77 1.77 1.89 Copiere D→H 55.8 7.46 1.85

7.3.3 Compararea performanţelor obţinute cu algoritmii SCO, PHCG şi PGO

Tabelul 7.3 prezintă timpii de execuţie obţinuţi pentru cele şase configuraţii de simulare prezentate în tabelul 7.1. Timpii de execuţie corespund unor simulări cu zece cicluri cardiaceşi sunt scoase în evidenţă cele mai bune grade de accelerare a vitezei de execuţie. Spaţiul dintre două puncte ale grid-ului a fost de 0.1 cm iar diferenţa dintre două momente de timp discrete de

s. Valorile sunt preluate din literatura de specialitate şi se bazează pe restricţia CFL. 510555 .Gradul de accelerare variază între 4.31x şi 6.72x. În cazul în care se foloseşte o condiţie la

frontieră de tipul WK, algoritmul PHCGCC are cele mai bune performanţe, iar în cazul în care se foloseşte o condiţie de frontieră de tip ST, algoritmul PGO conduce la cel mai mare grad de accelerare. Timpii de execuţie sunt mai mari în cazul folosirii condiţiei de frontieră de tip ST deoarece trebuie calculată integrala de convoluţie (ecuaţia (3.54)).

Deoarece această operaţie poate fi realizată mai eficient pe GPU decât pe CPU, algoritmul PGO este superior algoritmului PHCGCC atunci când se foloseşte o condiţie de frontieră de tip ST.

Tab. 7.3 Timpii de execuţie şi gradul de accelerare obţinuţi pentru simularea a zece cicluri cardiace cu algoritmii SCO, PHCGCC şi PGO. Primele patru cazuri corespund unui perete arterial elastic cu schemele numerice Lax-Wendroff (LW) sau descompunere în serie Taylor (TS) şi cu condiţii de frontieră de tip Windkessel (WK) sau arbore structurat (ST). Ultimele două cazuri corespund unui perete arterial vâscoelastic cu o schema numerică TS şi cu condiţii la frontiere de tip WK sau ST.

PHCGCC PGO Caz Schemă numerică

Perete arterial

Cond. de frontieră

SCO [s] Timp

[s] Accelerare Timp

[s] Accelerare

1 LW Elastic WK 273.4 63.4 4.31x 121.6 2.25x 2 LW Elastic ST 673.7 464.6 1.45x 135.6 4.97x 3 TS Elastic WK 396.3 64.1 6.18x 132.2 3.00x 4 TS Elastic ST 797.2 463.0 1.62x 146.5 5.44x 5 TS Vâscoel. WK 774.4 115.2 6.72x 249.2 3.11x 6 TS Vâscoel. ST 1179.6 531.7 2.22x 272.8 4.32x

Pentru un perete arterial elastic, în cazul folosirii algoritmului PHCGCC, timpii de execuţie sunt comparabili pentru schemele numerice LW şi TS (pentru ambele tipuri de condiţii de frontieră), cu uşoare avantaje pentru schema LW. În cazul algoritmului SCO, schema LW este superioară schemei TS.

Figura 7.7 prezintă o comparaţie a timpilor de execuţie obţinuţi cu algoritmul SCO şi cu algoritmul paralel care are cea mai bună performanţă pentru fiecare configuraţie de calcul. Cele patru configuraţii de calcul au fost obţinute prin combinarea diferitelor modele de pereţi arteriali

57


precum şi a diferitelor tipuri de condiţii de frontieră de ieşire şi prin alegerea schemei numerice cu cea mai bună performanţă (corespunzător rezultatelor din tab. 7.3).

Fig. 7.7 Timpii de execuţie obţinuţi cu algoritmul SCO şi cu algoritmul paralel care

conduce la cea mai bună performanţă pentru fiecare configuraţie de calcul.

Pentru a analiza efectul parametrilor simulării asupra gradului de reducere a timpului de execuţie, se prezintă în figura 7.8 gradele de accelerare obţinute pentru diferite valori ale spaţiului dintre punctele grid-ului: Δx = 0.25 cm (5486 grade de libertate, 8000 momente discrete de timp pe ciclu); Δx = 0.15 cm (9144 grade de libertate, 12500 momente discrete de timp pe ciclu); Δx = 0.1 cm (13716 grade de libertate, 18000 momente discrete de timp pe ciclu); Δx = 0.05 cm (27432 grade de libertate, 37000 momente discrete de timp pe ciclu).

Fig. 7.8 Gradul de accelerare obţinut pentru diferite

configuraţii ale spaţiului dintre punctele grid-ului.

Numărul momentelor de timp discrete pe ciclu este ales astfel încât să fie satisfăcută condiţia CFL. În fiecare caz, schema numerică şi algoritmul paralel utilizat la procesare corespund celui mai mare grad de accelerare obţinut pentru un spaţiu de 0.1 cm între punctele grid-ului. Figura 7.8 indică o creştere aproximativ liniară a gradului de accelerare, sugerând astfel faptul că GPU-ul nu este folosit la întreaga sa capacitate de calcul pentru niciuna dintre configuraţiile de calcul, pentru un spaţiu mai mare de 0.05 cm între punctele grid-ului.

58


59

7.4 CONCLUZII

Pentru testarea capacităţii de reducere a timpului de execuţie a unei implementări pe GPU a modelului unidimensional, s-a folosit un model arterial complet al circulaţiei umane, cuplat cu două tipuri diferite de condiţii de frontieră alese pe motive fiziologice (model Windkessel cu trei elemente aplicat ca şi condiţie de frontieră aperiodică şi arbore structurat aplicat ca şi condiţie de frontieră periodică) şi cu două modele diferite pentru pereţii arteriali. Gradele de reducere a timpului de execuţie variază între 4.31x şi 6.72x, confirmându-se astfel potenţialul sporit de reducere a timpului de execuţie pentru implementară bazată pe GPU.

Rezultatele au indicat faptul că, dacă se foloseşte o condiţie de frontieră aperiodică, algoritmul PHCGCC conduce la cele mai bune performanţe, în timp ce pentru o condiţie de frontieră periodică, algoritmul PGO conduce la cele mai bune performanţe. Aceste rezultate sunt motivate de mai multe aspecte. În primul rând, algoritmul PGO conduce la reducerea timpului de execuţie nu numai prin paralelizarea procesării pe GPU, dar şi prin paralelismul la nivel de task dintre CPU şi GPU. Deoarece în cazul algoritmului PHCGCC frontierele de ieşire sunt procesate pe CPU, gradul de accelerare obţinut pentru o condiţie de frontieră periodică este limitat. În acest caz, algoritmul PGO are performanţe mai bune, deşi nu foloseşte un paralelism la nivel de task, deoarece operaţia de scanare înmulţire-adunare necesară pentru calculul integralei de convoluţie este executată mai eficient pe GPU (mai ales atunci când numărul momentelor de timp discrete este mare, aşa cum este în cazul de faţă).

8. CONCLUZII FINALE

8.1 CONCLUZII

Cercetările doctorale au avut ca principal scop dezvoltarea, implementarea, testarea şi validarea unui model multiscalar al circulaţiei coronariene pentru diagnosticarea non-invazivă a stenozelor coronariene. Deoarece modelul multiscalar dezvoltat are un caracter general, obiectivul secundar al tezei a fost acela de a dezvolta, implementa, testa şi valida un model multiscalar pentru diagnosticarea non-invazivă a coarctaţiilor aortice. Având în vedere că aceste modele au fost dezvoltate pentru a fi utilizate în clinici de specialitate pentru diagnosticarea sistemului cardiovascular al pacienţilor care suferă de aceste tipuri de patologii, un timp de execuţie redus este crucial. În acest sens, un alt obiectiv secundar a fost acela de a reduce timpul de execuţie al modelului cvasi unidimensional, care reprezintă componenta de bază a modelelor multiscalare dezvoltate. Prin urmare, algoritmii corespunzători metodelor numerice de rezolvare a ecuaţiilor diferenţiale care stau la baza modelelor utilizate, au fost paralelizaţi în vederea execuţiei lor pe un procesor grafic. Implementările accelerate ale modelelor multiscalare sunt necesare şi atunci când modelele multiscalare trebuie să furnizeze rezultate în timp real în timpul intervenţiilor chirurgicale.

După cum a fost menţionat şi anterior, afecţiunile cardiovasculare reprezintă principală cauză de deces la nivel mondial. Dintre multitudinea de boli cardiovasculare, afecţiunile circulaţiei coronariene reprezintă cauza a aproape 50% dintre decese [***WHO, 2011]. Prin urmare, dezvoltarea unor metode non-invazive de diagnosticare a acestor afecţiuni este de interes maxim, atât din punctul de vedere al reducerii riscurilor asociate metodelor invazive de diagnosticare, cât şi din punctul de vedere al reducerii costurilor serviciilor medicale.

Noţiunile fundamentale privind simularea numerică multiscalară a hemodinamicii arteriale au fost prezentate în prima parte lucrării. Această abordare s-a impus în literatura de specialitate deoarece, pe de o parte, sistemul cardiovascular are un caracter de circuit închis şi, ca urmare, trebuie modelate toate componentele sale esenţiale iar, pe de altă, parte resursele necesare pentru simulările tridimensionale sunt prea mari pentru ca întregul circuit să fie modelat la acelaşi nivel de detaliu. Principalele tipuri de modele utilizate sunt modelele tridimensionale, modelele cvasi unidimensionale şi modelele cu parametrii distribuiţi. Deoarece condiţiile de frontieră joacă un rol foarte important în cadrul simulărilor cardiovasculare, s-au prezentat principalele modele


utilizate pentru impunerea acestor condiţii (modelul Windkessel cu trei elemente şi arborele structurat). Deoarece, timpul de execuţie este considerabil chiar şi în cazul modelelor multiscalare, s-a realizat şi o scurtă introducere în tehnicile de procesare paralelă prin utilizarea unui procesor grafic, care sunt utilizate pentru reducerea timpului de execuţie.

60

Activităţile de cercetare desfăşurate de-a lungul studiilor doctorale pot fi organizate pe trei mari direcţii: modele/metodologii generale care pot fi aplicate atât în cazul geometriilor arteriale

generice cât şi în cazul geometriilor arteriale achiziţionate de la pacienţi; modele/metodologii destinate în special geometriilor achiziţionate de la pacienţi; algoritmi pentru reducerea timpului de execuţie al modelelor utilizate pentru simulare

sistemului cardiovascular. Componenta de bază a metodologiilor prezentate în cursul tezei o reprezintă modelul cvasi

unidimensional. Acesta este dedus pe baza ecuaţiilor Navier-Stokes tridimensionale prin intermediul unor ipoteze simplificatoare. În literatura de specialitate s-a demonstrat că acest model poate fi utilizat pentru estimarea presiunilor şi debitelor arteriale în geometrii arteriale sănătoase. În acest sens, în capitolul 3 s-au prezentat rezultatele simulărilor hemodinamice obţinute pentru o geometrie arterială generică. S-au utilizat diferite condiţii de frontieră de ieşire (modelul Windkessel şi arborele structurat) precum şi diferite modele ale peretelui arterial (un model elastic şi două modele vâscoelastice). În vederea implementării corecte a modelelor vâscoelastice, atât în interiorul segmentelor vasculare cât şi la joncţiunea cu alte segmente, s-a dezvoltat un nou algoritm de cuplare. Rezultatele au arătat că utilizarea unui model vâscoelastic introduce diferenţe importante în comparaţie cu utilizarea unui model elastic pentru peretele arterial dar, pentru a se obţine rezultatele corecte, este necesar ca modelul vâscoelastic să fie implementat corect în toate punctele grid-ului.

Tot în prima direcţie de cercetare se încadrează şi aspectele prezentate în capitolul 4. S-a introdus o procedură evoluată de optimizare a simulărilor hemodinamice. Scopul principal al acestei proceduri este ca în cadrul simulărilor hemodinamice să se obţină anumite valori impuse de presiune şi debit. În comparaţie cu algoritmii introduşi anterior în literatura de specialitate, metodologia propusă în această lucrare permite reducerea numărului de iteraţii necesare pentru finalizarea procedurii de optimizare şi reduce riscul de eşuare al acesteia. Ambele aspecte sunt importante atunci când se are în vedere utilizarea procedurii într-un mediu clinic, în vederea realizării simulărilor hemodinamice pentru o geometrie specifică unui pacient. Trebuie subliniat faptul că procedura evoluată de optimizare nu a fost testată doar pentru geometrii arteriale sănătoase ci şi pentru geometrii arteriale patologice (geometrie cu coarctaţie aortică).

În continuare, atenţia a fost îndreptat spre simulările hemodinamice specifice pacienţilor (cea de-a doua direcţie de cercetare). Astfel, s-a introdus un model multiscalar al circulaţiei coronariene, utilizat pentru diagnosticarea non-invazivă a stenozelor coronariene. Acest model permite determinarea valorii indicatorului diagnostic FFR pe cale non-invazivă. În vederea aplicării modelului pentru geometriile coronariene achiziţionate de la pacient, s-au desfăşurat diferite activităţi de personalizare a modelelor. Astfel, s-a personalizat modelul inimii utilizat la frontiera de intrare a modelului multiscalar şi s-au personalizat modelele microvasculare utilizate la frontierele de ieşire ale vaselor arborelui coronarian, atât pentru starea de repaus cât şi pentru starea de hiperemie. Suplimentar, deoarece indicatorul diagnostic FFR se determină pe baza presiunilor arteriale, s-a dezvoltat o procedură de optimizare specifică circulaţiei coronariene, care permite obţinerea în cadrul simulării hemodinamice a aceloraşi presiuni sistolice, diastolice şi medii ca şi cele achiziţionate de la pacienţi. În vederea reducerii timpului de execuţie, care poate deveni substanţial în cazul în care această procedură de optimizare necesită multe iteraţii, procedura este iniţial aplicată unui model 0D detaliat, apoi unui model 1D cu geometrie redusă şi doar în final geometriei 1D complete. Modelul multiscalar a fost validat cu ajutorul a patru pacienţi, diferenţa dintre valoarea FFR măsurată şi cea simulată fiind mai mică sau egală cu 0.05, clasificarea stenozelor fiind corectă pentru toţi cei patru pacienţi.

Tot în cea de-a doua direcţie de cercetare se încadrează şi modelul multiscalar propus


pentru diagnosticarea non-invazivă a coarctaţiilor aortice. Şi în acest caz, diagnosticarea se realizează pe baza presiunilor aortice, şi anume folosind căderea de presiune maximă de-a lungul coarctaţiei în timpul unui ciclu cardiac. Componentele esenţiale ale modelului propus sunt reprezentate de modelul cvasi unidimensional, un model comprehensiv de cădere de presiune şi o metodologie de personalizare pentru modelul peretelui arterial şi pentru modelele Windkessel utilizate la frontierele de ieşire ale geometriei aortice. Validarea modelului propus a fost realizată cu ajutorul a patru pacienţi, eroarea absolută medie dintre căderea de presiune măsurată cu un cateter şi căderea de presiune simulată fiind de 1.45 mmHg.

61

De asemenea, în anexa 1 se prezintă o procedură automată pentru adaptarea condiţiei de frontieră de tip arbore structurat. Procedura se concentrează asupra rezistenţei şi complianţei acestui tip de condiţie de frontieră, valorile celor două proprietăţi fiind adaptate prin intermediul parametrilor arborelui structurat. Suplimentar, în anexa 2 se prezintă o metodologie pentru estimarea presiunilor aortice în cazul geometriilor sănătoase, prin utilizarea unor condiţii de frontieră de tip arbore structurat. În final, atenţia a fost îndreptată spre reducerea timpilor de execuţie ai modelelor utilizate în simulările hemodinamice prin utilizarea procesării paralele pe un procesor grafic. Deoarece, după cum a fost precizat, modelul cvasi unidimensional reprezintă componenta principală a modelelor multiscalare dezvoltate în cadrul tezei, la finalul tezei se evaluează potenţialul de reducere a timpului de execuţie prin execuţia modelului pe un procesor grafic GTX460. S-au dezvoltat doi algoritmi, unul hibrid CPU-GPU şi unul exclusiv GPU, aceştia fiind evaluaţi pentru diferite configuraţii (două scheme numerice diferite, condiţii de frontieră periodice şi aperiodice, modele elastice şi vâscoelastice pentru peretele arterial). S-a observat astfel că algoritmul hibrid conduce la cele mai bune performanţe în cazul folosirii unui model elastic pentru peretele arterial, în timp ce, algoritmul exclusiv GPU conduce la cele mai bune performanţe în cazul folosirii unui model vâscoelastic pentru peretele arterial.

De asemenea, în anexa 3 se prezintă implementarea optimizată pe GPU a algoritmilor bazaţi pe şabloane. Simulările au fost efectuate pe baza unei probleme de transfer de căldură staţionar în două dimensiuni iar schema de calcul folosită la rezolvarea numerică a ecuaţiei este larg răspândită în cadrul algoritmilor paraleli de rezolvare a ecuaţiilor diferenţiale cu derivate parţiale de acest tip (schema PSOR roşu-negru). În final, în anexa 4 se prezintă implementarea optimizată pe GPU a ecuaţiilor Navier-Stokes. Pentru rezolvarea numerică a acestor ecuaţii s-a folosit metoda compresibilităţii artificiale, calculele fiind efectuate în două kernel-uri separate, unul pentru actualizarea valorilor celor două componente ale vitezei şi unul pentru calculul noilor valori ale presiunii.

8.2 CONTRIBUŢII ORIGINALE

Prezenta lucrare de doctorat se înscrie în sfera preocupărilor de interes general din domeniul simulărilor hemodinamice multiscalare. Subiectele tratate reprezintă un demers ştiinţific complex şi de strictă actualitate, cu un caracter interdisciplinar important. Abordarea direcţiilor de cercetare urmate a presupus utilizarea, în cadrul proiectului, a unor cunoştinţe aprofundate din domenii ştiinţifice diverse precum matematica, dinamica fluidelor, teoria sistemelor şi, nu în ultimul rând, tehnologia informaţiilor. În acest subcapitol sunt evidenţiate contribuţiile originale pe care teza de doctorat le propune în domeniul simulărilor hemodinamice multiscalare. Contribuţiile aduse se încadrează în următoarele arii de cercetare: modele cardiovasculare; diagnosticarea non-invazivă a patologiilor cardiovasculare; proceduri de optimizare pentru personalizarea simulărilor hemodinamice; accelerarea simulărilor hemodinamice prin utilizarea unui procesor grafic.

8.2.1 Modele cardiovasculare

1. Dezvoltarea unui model multiscalar de ordin redus al circulaţiei coronariene. Acest model permite determinarea caracteristicilor hemodinamice principale, presiune şi debit, în circulaţia


coronariană atât în starea de repaus cât şi în starea de hiperemie intracoronariană. Una dintre caracteristicile principale ale modelului propus este cuplarea la frontiera de intrare a modelului cvasi unidimensional a unui model simplificat al inimii, care permite evaluarea presiunii intraventriculare. Această presiune este, la rândul ei, folosită în cadrul modelelor microvasculare ale circulaţiei coronariene pentru a modela efectul contracţiilor cardiace asupra circulaţiei coronariene.

62

2. Dezvoltarea unui algoritm pentru implementarea corespunzătoare a modelului vâscoelastic la frontierele de intrare şi de ieşire ale segmentelor vasculare. Algoritmul propus prezintă avantajul că poate fi folosit nu numai pentru cuplarea unor domenii 1D ci şi a unor domenii 3D. Acesta a fost testat pe baza unui arbore arterial complet, demonstrându-se că vâscoelasticitatea conduce la modificări de până la 6% pentru presiune şi până la 20% pentru debit. Rezultatele demonstrează că este important să se considere şi să se implementeze componenta vâscoasă a relaţiei arie-presiune la toate punctele grid-ului, inclusiv la punctele de cuplare dintre vase sau la frontierele de intrare/ieşire ale modelului. 3. Dezvoltarea unui model pentru estimarea presiunilor aortice pe baza datelor achiziţionate prin rezonanţă magnetică. Modelul detaliat în anexa 2 prezintă următoarele particularităţi: foloseşte condiţii de frontieră de tip arbore structurat, care reduce numărul parametrilor de estimat în comparaţie cu condiţiile de frontieră de tip Windkessel. Suplimentar, pentru estimarea proprietăţilor pereţilor arteriali se foloseşte o metodă bazată pe minimizarea erorii pătratice. Rezultatele demonstrează că modelul propus permite estimarea corectă nu doar a presiunilor ci şi a ariilor secţiunilor transversale de-a lungul unui ciclu cardiac. 4. Dezvoltarea unui algoritm pentru modelarea autoreglării coronariene. În vederea simulării corecte a circulaţiei coronariene în starea de repaus, este necesar să se modeleze corect fenomenul de autoreglare, care asigură menţinerea constantă a debitului coronarian chiar şi în cazul unor stenoze coronariene pronunţate. În acest sens s-a dezvoltat un algoritm recursiv, care permite adaptarea de la un ciclu cardiac la altul a rezistenţelor modelelor circulaţiei microvasculare astfel încât debitul să fie menţinut la un nivel constant.

8.2.2 Diagnosticarea non-invazivă a patologiilor cardiovasculare

5. Dezvoltarea, implementarea şi validarea unei metodologii pentru diagnosticarea non-invazivă a stenozelor coronariene. Indicatorul diagnostic principal care permite evaluarea funcţională a severităţii unei stenoze este indicatorul FFR. În vederea evaluării valorii acestui indicator pentru geometriile coronariane achiziţionate de la pacienţi, este nevoie ca modelul multiscalar de ordin redus al circulaţiei coronariene să fie personalizat. Activităţile principale de personalizare ale metodologiei propuse sunt personalizarea modelului simplificat al inimii, precum şi personalizarea modelelor circulaţiei coronariene microvasculare pentru stările de repaus şi de hiperemie ale pacientului. Metodologia propusă a fost validată pe baza datelor achiziţionate de la patru pacienţi, rezultatele obţinute fiind foarte bune. 6. Dezvoltarea, implementarea şi validarea unei metodologii pentru diagnosticarea non-invazivă a coarctaţiei aortice. Indicatorul diagnostic principal care permite evaluarea funcţională a severităţii unei coarctaţii aortice este căderea de presiune maximă de-a lungul coarctaţiei, în acest sens, fiind introdus un model comprehensiv de cădere de presiune. În vederea evaluării indicatorului diagnostic pentru geometria aortică achiziţionată de la pacienţi, este nevoie ca modelul cvasi unidimensional al circulaţiei aortice să fie personalizat. Activităţile principale de personalizare ale metodologiei propuse sunt personalizarea parametrilor modelului peretelui arterial precum şi a parametrilor modelelor Windkessel utilizate la cele patru frontiere de ieşire ale geometriei aortice. Metodologia propusă a fost validată pe baza datelor achiziţionate de la patru pacienţi, rezultatele obţinute fiind foarte bune.

8.2.3 Proceduri de optimizare pentru personalizarea simulărilor hemodinamice

7. Dezvoltarea, implementarea şi testarea unei proceduri de optimizare evoluate pentru personalizarea simulărilor hemodinamice. Principalele caracteristici ale procedurii propuse sunt


estimarea rezistenţelor şi complianţelor modelului multiscalar, precum şi compensarea adecvată a acestor mărimi atunci când se trece de la modelul cu parametrii distribuiţi la modelul multiscalar. Pentru compensarea celor două mărimi au fost dezvoltaţi algoritmi specializaţi. Avantajele principale ale procedurii propuse sunt posibilitatea aplicării acesteia în situaţii patologii (coarctaţie aortică), reducerea riscului de eşuare a procedurii de optimizare, precum şi reducerea timpului total de execuţie prin reducerea numărului de iteraţii. Procedura propusă a fost testată cu ajutorul a cinci configuraţii arteriale, rezultatele fiind foarte bune.

63

8. Dezvoltarea, implementarea şi testarea unei proceduri de optimizare pentru personalizarea simulărilor hemodinamicec coronariene. În vederea obţinerii în cadrul simulărilor hemodinamice ale circulaţiei arteriale coronariene a aceloraşi caracteristici hemodinamice (presiune, debit) ca şi cele ale pacientului s-a dezvoltat o procedură de optimizare complexă, care constă din mai multe etape: iniţial, procedura de optimizare este aplicată pentru un model detaliat cu parametrii distribuiţi, apoi pentru un model 1D cu geometrie redusă şi, în final, pentru modelul 1D cu geometrie complexă. În acest fel, numărul de iteraţii precum şi timpul total de execuţie necesar pentru modelul cu geometrie 1D completă este redus, în timp ce, presiuniile medii, sistolice şi diastolice corespund celor măsurate/estimate la pacient cu o eroare mai mică de 0.1 mmHg. 9. Dezvoltarea, implementarea şi testarea unei proceduri de optimizare pentru condiţia de frontieră de tip arbore structurat. Procedura prezentată în anexa 1 permite adaptarea parametrilor condiţiei de frontieră de tip arbore structurat în vederea obţinerii unor valori impuse de rezistenţă şi complianţă totală. În acest fel, condiţia de frontieră de tip arbore structurat poate fi folosită pentru simulări hemodinamice specifice pacienţilor.

8.2.4 Accelerarea simulărilor hemodinamice prin utilizarea unui procesor grafic

10. Dezvoltarea, implementarea şi testarea unor algoritmi pentru accelerarea timpului de execuţie al modelului cvasi unidimensional prin utilizarea unui procesor grafic. S-au dezvoltat doi algoritmi, un algoritm hibrid CPU-GPU (PHCG) şi un algoritm exclusiv GPU (PGO), care sunt comparaţi cu algoritmul secvenţial bazat pe CPU (SCO). Gradele de reducere a timpului de execuţie variază între 4.31x şi 6.72x, confirmându-se astfel potenţialul sporit de reducere a timpului de execuţie pentru implementarea bazată pe GPU 11. Dezvoltarea, implementarea şi testarea unui algoritm pentru reducerea timpului de execuţie a algoritmilor bazaţi pe şabloane prin utilizarea unui procesor grafic. Simulările au fost efectuate pe baza unei probleme de transfer de căldură staţionar în două dimensiuni, calculele efectuate pe GPU fiind împărţite în două kernel-uri. Abordarea optimă, prezentată în anexa 3, nu poate fi aplicată doar în cadrul rezolvării numerice a ecuaţiilor eliptice sau a ecuaţiilor cu derivate parţiale, ci şi în cazul oricărui kernel GPU ale cărui blocuri de fire trebuie să acceseze valorile învecinate ale unui domeniu de calcul 2D, sau 3D. 12. Dezvoltarea, implementarea şi testarea unui algoritm pentru reducerea timpului de execuţie a soluţiei numerice a ecuaţiilor Navier-Stokes prin utilizarea unui procesor grafic. Calculele au fost efectuate în două kernel-uri separate, unul pentru actualizarea valorilor celor două componente ale vitezei şi unul pentru calculul noilor valori ale presiunii. Versiunile optimizate ale kernel-urilor au fost folosite pentru a se realiza o comparaţie a algoritmului hibrid CPU-GPU cu algoritmul clasic bazat pe CPU. Comparaţiile au fost efectuate pe trei grid-uri cu densitate diferită, gradul de îmbunătăţire al performanţelor variind de la sub un ordin de mărime până la aproximativ două ordine de mărime.

8.3 DISEMINAREA REZULTATELOR

Cercetările întreprinse în perioada de studii doctorale au permis elaborarea şi publicarea, în calitate de prim autor sau coautor, a unui număr de 20 articole ştiinţifice în reviste şi buletine de specialitate, respectiv în volumele unor conferinţe ştiinţifice internaţionale, după cum urmează: articole publicate în reviste de specialitate, indexate ISI cu factor de impact: 3; articole publicate în reviste de specialitate, indexate în baze de date internaţionale: 3;


64

articole publicate în volumele unor manifestări ştiinţifice internaţionale, indexate în baze de date internaţionale: 11;

articole publicate în volumele unor manifestări ştiinţifice internaţionale: 3.

Reviste indexate ISI cu factor de impact

Itu, L. M., Sharma, P., Gulsun, M. A., Mihalef, V., Kamen, A., Greiser, A., “Determination of Time-varying Pressure Field from Phase Contrast MRI Data”, Journal of Cardiovascular Magnetic Resonance, Society of Cardiovascular Magnetic Resonance, USA, Vol. 14, Suppl. 1, pp. 36, 2012.

Itu, L. M., 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, Vol. 41, No. 4, pp. 669-681, 2013.

Itu, L. M., Sharma, P., Kamen, A., D., Suciu, C., Comaniciu, D., “Graphics Processing Unit Accelerated One-Dimensional Blood Flow Computation in the Human Arterial Tree”, International Journal on Numerical Methods in Biomedical Engineering, published through early view, doi: 10.1002/cnm.2585.

Reviste indexate BDI

Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A. “Comparison of Single and Double Floating Point Precision Performance for Tesla Architecture GPUs”, Bulletin of the Transilvania University of Brasov - Series I, Engineering Sciences, Vol. 53, pp. 131-138, 2011.

Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “GPU Enhanced Stream-Based Matrix Multiplication”, Bulletin of the Transilvania University of Brasov - Series I Engineering Sciences, Vol. 54, pp. 79-86, 2012.

Niţă, C., Itu, L. M., Suciu, C. “GPU Accelerated Fluid Flow Computations using the Lattice Boltzmann Method” Bulletin of the Transilvania University of Brasov - Series I, Engineering Sciences, Vol. 55, pp. 67-74, 2013.

Conferinţe internaţionale cu volume indexate BDI

Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A. “Optimized GPU Based Simulation of the Incompressible Navier-Stokes Equations on a MAC Grid”, Proc. of the 10th RoEduNet Inter. Conf. - RoEduNet, Iaşi, Romania, June 2011, pp. 5-8.

Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A. “GPU Optimized Computation of Stencil Based Algorithms”, Proc. of the 10th RoEduNet Inter. Conf. - RoEduNet, Iaşi, Romania, June 2011, pp. 1-4.

Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A. “GPU Accelerated Simulation of Elliptic Partial Differential Equations”, Proc. of the 6th IEEE Inter. Conf. on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications - IDAACS, Prague, Czech Republic, Sept. 2011, pp. 238-242.

Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A. “GPU Optimized Computation of the Artificial Compressibility Method”, Proc. of the 15th Inter. Conf. on System Theory, Control and Computing - ICSTCC, Sinaia, Romania, Oct. 2011, pp. 1-6.

Itu, L. M., Suciu, C., Postelnicu, A., Moldoveanu, F. “Analysis of Outflow Boundary Condition Implementations for 1D Blood Flow Models”, Proc. of IEEE Inter. Conf. on e-Health and Bioengineering - EHB, Iasi, Romania, Nov. 2011, pp. 1-4.

Itu, L. M., Sharma, P., Mihalef, V., Kamen, A., Suciu, C., Comaniciu, D. “A Patient-


specific Reduced-order Model for Coronary Circulation”, Proc. of the IEEE Inter. Symp. On Biomedical Imaging - ISBI, Barcelona, Spain, May 2012, pp. 832-835.

65

Itu, L. M., Sharma, P., Kamen, A., Suciu, C., Moldoveanu, F., Postelnicu, A. “GPU Accelerated Simulation of the Human Arterial Circulation”, Proc. of the 13th Inter. Conf. on Optimization of Electrical and Electronic Equipment - OPTIM, Brasov, Romania, May 2012, pp. 1478-1485.

Sharma, P., Itu, L. M., Zheng, X., Kamen, A., Bernhardt, D., Suciu, C., Comaniciu, D. “A Framework for Personalization of Coronary Flow Computations During Rest and Hyperemia”, Proc. of the 34th Annual Inter. Conf. of the IEEE Engineering in Medicine & Biology Society - EMBC, San Diego, California, USA, Aug. 2012, pp. 6665 - 6668.

Ralovich, K., Itu, L. M., Mihalef, V., Sharma, P., Ionasec, R., Vitanovski, D., Krawtschuk, W., Everett, A., Ringel, R., Navab, N., Comaniciu, D. “Hemodynamic Assessment of Pre- and Post-Operative Aortic Coarctation from MRI”, Proc. of the 15th Inter. on Medical Image Computing and Computer Assisted Intervention - MICCAI, Nice, France, Oct. 1-5, 2012, published in “Lecture Notes in Computer Science”, Ed. Springer, Vol. 7511, 2012, pp. 486-493.

Itu, L. M., Sharma, P., Kamen, A., D., Suciu, C., Comaniciu, D. “A Novel Coupling Algorithm for Computing Blood Flow in Viscoelastic Arterial Models”, Proc. of the 35th Annual Inter. Conf. of the IEEE Engineering in Medicine & Biology Society - EMBC 2013, Osaka, Japan, Jul. 2013.

Niţă, C., Itu, L. M., Suciu, C. “GPU Accelerated Blood Flow Computation using the Lattice Boltzmann Method”, 17th IEEE High Performance Extreme Computing Conference, Waltham, MA, USA, Sept. 2013.

Conferinţe internaţionale

Ralovich, K., Mihalef, V., Sharma, P., Itu, L. M., Vitanovski, D., Ionasec, R., Suehling, M., Everett, A., Pongiglione, G., Comaniciu, D., Navab, N. “Modeling and Simulation Framework for Hemodynamic Assessment of Aortic Coarctation Patients”, Proc. of the 20th Annual ISMRM Meeting & Exhibition – ISMRM, Melbourne, Australia, May 2012.

Itu, L. M., Sharma, P., Zheng, X., Mihalef, V., Kamen, A., Suciu, C. “Patient-Specific Modeling and Hemodyanmic Simulation in Healthy and Diseased Coronary Arteries”, Proc. of the ASME 2012 Summer Bioengineering Conference - SBC, Fajardo, Puerto Rico, June 2012.

Ralovich, K., Itu, L. M., Mihalef, V., Sharma, P., Ionasec, R., Vitanovski, D., Everett, A., Krawtschuk, W., Suehling, M., Navab, N., Comaniciu, D. “Non-invasive Assessment of Aortic Coarctation through Blood Flow Computation and MRI”, Proc. of the Virtual Physiological Human - VPH, London, UK, Sept. 2012.

BIBLIOGRAFIE SELECTIVĂ [Alastruey et al., 2011] Alastruey, J., Khir, A., Matthys, K., Segers, P., Sherwin, S., Verdonck, P., Parker, K., Peiro, J. “Pulse Wave Propagation in a Model Human Arterial Network: Assessment of 1-D Visco-elastic Simulations Against in Vitro Measurements”, Journal of Biomechanics, Vol. 44, pp. 2250–2258, 2011.

[Bessems, 2007] Bessems, D. “On the Propagation of Pressure and Flow Waves through the Patient-Specific Arterial System”, Ph.D. Thesis, TU Eindhoven, Netherlands, 2007.

[Chorin, 1967] Chorin, A. J. “A Numerical Method for Solving Incompressible Viscous Problems”, Journal of Computational Physics, Vol. 2, pp. 12–26, 1967.

http://link.springer.com/bookseries/558


[Coogan et al., 2011] Coogan, J. S., Chan, F. P., Taylor, C. A., Feinstein, J. A. “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, Vol. 77, pp. 680–691, 2011.

66

[Figueroa et al., 2006] Figueroa, C. A., Vignon-Clementel, I. E., Jansen, K. C., Hughes, T. J. R., Taylor, C. A. “A Coupled Momentum Method for Modeling Blood Flow in Three-Dimensional Deformable Arteries”, Computer Methods in Applied Mechanics and Engineering, Vol. 195, pp. 5685–5706, 2006.

[Formaggia et al., 2003] Formaggia, L., Lamponi, D., Quarteroni, A. “One Dimensional Models for Blood Flow in Arteries”, Journal of Engineering Mathematics, Vol. 47, pp. 251–276, 2003.

[Huo et al., 2012] Huo, Y., Svendsen, M., Choy, J. S., Zhang, Z. D., Kassab, G. S. “A Validated Predictive Model of Coronary Fractional Flow Reserve”, Journal of the Royal Society Interface, Vol. 9, pp. 1325–38, 2012.

[Ismail et al., 2013] Ismail, M., Wall, W. A., Gee, M. W. “Adjoint-based Inverse Analysis of Windkessel Parameters for Patient-Specific Vascular Models”, Journal of Computational Physics, Vol. 244, pp. 113–130, 2013.

[Itu et al., 2011(a)] Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “Optimized GPU Based Simulation of the Incompressible Navier-Stokes Equations on a MAC Grid”, Proc. of the 10th RoEduNet Inter. Conf. - RoEduNet, Iaşi, Romania, June 2011, pp. 5–8.

[Itu et al., 2011(b)] Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “GPU Optimized Computation of Stencil Based Algorithms”, Proc. of the 10th RoEduNet Inter. Conf. - RoEduNet, Iaşi, Romania, June 2011, pp. 1–4.

[Itu et al., 2011(c)] Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “GPU Accelerated Simulation of Elliptic Partial Differential Equations”, Proc. of the 6th IEEE Inter. Conf. on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications - IDAACS, Prague, Czech Republic, Sept. 2011, pp. 238–242.

[Itu et al., 2011(d)] Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “GPU Optimized Computation of the Artificial Compressibility Method”, Proc. of the 15th Inter. Conf. on System Theory, Control and Computing - ICSTCC, Sinaia, Romania, Oct. 2011, pp. 1–6.

[Itu et al., 2011(e)] Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “Comparison of Single and Double Floating Point Precision Performance for Tesla Architecture GPUs”, Bulletin of the Transilvania University of Brasov - Series I Engineering Sciences, Vol. 53, pp. 131–138, 2011.

[Itu et al., 2011(f)] Itu, L. M., Suciu, C., Postelnicu, A., Moldoveanu, F., “Analysis of Outflow Boundary Condition Implementations for 1D Blood Flow Models”, Proc. of IEEE Inter. Conf. on e-Health and Bioengineering - EHB, Iasi, Romania, Nov. 2011, pp. 1–4.

[Itu et al., 2012(a)] Itu, L. M., Sharma, P., Gulsun, M. A., Mihalef, V., Kamen, A., Greiser, A., “Determination of Time-varying Pressure Field from Phase Contrast MRI Data”, Journal of Cardiovascular Magnetic Resonance, Society of Cardiovascular Magnetic Resonance, USA, Vol. 14, Suppl. 1, pp. 36, 2012.

[Itu et al., 2012(b)] Itu, L. M., Sharma, P., Mihalef, V., Kamen, A., Suciu, C., Comaniciu, D., “A Patient-specific Reduced-order Model for Coronary Circulation”, Proc. of the IEEE Inter. Symp. On Biomedical Imaging - ISBI, Barcelona, Spain, May 2012, pp. 832–835.

[Itu et al., 2012(c)] Itu, L. M., Sharma, P., Kamen, A., Suciu, C., Moldoveanu, F., Postelnicu, A., “GPU Accelerated Simulation of the Human Arterial Circulation”, Proc. of the 13th Inter. Conf. on Optimization of Electrical and Electronic Equipment - OPTIM, Brasov, Romania, May 2012, pp. 1478–1485.

[Itu et al., 2012(d)] Itu, L. M., Sharma, P., Zheng, X., Mihalef, V., Kamen, A., Suciu, C., “Patient-specific Modeling and Hemodyanmic Simulation in Healthy and Diseased Coronary Arteries”, Proc. of the ASME 2012 Summer Bioengineering Conference - SBC, Fajardo, Puerto Rico, June 2012.

[Itu et al., 2012(e)] Itu, L. M., Suciu, C., Moldoveanu, F., Postelnicu, A., “GPU Enhanced Stream-Based


Matrix Multiplication”, Bulletin of the Transilvania University of Brasov - Series I Engineering Sciences, Vol. 54, pp. 79–86, 2012.

67

[Itu et al., 2013(a)] Itu, L. M., 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, Vol. 41, No. 4, pp. 669–681, 2013.

[Itu et al., 2013(b)] Itu, L. M., Sharma, P., Kamen, A., D., Suciu, C., Comaniciu, D. “A Novel Coupling Algorithm for Computing Blood Flow in Viscoelastic Arterial Models”, Proc. of the 35th Annual Inter. Conf. of the IEEE Engineering in Medicine & Biology Society - EMBC 2013, Osaka, Japan, Jul. 2013.

[Itu et al., 2013(c)] Itu, L. M., Sharma, P., Kamen, A., D., Suciu, C., Comaniciu, D. “GPU Accelerated One-Dimensional Blood Flow Computation in the Human Arterial Tree”, International Journal on Numerical Methods in Biomedical Engineering, published through early view, doi: 10.1002/cnm.2585.

[Kassab et al., 2006] Kassab, G. S. “Scaling Laws of Vascular Trees: of Form and Function”, American Journal of Physiology: Heart and Circulatory Physiology, Vol. 290, pp. 894–903, 2006.

[Kim et al., 2010] Kim, H. J., Vignon-Clementel, I. E., Coogan, J. S., Figueroa, C. A., Jansen, K. E., Taylor, C. A. “Patient-Specific Modeling of Blood Flow and Pressure in Human Coronary Arteries”, Annals of Biomedical Engineering, Vol. 38, pp. 3195–3209, 2010.

[Kumar et al., 2003] Kumar, R., Quateroni, A., Formaggia, L., Lamponi, D. “On Parallel Computation of Blood Flow in Human Arterial Network Based on 1-D Modelling”, Computing, Vol. 71, pp. 321–351, 2003.

[Mantero et al., 1992] Mantero, S., Pietrabissa, R., Fumero, R. “The Coronary Bed and its Role in the Cardiovascular System: a Review and an Introductory Single-Branch Model”, Journal of Biomedical Engineering, Vol. 14, pp. 109–116, 1992.

[Murray, 1926(a)] 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.

[Olufsen et al., 2000] Olufsen, M., Peskin, C., Kim, W. Y., Pedersen, E. “Numerical Simulation and Experimental Validation of Blood Flow in Arteries with Structured-Tree Outflow Conditions”, Annals of Biomedical Engineering, Vol. 28, pp. 1281–1299, 2000.

[Passerini, 2009(a)] Passerini, T. “Computational Hemodynamics of the Cerebral Circulation: Multiscale Modeling from the Circle of Willis to Cerebral Aneurysms”, Ph.D. Thesis, Politecnico di Milano, Italy, 2009.

[Pijls et al., 2000] Pijls, N. H., De Bruyne, B. “Coronary Pressure”, Kluwer Academic Publishing, London, UK, 2000.

[Reymond et al., 2011] Reymond, P., Bohraus, Y., Perren, F., Lazeyras, F., Stergiopulos, N. “Validation of a Patient-Specific One-Dimensional Model of the Systemic Arterial Tree”, American Journal of Physiology - Heart and Circulatory Physiology, Vol. 301, pp. 1173–1182, 2011.

[Stergiopulos et al., 1992] Stergiopulos, N., Young, D., Rogge, T. “Computer Simulation of Arterial Flow with Applications to Arterial and Aortic Stenosis”, Journal of Biomechanics, Vol. 25, pp. 1477–1488, 1992.

[Taylor et al., 2010] Taylor, C. A., Steinman D. A. “Image-Based Modeling of Blood Flow and Vessel Wall Dynamics: Applications, Methods and Future Directions”, Annals of Biomedical Engineering, Vol. 38, pp. 1188–1203, 2010.

[Vignon-Clementel et al., 2006] Vignon-Clementel, I., Figueroa, C. A., Jansen, K. E., Taylor, C. A. “Outflow Boundary Conditions for Three-Dimensional Finite Element Modeling of Blood Flow and Pressure in Arteries”, Comput. Methods Applied Mechanical Engineering, Vol. 195, pp. 3776–3796, 2006.

[Wilson et al., 1990] Wilson, R. F., Wyche, K., Christensen, B. V., Zimmer, S., Laxson, D. D. “Effects of Adenosine on Human Coronary Arterial Circulation”, Circulation, Vol. 82, pp. 1595–1606, 1990.


68

REZUMAT Simulările hemodinamice reprezintă o abordare modernă pentru analiza comportamentului sistemului cardiovascular. Acestea permit evaluarea non-invazivă a indicatorilor de diagnoză corespunzători unor diverse stări patologice. În lucrarea de faţă se prezintă dezvoltarea, implementarea, testarea şi validarea unui model multiscalar al circulaţiei coronariene pentru diagnosticarea non-invazivă a stenozelor coronariene. Având în vedere caracterul general al modelului multiscalar dezvoltat, obiectivul secundar al tezei a fost acela de a se dezvolta, implementa, testa şi valida un model multiscalar pentru diagnosticarea non-invazivă a coarctaţiilor aortice. Modelele propuse au fost dezvoltate cu scopul de fi utilizate în clinici de specialitate pentru diagnosticarea non-invazivă a afecţiunilor cardiovasculare ale pacienţilor. În acest sens este nevoie ca rezultatele să fie obţinute într-un timp cât mai scurt, pentru a putea diagnostica pacientul mai rapid, dar şi pentru a putea rula mai multe simulări într-un anumit interval de timp. Prin urmare, în vederea reducerii timpului de execuţie, s-a realizat implementarea pe un procesor paralel (procesor grafic – GPU) a soluţiei numerice a modelului cvasi unidimensional. În vederea realizării unor simulări hemodinamice personalizate, este necesar ca nu doar geometria arterială să fie achiziţionată de la pacient, ci şi parametri simulărilor să fie personalizaţi. Astfel, în cazul simulărilor specifice circulaţiei coronariene s-au estimat, pe de o parte, valorile parametrilor modelelor microvasculare, iar, pe de altă parte, valorile parametrilor modelului simplificat al inimii. În cazul evaluării non-invazive a coarctaţiilor aortice, pe lângă parametri condiţiilor de frontieră de ieşire, s-au personalizat şi proprietăţile pereţilor arteriali. Suplimentar, în vederea obţinerii în cadrul simulărilor hemodinamice a unor valori de presiune/debit egale cu cele achiziţionate pe cale non-invazivă de la pacient, s-au dezvoltat proceduri de optimizare iterative, complet automatizate. Atât modelul multiscalar al circulaţiei coronariene cât şi modelul coarctaţiei aortice au fost validate pe baza datelor achiziţionate de la pacienţi.

ABSTRACT Hemodynamic computations represent a modern approach for the analysis of the cardiovascular system. They enable the non-invasive assessment of diagnostic parameters corresponding to various pathologic conditions. The present work focuses on the development, implementation, testing and validation of a reduced-order multiscale model of the coronary circulation for the non-invasive diagnosis of coronary stenoses. Based on the properties of this multiscale model, a secondary goal of the thesis was to develop, implement, test and validate a multiscale model for the non-invasive diagnosis of aortic coarctations. The final goal of the developed models is for them to be applied in a clinical setting for the non-invasive, patient-specific assessment of cardiovascular pathologies. Thus, execution time becomes a crucial aspect, on one hand to diagnose a patient faster, and, on the other hand, to run more computations in a certain amount of time. Hence, to reduce the execution time, the numerical solution of the quasi one-dimensional model, which represents the main component of the multiscale models developed herein, was implemented on a parallel processor (graphics processing unit). To perform patient-specific hemodynamic computations, it is of outmost importance to not only acquire the patient-specific arterial geometry, but to also personalize the parameters of the computations. For the hemodynamic computations of the coronary circulation, both the values of the parameters of the microvascular models and the values of the parameters of the reduced-order heart model were estimated. For the non-invasive assessment of aortic coarctations, besides the parameters of the outlet boundary conditions, also the parameters of the arterial wall model were personalized. Additionally, to match non-invasively acquired patient-specific pressure and flow rate values, fully automated, iterative parameter estimation methods were developed. Both the multiscale model of the coronary circulation, and the multiscale model of aortic coarctation were validated using patient-specific data.


69

Curriculum vitae

Educaţie 2010-2013: studii de doctorat


Facultatea de Inginerie Electrică şi Ştiinţa Calculatoarelor

Conducător ştiinţific: prof. dr. ing. Florin MOLDOVEANU

2004-2009: studii de licenţă


Facultatea de Inginerie Electrică şi Ştiinţa Calculatoarelor

Specializarea: Automatică şi Informatică Aplicată

Experienţă profesională

Februarie ÷ aprilie 2011, august ÷ noiembrie 2011, iulie ÷ septembrie

2012

Stagiu de pregătire la centrul de cercetarea Siemens Corporate

Research, Princeton, New Jersey, USA, Departamentul de Imaging and

Computer Vision.

Ianuarie 2011 ÷ prezent

Bursă de cercetare din partea Siemens SRL, Braşov, România,

Departamentul Corporate Technology

Membru în societăţi profesionale

Societatea Română de Automatică şi Informatică Tehnică

Limbi străine Engleză – avansat

Germană – avansat

Domenii de interes Simulări hemodinamice

Modelare multiscalară de ordin redus

Modelarea circulaţiei coronariane

Evaluarea non-invazivă a patologiilor cardiovasculare

Procesare paralelă

Competenţe pe

calculator

Programare: C/C++, Java, Microsoft Visual Studi, NI

LabWindows/CVI, Eclipse IDE.

Simulare: MATLAB/SIMULINK, ANSYS Fluent.

Programare PLC: STEP 7, CoDeSys, AC1131, GX IEC Developer.


70

Curriculum vitae

Education 2010-2013: Ph. D.

Transilvania University of Braşov

Faculty of Electrical Engineering and Computer Science

Ph. D. supervisor: prof. dr. eng. Florin MOLDOVEANU

2004-2009: long term higher education

Transilvania University of Brasov

Faculty of Electrical Engineering and Computer Science

Specialization: Computer and System Engineering

Professional

experience

February ÷ April 2011, August ÷ November 2011, July ÷ September

2012

Internship at Siemens Corporate Research, Princeton, New Jersey, USA,

Departament of Imaging and Computer Vision.

January 2011 ÷ present

Research scholarship from Siemens SRL, Braşov, România,

Departament of Corporate Technology

Professional societies Romanian Automation and Technical Informatics Society

Foreign languages English – advanced

German – advanced

Research interests Hemodynamic computations

Reduced-order multiscale modeling

Modeling of the coronary circulation

Non-invasive assessment of cardiovascular diseases

Parallel Processing

Computer skills Programming: C/C++, Java, Microsoft Visual Studi, NI

LabWindows/CVI, Eclipse IDE.

Simulation: MATLAB/SIMULINK, ANSYS Fluent.

PLC Programming: STEP 7, CoDeSys, AC1131, GX IEC Developer.

Universitatea din Craiova - old.unitbv.roold.unitbv.ro/Portals/31/Sustineri de...

Documents

Transcript of Universitatea din Craiova - old.unitbv.roold.unitbv.ro/Portals/31/Sustineri de...