Kolonmodell
Den finita elementprogramvaran Abaqus (version 6.13, SIMULIA, Providence, RI) användes för att utveckla modellerna och genomföra tillhörande simuleringar. Alla simuleringar utfördes på en Dell Workstation T7810 med två Intel Xeon-processorer och 32 GByte minne. I vår tidigare studie av vävnadsbeteende hos svin upprättade vi en konstitutiv modell av anisotropt material för de nedåtgående och spiralformade regionerna i tjocktarmen hos svin, baserad på samtidiga inflations- och utvidgningstester15. Denna modell har i flera tidigare studier visat sig på ett adekvat sätt fånga kolonets anisotropa beteende19,34,35,36. De materialparametrar som fastställdes i vår tidigare studie bibehölls eftersom det, såvitt vi vet, är den enda studien som modellerar ett stort djurs tjocktarm under uppblåsning och utvidgning, de två viktigaste deformationsformerna i tjocktarmen. Dessutom har det konstaterats att tjocktarmen hos svin är relevant för studier av divertikulosor33,37. Den valda materialmodellen och de valda parametrarna är därför lämpliga val för denna studie. Eftersom spiralkolonets anatomi inte är relevant för människor användes här resultat från prover från det nedåtgående kolonet för att karaktärisera kolonvävnadsmaterialet. Vi har i vårt tidigare arbete rapporterat att det finns en försumbar variation i diametern lokalt längs längden på ett nedåtgående kolon-segment från svin och därför modellerades det normala kolonet som ett cylindriskt segment (fig. 1a). Den inre och yttre radien sattes till \({R}_{i}=11,5\) mm respektive \({R}_{o}=12,7\) mm, baserat på genomsnittliga värden som uppmätts i vårt tidigare arbete. Den axiella längden sattes till \(L=10\) cm så att den är flera gånger större än radien. Kolonvävnaden antogs vara ett homogent, anisotropt och inkompressibelt hyperelastiskt material vars mekaniska beteende beskrivs av följande anisotropa töjningsenergifunktion \(\(\bar{W}):\)15,35
Den första termen representerar en Neo-Hookean respons som karakteriserar beteendet hos de icke-kollagena beståndsdelarna i tjocktarmsvävnaden. Mängden \({I}_{1}={\lambda }_{z}^{2}+{\lambda }_{\theta }^{2}+\frac{1}{{{\lambda }_{z}^{2}{\lambda }_{\theta }^{2}}}\) representerar det första invariant av deformationstensorn med \({\lambda }_{z}\) och \({\lambda }_{\theta }\) som hänvisar till den axiella och cirkumferentiella sträckningen, respektive. De efterföljande termerna karakteriserar bidraget från kollagenfibrerna, med \({k}_{1}\) och \({k}_{2}\) som ett mått på fiberstyvhet. Den andra termen (med översatt \(l\\)) karakteriserar beteendet hos kollagenfibrerna som är inriktade längs den longitudinella (axiella) riktningen, medan den tredje termen (med översatt \(s\)) karakteriserar responsen hos kollagenfibrerna som är utspridda i två preferentiella symmetriska riktningar med avseende på den cirkumferentiella riktningen. Mängden \({I}_{4}\) kallas pseudoinvariant och karakteriserar den mekaniska responsen hos fibrerna längs de preferentiella riktningarna. För den andra termen gäller \({I}_{4}^{l}={\lambda }_{z}^{2}\) eftersom den tar hänsyn till fibrerna längs den longitudinella riktningen. Den tredje termen \({I}_{4}^{s}\) uttrycks som:
Här, \(\pm {\gamma }^{s}\) anger vinkeln mellan de två symmetriska fiberorienteringarna i förhållande till omkretsriktningen. Fem uppsättningar materialparametrar, som tidigare tagits fram för fem prover av nedåtgående kolon hos svin, användes i detta arbete för materialparametrarna \({C}_{10},\,{k}_{1}^{l},\,{k}_{2}^{l},\,{k}_{1}^{s},\,{k}_{2}^{s},\) och \({\gamma }^{s}\). Denna funktion för töjningsenergi är inte direkt tillgänglig i Abaqus. Den integrerades därför för första gången genom att implementera en Fortranbaserad användarunderrutin kallad UANISOHYPER_INV, som gör det möjligt att definiera ett anisotropt hyperelastiskt materialbeteende med hjälp av den invarianta formuleringen. Denna underrutin krävde att man definierade den första och andra derivatan av töjningsenergifunktionen (ekv. 1) med avseende på de skalära invarianterna \({I}_{1}\) och \({I}_{4}\). För att verifiera att underrutinen implementerades korrekt och validera beräkningsmodellen för tjocktarmsvävnaden jämförde vi modellskattade spänningsvärden för isokorisk cirkulationsspänning \({\bar{\sigma }}_{\theta }^{m}\) och isokorisk axialspänning \({\bar{\sigma }}_{z}^{m}}\) på tjocktarmens yttre yta med värden som erhållits analytiskt baserat på töjningsenergifunktionen och som beräknats på följande sätt:
Determinationskoefficienten \({R}^{2}\) beräknades på följande sätt för att kvantifiera avvikelsen mellan beräknings- och analysvärden:
där \(A={{\bar{\sigma }}_{z},\,{\bar{\sigma }}_{\theta }\), och subscript ”\(avg\)” anger medelvärdet av de analytiska värdena över alla \(n\) datapunkter. Ett värde på \({R}^{2}\) nära 1 indikerar att en god korrelation erhålls globalt.
Divertikelmodell
För modellerna av den sjuka tjocktarmen modellerades en divertikelliknande struktur som skärningspunkten mellan en sfär med diametern \(D\) och den cylindriska formen av den normala tjocktarmen på en höjd \(H\) över ytan. Endast en påse beaktades i denna studie för att fokusera på effekten av en enda påse utan ytterligare variationer som skulle kunna införlivas av flera påsar (antal påsar, positioner längs tjocktarmen, relativa positioner i förhållande till varandra, variationer i individernas storlek etc.). Såvitt vi vet finns det ingen studie som systematiskt rapporterar divertiklarnas geometri. Det nämns ofta att de vanligtvis är mellan 5 och 10 mm i diameter utan att det anges ordentligt var diametern mäts38,39. För att täcka in en mängd olika storlekar på poucherna beaktades tre värden i simuleringarna för den initiala pouchhöjden \(H\) och den initiala pouchdiametern \(D\): \(H=2\) mm, 4 mm och 6 mm och \(D\) = 8 mm, 10 mm och 12 mm. Kombinationer av dessa värden för D och H ledde till nio sjuka modeller med olika storlekar på pungen. En rund filé med radie \({r}_{f}=2\) mm inkluderades i skärningspunkten mellan cylindern som representerar den normala tjocktarmen och klotet som representerar påsen för att avlägsna skarpa kanter som kanske inte förekommer i den verkliga vävnaden och som kan förändra spänningsvärdena. Divertikelvävnadens tjocklek antogs vara densamma som tjocktarmen (dvs. \({R}_{o}-{R}_{i}\)). Det har gjorts försök att karakterisera kolonets mekanik vid divertikulos23,30,40,41. Inga materialegenskaper har dock fastställts specifikt för divertikelvävnaden (dvs. karakterisering av den isolerade divertikelpåsen). Egenskaper som liknar den normala vävnaden ansågs därför gälla för pouchsektionen.
Gränsvillkor
Gränsvillkoren infördes baserat på vår tidigare studie av tjocktarmen hos svin15. Vi observerade att luminaltryck (tryck på vävnadens innervägg med yttre tryck lika med noll) värden över 1,5 kPa ledde till permanenta skador på vävnaden under passiva förhållanden. I beräkningssimuleringarna ökades detta värde ytterligare för att visualisera trenden för utvecklingen av stressvärdena, och ett luminalt tryck på upp till 5 kPa påtvingades den inre ytan av tjocktarmen, inklusive pouch. Det yttre väggtrycket hölls lika med 0. Eftersom vi arbetar med en tunnväggig vävnad (förhållandet mellan innerdiameter och tjocklek ligger nära 20) kan det luminala trycket således likställas med deltatrycket mellan det inre trycket (t.ex. det tryck som påförs inuti vävnadens kanal av fekalier eller tarmgaser) och det yttre trycket (t.ex. trycket i buken). Typiskt sett krävs ett litet tryck för att öppna upp tjocktarmsvävnaden som annars skulle kollapsa. Detta tryck (~0,2 kPa) är dock mycket litet jämfört med det område som beaktas här, så det ignoreras. Trycket applicerades i kvasistatiska steg. En axial sträckning på cirka 10 % observerades vanligtvis in vivo i tjocktarmen under vår tidigare studie. En 10-procentig axiell sträckning infördes därför i \(Z\)-riktningen innan trycket sattes på. Vi har funnit mycket små öppningsvinkelvärden i tjocktarmsvävnaden hos svin, och resterande cirkumferentiella spänningar försummades därför i modellen. För både normala och sjuka modeller använde vi symmetri med avseende på \(YZ\)-planet och \(XY\)-planet och löste endast en fjärdedel av modellen för att uppnå högre effektivitet, innan vi rekonstruerade hela modellen för visualiseringsändamål.
Mätning
Att lämplig geometri sektionering gjordes i varje enskilt fall för att säkerställa att nätning är möjlig med hjälp av linjära hexaederelement. Hybridelementformulering valdes för att framtvinga materialets inkompressibilitet (dvs. element C3D8H från Abaqus-biblioteket). En undersökning av nätkonvergens genomfördes för den normala modellen genom att minska referensmaskstorleken från den standardstorlek som föreslås automatiskt av Abaqus (1,3 mm) tills variationen av det maximala värdet av de maximala huvudspänningarna mellan två på varandra följande maskor var mindre än 10 %. Den konvergerande maskstorleken användes sedan även i de sjuka modellerna, med några ytterligare förfiningar runt påsarna. Totalt användes 69608 element och 88350 noder för den normala modellen med den finaste maskstorlek som behölls för analys (en fjärdedel av geometrin). Mellan 69764 element och 80656 element (88605 respektive 102440 noder) användes för de nio sjukdomsmodellerna (en fjärdedel av geometrin). En lokal materialorientering specificerades med hjälp av det diskreta formuleringsalternativet som finns tillgängligt i Abaqus för att korrekt definiera fiberorienteringen i varje element så som krävs i töjningsenergifunktionen.
Resultatanalys
Totalt 10 modeller (en normal och nio sjuka med olika storlekar på på påsarna) beaktades. Varje modell simulerades med 5 olika uppsättningar materialparametrar vilket resulterade i totalt 50 simulerade fall. Ett mått på stressnivå är nödvändigt för att kunna jämföra mellan modellerna. Eftersom vi inte fastställer ett kriterium för brott utan bara tittar på variationer under olika förhållanden, skulle ett mått som maximal huvudspänning eller von Mises-spänning vara lämpligt. Resultaten visas här endast med avseende på maximal huvudspänning för att undvika överflödig information, eftersom de övergripande observationerna med von Mises-spänning visade sig vara likartade. Analyser och resultat med von Mises-spänning ingår i de uppgifter som är kopplade till denna undersökning om läsaren vill ta del av dem. För varje simulerat fall observerades den maximala huvudspänningen (vid centrum av ett nätelement), kallad \({\sigma }_{MPS}\), i ett tryckområde från 0 till 5 kPa med en ökning på 0,25 kPa. Ett Python-skript utvecklades för att automatiskt söka efter det maximala värdet av \({\sigma }_{MPS}\), kallat \({\sigma }_{MPS,max}\) vid de olika trycken för varje simulerat fall. För varje modell beräknades det medelvärde som erhölls för varje simulerat fall (dvs. för alla fem uppsättningar av materialparametrar), och det behölls för efterföljande analys. Dessa värden kallas \({\sigma }_{MPS,\,max}^{avg}\) och \({\sigma }_{MPS,\,max}^{avg,\,n}\) för sjuka respektive normala modeller.
Värdena \({\sigma }_{MPS,\,max}^{avg}\) korrelerades, för tryck mellan 1 och 5 kPa (i steg om 0,25 kPa), med olika parametrar i påsens geometri (Fig. 1c) för att identifiera sambandet mellan påsens geometri och hög belastning: påsens halsbredd längs \(x\)-axeln \({D}_{x}\), påsens halsbredd längs \(z\)-axeln \({D}_{z}\), arean vid basen av påsens hals \({A}_{b}\), och påsens yta på lumensidan \({S}_{p}\). Pearsons korrelationskoefficient \({r}^{2}\) beräknades för varje parameter. Det tvåsidiga p-värdet beräknades också för att utvärdera korrelationens betydelse.
För övrigt beräknade vi trycket \({P}_{d}\) i de sjuka modellerna så att för ett givet tryck \({P}_{n}\) i den normala tjocktarmen har vi:
Följande tryckfaktor beräknades därefter som:
Denna faktor anger den procentandel med vilken trycket \({P}_{n}\)måste sänkas i den sjuka tjocktarmen så att \({\sigma }_{MPS,max}\)\) liknar \({\sigma }_{MPS,\,max}^{n}\), dvs.e maximala värdet av den maximala huvudspänningen i den sjuka modellen inte överstiger motsvarande värde i den normala modellen. Det reducerade trycket betecknas då som \({P}_{d}\). För varje simulerat fall minskades \({P}_{n}\) med 0,25 kPa tills sambandet från ekv. 6 uppfylldes, vilket gav motsvarande värde för \({P}_{d}\) och tryckfaktorn \(f\) beräknades. Detta gjordes för \({P}_{n}\) i intervallet 2 till 5 kPa (i steg om 1 kPa). Det medelvärde som erhölls för modellerna över alla simulerade fall beräknades och behölls för efterföljande analys.
Slutligt fastställde vi påverkanszonen för en påse, dvs. avståndet runt en påse där spänningen är förhöjd. För att åstadkomma detta upprättade vi för varje sjukdomsfall en karta över den maximala huvudspänningen på tjocktarmens lumen sidoyta (eftersom spänningen är högst på denna yta) längs maskknutpunkter som är belägna på de längsgående och cirkumferentiella banorna som går genom pouchens centrum. Utvecklingen av den maximala huvudspänningen längs dessa banor observerades som en funktion av avstånden \({d}_{l}\) och \({d}_{c}\) mätt i förhållande till pouchens centrum. Därefter beräknade vi de längsta avstånden \({d}_{l}^{inf}\) längs den längsgående banan och \({d}_{c}^{inf}\) längs den cirkumferentiella banan efter vilka den maximala huvudspänningen sjunker med mindre än 10 % av motsvarande värde i den normala modellen längs samma banor. Dessa avstånd utgjorde vårt mått på påverkanszonen. Eftersom de når ett platåvärde efter ett visst tryck korrelerades deras medelvärden över alla fall för varje modell vid 5 kPa med olika parametrar för påsens geometri, på samma sätt som \({\sigma }_{MPS,\,max}^{avg}\).
Alla beräkningar utfördes i programmeringsspråket Python genom att utveckla Jupyter-notebooks42. Data lagrades och manipulerades med Pandas-biblioteket i Python43. Grundläggande matematiska operationer gjordes med biblioteket NumPy44. Korrelationskoefficienten \({r}^{2}\) och det tvåsidiga p-värdet beräknades med funktionen ”stats” i Scipy-biblioteket45. Alla diagram realiserades med biblioteket Matplotlib46.