Model tlustého střeva
Pro vývoj modelů a související simulace byl použit software konečných prvků Abaqus (verze 6.13, SIMULIA, Providence, RI). Všechny simulace byly provedeny na pracovní stanici Dell T7810 se dvěma procesory Intel Xeon a 32 GB paměti. V naší předchozí studii chování prasečí tkáně jsme vytvořili anizotropní materiálový konstitutivní model pro sestupnou a spirální oblast prasečího tlustého střeva na základě současných testů nafukování a roztahování15. V několika předchozích studiích se ukázalo, že tento model adekvátně zachycuje anizotropní chování tlustého střeva19,34,35,36. Materiálové parametry stanovené v naší předchozí studii byly zachovány, protože se podle našich znalostí jedná o jedinou studii modelující tlusté střevo velkých zvířat při nafukování a roztahování, což jsou dva hlavní způsoby deformace tlustého střeva. Navíc bylo zjištěno, že prasečí tlusté střevo je vhodné pro studium divertikulózy33,37. Zvolený materiálový model a parametry jsou tedy pro tuto studii vhodnou volbou. Vzhledem k tomu, že anatomie spirálního tračníku není pro člověka relevantní, byly zde pro charakterizaci materiálu tkáně tlustého střeva použity výsledky ze vzorků sestupného tračníku. V naší předchozí práci jsme uvedli, že lokálně podél délky sestupného segmentu tlustého střeva prasete dochází k zanedbatelné variabilitě průměru, proto bylo normální tlusté střevo modelováno jako válcový segment (obr. 1a). Vnitřní a vnější poloměr byl stanoven na \({R}_{i}=11,5\) mm a \({R}_{o}=12,7\) mm na základě průměrných hodnot naměřených v naší předchozí práci. Axiální délka byla stanovena na \(L=10\) cm tak, aby byla několikanásobně větší než poloměr. Konkrétně se předpokládalo, že tkáň tlustého střeva je homogenní, anizotropní a nestlačitelný hyperelastický materiál, jehož mechanické chování popisuje následující anizotropní funkce deformační energie \(\bar{W}:\)15,35
První člen představuje neo-.Hookovu odezvu charakterizující chování nekolagenních složek tkáně tlustého střeva. Veličina \({I}_{1}={\lambda }_{z}^{2}+{\lambda }_{\theta }^{2}+\frac{1}{{\lambda }_{z}^{2}{\lambda }_{\theta }^{2}}}) představuje první. deformačního tenzoru, přičemž \({\lambda }_{z}\) a \({\lambda }_{\theta }\) se vztahují k osovému a obvodovému protažení, v tomto pořadí. Další členy charakterizují příspěvek kolagenových vláken, přičemž \({k}_{1}\) a \({k}_{2}\) jsou mírou tuhosti vláken. Konkrétně druhý člen (s horním indexem \(l\)) charakterizuje chování kolagenových vláken uspořádaných v podélném (axiálním) směru, zatímco třetí člen (s horním indexem \(s\)) charakterizuje odezvu kolagenových vláken rozptýlených ve dvou preferenčních symetrických směrech vzhledem k obvodovému směru. Veličina \({I}_{4}\) se označuje jako pseudoinvariantní a charakterizuje mechanickou odezvu vláken podél preferenčních směrů. Pro druhý člen platí, že \({I}_{4}^{l}={\lambda }_{z}^{2}\), protože zohledňuje vlákna podél podélného směru. Třetí člen \({I}_{4}^{s}\) je vyjádřen jako:
Tady, \(\pm {\gamma }^{s}\) označuje úhel dvou symetrických úhlů orientace vlákna vzhledem k obvodovému směru. V této práci bylo pro materiálové parametry \({C}_{10},\,{k}_{1}^{l},\,{k}_{2}^{l},\,{k}_{1}^{s},\,{k}_{2}^{s},\) a \({\gamma }^{s}\) použito pět sad materiálových parametrů, které byly dříve odvozeny pro pět vzorků sestupného tračníku prasat. Tato funkce deformační energie není v programu Abaqus přímo k dispozici. Proto byla poprvé integrována implementací uživatelského podprogramu v jazyce Fortran s názvem UANISOHYPER_INV, který umožňuje definovat anizotropní hyperelastické chování materiálu pomocí invariantní formulace. Tento podprogram vyžadoval definování první a druhé derivace deformačně-energetické funkce (rovnice 1) vzhledem ke skalárním invariantům \({I}_{1}\) a \({I}_{4}\). Abychom ověřili, že podprogram byl implementován správně, a ověřili platnost výpočetního modelu tkáně tlustého střeva, porovnali jsme modelově odhadnuté hodnoty izochorického obvodového napětí \({\bar{\sigma }}_{\theta }^{m}\) a izochorického axiálního napětí \({\bar{\sigma }}_{z}^{m}\) na vnějším povrchu tlustého střeva s hodnotami získanými analyticky na základě funkce deformační energie a vypočtenými takto:
Koeficient determinace \({R}^{2}\) byl vypočten následujícím způsobem pro kvantifikaci odchylky mezi výpočetními a analytickými hodnotami:
kde \(A={\bar{\sigma }}_{z},\,{\bar{\sigma }}_{\theta }\) a index „\(avg\)“ označuje průměr analytických hodnot ze všech datových bodů \(n\). Hodnota \({R}^{2}\) blízká 1 znamená, že je globálně dosaženo dobré korelace.
Model divertiklu
Pro modely nemocného tlustého střeva byla struktura podobná divertiklu modelována jako průsečík koule o průměru \(D\) s válcovým tvarem normálního tlustého střeva ve výšce \(H\) nad povrchem. V této studii byl uvažován pouze jeden vak, aby bylo možné se zaměřit na vliv jediného vaku bez dalších variabilit, které by byly zahrnuty z více vaků (počet vaků, pozice podél tlustého střeva, relativní poloha vůči sobě navzájem, variabilita v individuální velikosti atd.). Pokud je nám známo, neexistuje žádná studie, která by systematicky uváděla geometrii divertiklů. Běžně se uvádí, že jejich průměr se obvykle pohybuje mezi 5 až 10 mm, aniž by bylo řádně uvedeno, kde se průměr měří38,39. Aby byly pokryty různé velikosti váčků, byly při simulacích uvažovány tři hodnoty počáteční výšky váčku \(H\) a počátečního průměru váčku \(D\): \(H=2\) mm, 4 mm a 6 mm a \(D\) = 8 mm, 10 mm a 12 mm. Kombinace těchto hodnot D a H vedly k 9 modelům nemocných s různými velikostmi váčků. Na průsečíku mezi válcem představujícím normální tlusté střevo a koulí představující vak byl zahrnut kruhový filet o poloměru \({r}_{f}=2\) mm, aby se odstranily ostré hrany, které se ve skutečné tkáni nemusí vyskytovat a mohou měnit hodnoty napětí. Předpokládalo se, že tloušťka tkáně divertiklu je stejná jako tloušťka tlustého střeva (tj. \({R}_{o}-{R}_{i}\)). Byly učiněny pokusy charakterizovat mechaniku tlustého střeva při divertikulóze23,30,40,41. Nebyly však stanoveny žádné materiálové vlastnosti specificky pro tkáň divertiklu (tj. charakterizace izolovaného divertikulárního vaku). Pro část vaku byly proto uvažovány vlastnosti podobné normální tkáni.
Ohraničující podmínky
Ohraničující podmínky byly stanoveny na základě naší předchozí studie prasečího tlustého střeva15. Pozorovali jsme, že hodnoty luminálního tlaku (tlak na vnitřní stěnu tkáně při vnějším tlaku rovném nule) vyšší než 1,5 kPa vedly za pasivních podmínek k trvalému poškození tkáně. Ve výpočtových simulacích byla tato hodnota dále posunuta, aby byl znázorněn trend vývoje hodnot napětí, a na vnitřní povrch tlustého střeva včetně vaku byl vyvinut luminální tlak až 5 kPa. Tlak na vnější stěnu byl udržován na hodnotě 0. Vzhledem k tomu, že pracujeme s tenkostěnnou tkání (poměr vnitřního průměru a tloušťky se blíží 20), mohl být luminální tlak tedy asimilován jako delta mezi vnitřním tlakem (např. tlakem vyvolaným uvnitř průchodu tkáně fekáliemi nebo střevním plynem) a vnějším tlakem (např. břišním tlakem). Obvykle je k otevření tkáně tlustého střeva, která by jinak zkolabovala, zapotřebí malý tlak. Tento tlak (~0,2 kPa) je však ve srovnání se zde uvažovaným rozsahem velmi malý, takže je zanedbán. Tlak byl aplikován v kvazistatických přírůstcích. Během naší předchozí studie bylo v tlustém střevě in vivo obvykle pozorováno axiální roztažení o velikosti přibližně 10 %. Proto bylo před zvýšením tlaku zavedeno 10% axiální protažení ve směru \(Z\). Ve tkáni prasečího tlustého střeva jsme zjistili velmi malé hodnoty úhlu otevření, a proto byla v modelu zanedbána zbytková obvodová napětí. U normálních i nemocných modelů jsme použili symetrii vzhledem k rovině \(YZ\) a rovině \(XY\) a řešili jsme pouze čtvrtinu modelu, abychom dosáhli vyšší efektivity, a poté jsme celý model rekonstruovali pro účely vizualizace.
Mezování
V každém případě bylo provedeno vhodné rozdělení geometrie, aby bylo zajištěno, že je možné použít síťování pomocí lineárních šestistěnných prvků. Pro vynucení nestlačitelnosti materiálu byla zvolena formulace hybridních prvků (tj. prvek C3D8H z knihovny Abaqus). Pro normálový model byla provedena studie konvergence sítě zmenšováním velikosti referenční sítě počínaje výchozí velikostí automaticky navrženou programem Abaqus (1,3 mm), dokud odchylka maximální hodnoty maximálních hlavních napětí mezi dvěma po sobě jdoucími sítěmi nebyla menší než 10 %. Konvergující velikost sítě byla poté použita i v modelech nemocných s některými dalšími upřesněními v okolí vaku. Pro normální model bylo použito celkem 69608 prvků a 88350 uzlů při nejjemnější velikosti sítě zachované pro analýzu (čtvrtina geometrie). Pro devět modelů s nemocí bylo použito 69764 prvků a 80656 prvků (88605, resp. 102440 uzlů) (čtvrtina geometrie). Lokální orientace materiálu byla zadána pomocí možnosti diskrétní formulace dostupné v programu Abaqus, aby byla správně definována orientace vláken v každém prvku podle požadavků funkce deformační energie.
Výsledky analýzy
Bylo uvažováno celkem 10 modelů (jeden normální a devět nemocných s různou velikostí vaku). Každý model byl simulován s 5 různými sadami materiálových parametrů, což dalo celkem 50 simulovaných případů. Pro porovnání jednotlivých modelů je nutné stanovit míru napětí. Vzhledem k tomu, že nestanovujeme kritérium selhání, ale pouze sledujeme odchylky za různých podmínek, byla by vhodná jakákoli míra, například maximální hlavní napětí nebo von Misesovo napětí. Výsledky zde uvádíme pouze z hlediska maximálního hlavního napětí, abychom se vyhnuli nadbytečným informacím, protože celková pozorování s von Misesovým napětím byla shledána podobnými. Analýza a výsledky s von Misesovým napětím jsou zahrnuty v údajích spojených s touto studií, pokud si je čtenář přeje prohlédnout. Pro každý simulovaný případ bylo maximální hlavní napětí (v centroidu prvku sítě), označované jako \({\sigma }_{MPS}\), sledováno v rozsahu tlaků 0 až 5 kPa s přírůstkem 0,25 kPa. Byl vytvořen skript v jazyce Python, který automaticky vyhledává maximální hodnotu \({\sigma }_{MPS}\), označovanou jako \({\sigma }_{MPS,max}\) při různých tlacích pro každý simulovaný případ. Pro každý model byla vypočtena průměrná hodnota získaná pro každý simulovaný případ (tj. pro všech pět sad materiálových parametrů) a ponechána pro následnou analýzu. Tyto hodnoty se označují jako \({\sigma }_{MPS,\,max}^{avg}\) a \({\sigma }_{MPS,\,max}^{avg,\,n}\) pro nemocné a normální modely.
Hodnoty \({\sigma }_{MPS,\,max}^{avg}\) byly pro tlak v rozmezí 1 až 5 kPa (s přírůstkem 0,25 kPa) korelovány s různými parametry geometrie vaku (obr. 1c), aby se zjistil vztah mezi geometrií vaku a vysokým napětím: šířka hrdla vaku podél osy \(x\) \({D}_{x}\), šířka hrdla vaku podél osy \(z\) \({D}_{z}\), plocha na bázi hrdla vaku \({A}_{b}\) a povrch vaku na straně lumen \({S}_{p}\). Pro každý parametr byl vypočten Pearsonův korelační koeficient \({r}^{2}\). K vyhodnocení významnosti korelace byla rovněž vypočtena dvouvýběrová p-hodnota.
Dále jsme vypočítali tlak \({P}_{d}\) v modelech nemocných tak, že pro daný tlak \({P}_{n}\) v normálním tlustém střevě máme:
Následující tlakový faktor byl následně vypočten jako:
Tento faktor udává %, o které se musí snížit tlak \({P}_{n}\)v nemocném tlustém střevě, aby \({\sigma }_{MPS,byl podobný \({\sigma }_{MPS,\,max}^{n}\), tj.maximální hodnota maximálního hlavního napětí v nemocném modelu nepřesahuje odpovídající hodnotu v normálním modelu. Redukovaný tlak se pak označuje jako \({P}_{d}\). Konkrétně pro každý simulovaný případ bylo \({P}_{n}\) sníženo o 0,25 kPa, dokud nebyl splněn vztah z rovnice 6, čímž byla získána odpovídající hodnota \({P}_{d}\) a vypočten tlakový faktor \(f\). To bylo provedeno pro \({P}_{n}\) v rozsahu 2 až 5 kPa (po 1 kPa). Byla vypočtena průměrná hodnota získaná pro modely ve všech simulovaných případech a ponechána pro další analýzu.
Nakonec jsme stanovili zónu vlivu vaku, tj. vzdálenost kolem vaku, kde je zvýšené napětí. Abychom toho dosáhli, stanovili jsme pro každý nemocný případ mapu maximálního hlavního napětí na bočním povrchu lumen tlustého střeva (protože na tomto povrchu je napětí nejvyšší) podél uzlů sítě umístěných na podélných a obvodových drahách procházejících středem vaku. Vývoj maximálního hlavního napětí podél těchto cest byl sledován jako funkce vzdáleností \({d}_{l}\) a \({d}_{c}\) měřených vzhledem ke středu vaku. Poté jsme vypočítali nejvzdálenější vzdálenosti \({d}_{l}^{inf}\) podél podélné dráhy a \({d}_{c}^{inf}\) podél obvodové dráhy, po kterých maximální hlavní napětí klesne na 10 % odpovídající hodnoty v normálovém modelu podél stejných drah. Tyto vzdálenosti představují naši míru zóny vlivu. Vzhledem k tomu, že po dosažení určitého tlaku dosáhnou hodnoty plató, byly jejich průměrné hodnoty ve všech případech pro každý model při tlaku 5 kPa korelovány s různými parametry geometrie vaku, podobně jako \({\sigma }_{MPS,\,max}^{avg}\).
Všechny výpočty byly provedeny v programovacím jazyce Python pomocí notebooků Jupyter42. Data byla ukládána a manipulována pomocí knihovny Pandas jazyka Python43. Základní matematické operace byly provedeny pomocí knihovny NumPy44. Korelační koeficient \({r}^{2}\) a dvouvýběrová p-hodnota byly vypočteny pomocí funkce „stats“ knihovny Scipy45. Všechny grafy byly vytvořeny pomocí knihovny Matplotlib46.
.