Model jelita grubego
Do opracowania modeli i przeprowadzenia związanych z nimi symulacji wykorzystano oprogramowanie Abaqus (wersja 6.13, SIMULIA, Providence, RI). Wszystkie symulacje przeprowadzono na stacji roboczej Dell Workstation T7810 z dwoma procesorami Intel Xeon i 32 GB pamięci. W naszych poprzednich badaniach nad zachowaniem tkanki świńskiej, stworzyliśmy anizotropowy model konstytutywny materiału dla zstępujących i spiralnych regionów okrężnicy świńskiej w oparciu o jednoczesne testy inflacji i rozciągania15. W kilku poprzednich badaniach wykazano, że model ten odpowiednio oddaje anizotropowe zachowanie okrężnicy19,34,35,36. Parametry materiałowe ustalone w naszym poprzednim badaniu zostały zachowane, ponieważ, według naszej wiedzy, jest to jedyne badanie modelujące okrężnicę dużych zwierząt w warunkach inflacji i rozszerzania, dwóch głównych trybów deformacji okrężnicy. Ponadto ustalono, że okrężnica świń jest odpowiednia do badania uchyłkowatości33,37. Wybrany model materiału i parametry są zatem odpowiednie dla tego badania. Ponieważ anatomia okrężnicy spiralnej nie jest odpowiednia dla ludzi, wyniki z próbek okrężnicy zstępującej zostały użyte do scharakteryzowania materiału tkankowego okrężnicy. W naszej poprzedniej pracy wykazaliśmy, że lokalna zmiana średnicy na długości segmentu okrężnicy zstępującej u świń jest nieistotna, dlatego okrężnica normalna została zamodelowana jako segment cylindryczny (Rys. 1a). Promień wewnętrzny i zewnętrzny zostały ustalone odpowiednio na 11,5 mm i 12,7 mm, na podstawie średnich wartości zmierzonych w naszej poprzedniej pracy. Długość osiowa została ustalona na L=10 cm, tak aby była kilka razy większa od promienia. Założono, że tkanka jelita grubego jest jednorodnym, anizotropowym i nieściśliwym materiałem hiperelastycznym, którego zachowanie mechaniczne jest opisane następującą anizotropową funkcją energii odkształcenia \(\bar{W}:\)15,35
Pierwszy człon reprezentuje odpowiedź Neo-Hookean’a charakteryzujący zachowanie niekolagenowych składników tkanki okrężnicy. Wielkość \({I}_{1}={lambda }_{z}^{2}+{lambda }_{z}^{2}}{{lambda }_{z}^{2}}}}{{lambda }_{z}^{2}}}}}) reprezentuje pierwszy pierwszy niezmiennik tensora deformacji, przy czym \({{lambda }_{z}}) i \({{lambda }_{theta }}) odnoszą się odpowiednio do rozciągnięcia osiowego i obwodowego, odpowiednio. Kolejne terminy charakteryzują udział włókien kolagenowych, przy czym \({k}_{1}} i \({k}_{2}}) są miarą sztywności włókien. W szczególności, drugi człon (z indeksem górnym) charakteryzuje zachowanie włókien kolagenowych ułożonych wzdłuż kierunku podłużnego (osiowego), podczas gdy trzeci człon (z indeksem górnym) charakteryzuje odpowiedź włókien kolagenowych rozproszonych w dwóch preferencyjnych, symetrycznych kierunkach w odniesieniu do kierunku obwodowego. Wielkość \({I}_{4}}}jest określana jako pseudo niezmienna i charakteryzuje odpowiedź mechaniczną włókien wzdłuż kierunków preferencyjnych. W przypadku drugiego członu, \({I}_{4}^{l}={lambda }_{z}^{2}}, ponieważ uwzględnia on włókna wzdłuż kierunku podłużnego. Trzeci człon, \({I}_{4}^{s}}, wyraża się jako:
Tutaj, \^{2}{gamma }^{s}$) oznacza kąt dwóch symetrycznych orientacji włókien w odniesieniu do kierunku obwodowego. W niniejszej pracy wykorzystano pięć zestawów parametrów materiałowych, uzyskanych wcześniej dla pięciu próbek okrężnicy zstępującej świni, dla parametrów materiałowych \({C}_{10},\, {k}_{1}^{l},\, {k}_{2}^{l},\, {k}_{1}^{s},\, {k}_{2}^{s},\) i \({gamma }^{s}}). Ta funkcja energii odkształcenia nie jest bezpośrednio dostępna w programie Abaqus. W związku z tym, po raz pierwszy zintegrowano ją z programem poprzez zaimplementowanie podprogramu UANISOHYPER_INV, który pozwala na zdefiniowanie anizotropowego, hiperelastycznego zachowania materiału z wykorzystaniem sformułowania niezmienniczego. Podprogram ten wymagał zdefiniowania pierwszej i drugiej pochodnej funkcji odkształcenie-energia (Eq. 1) względem skalarnych niezmienników \({I}_{1}\) oraz \({I}_{4}\). Aby sprawdzić, czy podprogram został poprawnie zaimplementowany i zwalidować obliczeniowy model tkanki okrężnicy, porównaliśmy oszacowane przez model wartości naprężeń izochorycznego naprężenia obwodowego \({bar{sigma }}_{theta }^{m}}) i izochorycznego naprężenia osiowego \({bar{sigma }}_{z}^{m}}) na zewnętrznej powierzchni okrężnicy z wartościami uzyskanymi analitycznie na podstawie funkcji energii odkształcenia i obliczonymi w następujący sposób:
Współczynnik determinacji {R}^{2} obliczono w następujący sposób, aby określić odchylenie między wartościami obliczeniowymi i analitycznymi:
gdzie {A={bar{sigma }}_{z},\), a indeks dolny „avg” oznacza średnią wartości analitycznych dla wszystkich punktów danych. Wartość \({R}^{2}} bliska 1 wskazuje, że dobra korelacja jest uzyskiwana globalnie.
Model uchyłka
Dla modeli chorej okrężnicy, struktura przypominająca uchyłek była modelowana jako przecięcie kuli o średnicy \(D) z cylindrycznym kształtem normalnej okrężnicy na wysokości \(H) nad powierzchnią. W tym badaniu uwzględniono tylko jeden worek, aby skupić się na efekcie pojedynczego worka bez dodatkowych zmiennych, które mogłyby być włączone z wielu worków (liczba worków, położenie wzdłuż okrężnicy, względne położenie względem siebie, zmienność w wielkości osobników, itp.) Według naszej wiedzy, nie ma badań systematycznie opisujących geometrię uchyłków. Powszechnie mówi się, że mają one zwykle od 5 do 10 mm średnicy, bez właściwego wskazania miejsca pomiaru średnicy38,39. Aby uwzględnić różne rozmiary worków, w symulacjach uwzględniono trzy wartości początkowej wysokości worka (H) i początkowej średnicy worka (D): \(H=2\) mm, 4 mm i 6 mm, oraz \(D\) = 8 mm, 10 mm i 12 mm. Kombinacje tych wartości D i H pozwoliły uzyskać 9 chorych modeli o różnych rozmiarach torebki. Na przecięciu cylindra reprezentującego normalną okrężnicę i kuli reprezentującej worek umieszczono okrągłe wypełnienie o promieniu 2 mm w celu usunięcia ostrych krawędzi, które mogą nie występować w rzeczywistej tkance i mogą zmieniać wartości naprężeń. Przyjęto, że grubość tkanki uchyłka jest taka sama jak grubość okrężnicy (tj. \({R}_{o}-{R}_{i}\)). Podejmowano próby scharakteryzowania mechaniki okrężnicy w uchyłkowatości23,30,40,41. Nie ustalono jednak właściwości materiału specyficznych dla tkanki uchyłka (tj. charakterystyki izolowanego worka uchyłkowego). Właściwości podobne do normalnej tkanki były więc brane pod uwagę dla przekroju worka.
Warunki brzegowe
Warunki brzegowe zostały nałożone w oparciu o nasze wcześniejsze badania okrężnicy świńskiej15. Zaobserwowaliśmy, że ciśnienie luminalne (ciśnienie na wewnętrznej ścianie tkanki przy ciśnieniu zewnętrznym równym zeru) o wartościach powyżej 1,5 kPa prowadziło do trwałego uszkodzenia tkanki w warunkach pasywnych. W symulacjach obliczeniowych, w celu zobrazowania trendu ewolucji wartości naprężeń, wartość tę przesunięto jeszcze dalej i na wewnętrzną powierzchnię okrężnicy, łącznie z torebką, wywierano ciśnienie luminalne do 5 kPa. Ciśnienie w ścianie zewnętrznej utrzymywano na poziomie 0. Ponieważ pracujemy z tkanką cienkościenną (stosunek średnicy wewnętrznej do grubości bliski 20), ciśnienie luminalne może być asymilowane jako delta pomiędzy ciśnieniem wewnętrznym (np. ciśnienie wywierane wewnątrz kanału tkanki przez masy kałowe lub gaz jelitowy) a ciśnieniem zewnętrznym (np. ciśnienie brzuszne). Zazwyczaj do otwarcia tkanki okrężnicy, która w przeciwnym razie zapadłaby się, potrzebne jest niewielkie ciśnienie. Ciśnienie to (~0,2 kPa) jest jednak bardzo małe w porównaniu z zakresem rozpatrywanym tutaj, więc zostało zignorowane. Ciśnienie było stosowane w quasi-statycznych przyrostach. Rozciąganie osiowe o około 10% było typowo obserwowane in vivo w okrężnicy podczas naszych poprzednich badań. Tak więc, 10% osiowe rozciągnięcie zostało narzucone wzdłuż kierunku \(Z) przed ciśnieniem. Stwierdziliśmy bardzo małe wartości kąta rozwarcia w tkance okrężnicy świńskiej, a szczątkowe naprężenia obwodowe zostały w związku z tym pominięte w modelu. Dla obu modeli, normalnego i chorego, zastosowaliśmy symetrię względem płaszczyzny Y i X i rozwiązaliśmy tylko jedną czwartą modelu, aby osiągnąć wyższą wydajność, przed rekonstrukcją całego modelu w celu wizualizacji.
Meshing
W każdym przypadku wykonano odpowiedni przekrój geometrii, aby zapewnić, że meshing jest możliwy przy użyciu liniowych elementów heksaedrycznych. Wybrano hybrydową formułę elementu w celu wymuszenia nieściśliwości materiału (tj. element C3D8H z biblioteki Abaqus). Badanie zbieżności siatki przeprowadzono dla modelu normalnego, zmniejszając rozmiar siatki referencyjnej począwszy od domyślnego rozmiaru sugerowanego automatycznie przez program Abaqus (1.3 mm) do momentu, gdy zmiana maksymalnej wartości maksymalnych naprężeń głównych pomiędzy dwoma kolejnymi siatkami była mniejsza niż 10%. Zbieżny rozmiar siatki zastosowano następnie również w modelach chorych, z kilkoma dodatkowymi udoskonaleniami wokół woreczka. W modelu normalnym wykorzystano łącznie 69608 elementów i 88350 węzłów przy najdrobniejszym rozmiarze oczka sieci zachowanym do analizy (jedna czwarta geometrii). W przypadku dziewięciu modeli chorych użyto od 69764 elementów do 80656 elementów (odpowiednio 88605 i 102440 węzłów) (jedna czwarta geometrii). Lokalna orientacja materiału została określona przy użyciu opcji formowania dyskretnego dostępnej w programie Abaqus w celu prawidłowego zdefiniowania orientacji włókien w każdym elemencie, zgodnie z wymaganiami funkcji energii odkształcenia.
Analiza wyników
Rozważono w sumie 10 modeli (jeden normalny i dziewięć chorych z różnymi rozmiarami woreczków). Każdy model był symulowany z 5 różnymi zestawami parametrów materiałowych, co dało w sumie 50 symulowanych przypadków. Miara poziomu naprężenia jest niezbędna do porównania modeli. Ponieważ nie ustalamy kryterium zniszczenia, a jedynie analizujemy zmiany w różnych warunkach, każda miara, taka jak maksymalne naprężenie główne lub naprężenie von Misesa, będzie odpowiednia. Wyniki przedstawiono tutaj tylko w odniesieniu do maksymalnego naprężenia głównego, aby uniknąć zbędnych informacji, ponieważ ogólne obserwacje dotyczące naprężenia von Misesa okazały się podobne. Analiza i wyniki dotyczące naprężeń von Misesa są zawarte w danych związanych z tym badaniem, jeśli czytelnik chce się z nimi zapoznać. Dla każdego symulowanego przypadku zaobserwowano maksymalne naprężenie główne (w centroidzie elementu siatki), określane jako \(\sigma }_{MPS}}), w zakresie ciśnienia od 0 do 5 kPa z przyrostem 0,25 kPa. Opracowano skrypt Pythona do automatycznego wyszukiwania maksymalnej wartości \(Sygma }_{MPS}}, oznaczonej jako \(Sygma }_{MPS,max}}) przy różnych ciśnieniach dla każdego symulowanego przypadku. Dla każdego modelu, średnia wartość uzyskana dla każdego symulowanego przypadku (tj. dla wszystkich pięciu zestawów parametrów materiałowych), została obliczona i zachowana do dalszej analizy. Wartości te są określane jako \(\MPS,\,max}^{avg}} i \(\MPS,\,max}^{avg,\,n}} odpowiednio dla modelu chorego i normalnego.
Wartości \({sigma }_{MPS,\,max}^{avg}}, zostały skorelowane, dla ciśnienia pomiędzy 1 i 5 kPa (w przyroście 0,25 kPa), z różnymi parametrami geometrii woreczka (Rys. 1c) w celu określenia zależności pomiędzy geometrią woreczka a wysokim naprężeniem: szerokość szyjki woreczka wzdłuż osi \(x) \({D}_{x}}, szerokość szyjki woreczka wzdłuż osi \(z} \({D}_{z}}), powierzchnia u podstawy szyjki woreczka \({A}_{b}}) oraz powierzchnia woreczka od strony światła \({S}_{p}}). Dla każdego parametru obliczono współczynnik korelacji Pearsona \({r}^{2}\). W celu oceny istotności korelacji obliczono również dwuogonową wartość p-value. Następujący współczynnik ciśnienia został następnie obliczony jako: Współczynnik ten wskazuje %, o jaki ciśnienie \({P}_{n}} musi być zmniejszone w chorym jelicie grubym, aby \({sigma }_{MPS,max}} jest podobne do \(\MPS,\,max}^{n}}, tzn.maksymalna wartość maksymalnych naprężeń głównych w modelu chorym nie przekracza odpowiadającej jej wartości w modelu normalnym. Zmniejszone ciśnienie jest wtedy oznaczane jako \({P}_{d}}). Dla każdego symulowanego przypadku obniżano ciśnienie o 0,25 kPa, aż do spełnienia zależności z równania 6, uzyskując odpowiednią wartość ciśnienia i obliczając współczynnik ciśnienia. Czynność tę wykonano dla \({P}_{n}} w zakresie od 2 do 5 kPa (z przyrostem co 1 kPa). Obliczano średnią wartość uzyskaną dla modeli we wszystkich symulowanych przypadkach i zachowywano ją do dalszej analizy. Na koniec określono strefę wpływu woreczka, tj. odległość wokół woreczka, w której naprężenia są podwyższone. Aby to osiągnąć, ustaliliśmy dla każdego chorego przypadku mapę maksymalnego naprężenia głównego na bocznej powierzchni światła jelita grubego (ponieważ naprężenie jest największe na tej powierzchni) wzdłuż węzłów siatkowania zlokalizowanych na ścieżkach wzdłużnych i obwodowych przechodzących przez środek woreczka. Ewolucja maksymalnego naprężenia głównego wzdłuż tych ścieżek była obserwowana jako funkcja odległości \({d}_{l}} i \({d}_{c}}} mierzonych w odniesieniu do środka woreczka. Następnie obliczyliśmy najdalsze odległości \({d}_{l}^{inf}} wzdłuż ścieżki podłużnej i \({d}_{c}^{inf}} wzdłuż ścieżki obwodowej, po których maksymalne naprężenie główne spada w granicach 10% odpowiadającej mu wartości w modelu normalnym wzdłuż tych samych ścieżek. Odległości te stanowiły naszą miarę strefy wpływu. Obserwując, że osiągają one wartość plateau po pewnym ciśnieniu, ich średnie wartości dla wszystkich przypadków dla każdego modelu przy ciśnieniu 5 kPa skorelowano z różnymi parametrami geometrii woreczka, podobnie jak w przypadku ∗(sigma }_{MPS,∗,max}^{avg}}). Wszystkie obliczenia wykonano w języku programowania Python, tworząc notatniki Jupyter42. Dane były przechowywane i manipulowane za pomocą biblioteki Pandas języka Python43. Podstawowe operacje matematyczne wykonywano za pomocą biblioteki NumPy44. Współczynnik korelacji \({r}^{2}\) oraz dwuogonową wartość p-value obliczono za pomocą funkcji „stats” biblioteki Scipy45. Wszystkie wykresy zostały zrealizowane za pomocą biblioteki Matplotlib46. .