Colon model
Den finite element software Abaqus (version 6.13, SIMULIA, Providence, RI) blev brugt til at udvikle modellerne og udføre de tilhørende simuleringer. Alle simuleringer blev udført på en Dell Workstation T7810 med to Intel Xeon-processorer og 32 GByte hukommelse. I vores tidligere undersøgelse af svinevævets adfærd etablerede vi en anisotropisk materialekonstitutiv model for de nedadgående og spiralformede regioner af svinekolonet baseret på samtidige inflations- og udvidelsesforsøg15. Denne model har i flere tidligere undersøgelser vist sig at fange den anisotrope opførsel af tyktarmen på tilfredsstillende vis19,34,35,36. De materialeparametre, der blev fastlagt i vores tidligere undersøgelse, blev bibeholdt, da det, så vidt vi ved, er den eneste undersøgelse, der modellerer et stort dyrs tyktarm under inflation og udvidelse, de to vigtigste deformationsformer i tyktarmen. Desuden er det blevet fastslået, at svinekolon er relevant for undersøgelse af diverticulose33,37. Den valgte materialemodel og de valgte parametre er således passende valg for denne undersøgelse. Da spiralkolonets anatomi ikke er relevant for mennesker, blev resultaterne fra prøver fra det nedadgående kolon her anvendt til karakterisering af kolonvævsmaterialet. Vi har i vores tidligere arbejde rapporteret, at der er en ubetydelig variation i diameteren lokalt langs længden af et nedadgående kolon-segment fra svin, og derfor blev det normale kolon modelleret som et cylindrisk segment (fig. 1a). Den indre og ydre radius blev sat til henholdsvis \({R}_{i}=11,5\) mm og \({R}_{o}=12,7\) mm på grundlag af de gennemsnitlige værdier, der er målt i vores tidligere arbejde. Den aksiale længde blev sat til \(L=10\) cm, således at den er flere gange større end radius. Specifikt blev tyktarmsvævet antaget at være et homogent, anisotropt og inkompressibelt hyperelastisk materiale, hvis mekaniske opførsel beskrives ved følgende anisotropiske deformationsenergifunktion \(\(\bar{W}:\)15,35
Det første udtryk repræsenterer en neo-Hookean respons, der karakteriserer opførslen af de ikke-kollagene bestanddele af tyktarmsvævet. Mængden \({I}_{1}={{\lambda }_{z}^{2}+{\lambda }_{\theta }^{2}+\frac{1}{{{\lambda }_{z}^{2}{{\lambda }_{\theta }^{2}}}\) repræsenterer den første invariant af deformationstensoren med \({\lambda }_{z}\) og \({\lambda }_{\theta }\), der henviser til de aksiale og cirkumferentielle strækninger, henholdsvis. De efterfølgende termer karakteriserer bidraget fra kollagenfibrene, idet \({k}_{1}\) og \({k}_{2}\) er et mål for fiberstivhed. Mere specifikt karakteriserer det andet udtryk (med overskriften \(l\\)) kollagenfibrenes opførsel langs den langsgående (aksiale) retning, mens det tredje udtryk (med overskriften \(s\)) karakteriserer responsen fra kollagenfibrene spredt i to præferentielle symmetriske retninger i forhold til omkredsretningen. Mængden \({I}_{4}\) betegnes som pseudo-invariant og karakteriserer fibrenes mekaniske respons langs de præferentielle retninger. For det andet udtryk gælder \({I}_{4}^{l}}={\lambda }_{z}^{2}\), da det tager højde for fibrene langs længderetningen. For det tredje udtryk udtrykkes \({I}_{4}^{s}\) som:
Her, \(\pm {\gamma }^{s}\) angiver vinklen mellem de to symmetriske fiberorienteringer vinkel i forhold til omkredsretningen. Fem sæt materialeparametre, der tidligere er afledt for fem prøver af nedadgående tyktarm fra svin, blev anvendt i dette arbejde til materialeparametrene \({C}_{10},\,{k}_{1}^{l},\,{k}_{2}^{l},\,{k}_{1}^{s},\,{k}_{2}^{s},\) og \({\gamma }^{s}\). Denne strain-energifunktion er ikke direkte tilgængelig i Abaqus. Den blev derfor for første gang integreret ved at implementere en Fortran-baseret brugerunderrutine kaldet UANISOHYPER_INV, som gør det muligt at definere anisotropisk hyperelastisk materialeadfærd ved hjælp af den invariante formulering. Denne underrutine krævede, at man definerede første og anden afledte af strain-energifunktionen (Eq. 1) med hensyn til de skalare invarianter \({I}_{1}\) og \({I}_{4}\). For at verificere, at underprogrammet blev implementeret korrekt, og validere den beregningsmæssige model af tyktarmsvævet, sammenlignede vi de modelestimerede spændingsværdier for isokorisk omkredsspænding \({\bar{\sigma }}_{\theta }^{m}\) og isokorisk aksialspænding \({\bar{{\sigma }}_{z}^{m}\) på den ydre tyktarmsoverflade med værdier, der er opnået analytisk på grundlag af strain-energifunktionen og beregnet som følger:
Determinationskoefficienten \({R}^{2}\) blev beregnet som følger for at kvantificere afvigelsen mellem beregningsværdier og analytiske værdier:
hvor \(A={{\bar{\sigma }}_{z},\,{\bar{\sigma }}_{\theta }}\), og subscriptet “\(avg\)” angiver gennemsnittet af de analytiske værdier over alle \(n\) datapunkter. En værdi af \({R}^{2}\) tæt på 1 angiver, at der globalt set opnås en god korrelation.
Divertikulummodel
For modellerne af den syge tyktarm blev en divertikellignende struktur modelleret som skæringspunktet mellem en kugle med diameter \(D\) og den cylindriske form af den normale tyktarm i en højde \(H\) over overfladen. Der blev kun taget hensyn til én pose i denne undersøgelse for at fokusere på virkningen af en enkelt pose uden yderligere variabiliteter, der ville blive indarbejdet ved flere poser (antal poser, positioner langs tyktarmen, relativ position i forhold til hinanden, variabilitet i den individuelle størrelse osv.) ). Så vidt vi ved, findes der ingen undersøgelse, der systematisk rapporterer divertiklernes geometri. Det nævnes almindeligvis, at de typisk er mellem 5 og 10 mm i diameter uden ordentlig angivelse af, hvor diameteren er målt38,39. For at dække en række forskellige posestørrelser blev der i simuleringerne taget hensyn til tre værdier for den oprindelige posens højde \(H\) og den oprindelige posens diameter \(D\): \(H=2\) mm, 4 mm og 6 mm, og \(D\) = 8 mm, 10 mm og 12 mm. Kombinationer af disse værdier for D og H førte til 9 syge modeller med forskellige posestørrelser. Der blev medtaget en rund filet med radius \({r}_{f}=2\) mm i skæringspunktet mellem cylinderen, der repræsenterer den normale tyktarm, og kuglen, der repræsenterer posen, for at fjerne skarpe kanter, som ikke nødvendigvis forekommer i det faktiske væv og kan ændre spændingsværdierne. Tykkelsen af divertikelvævet blev antaget at være den samme som tykkelsen af tyktarmen (dvs. \({R}_{o}-{R}_{i}\)). Der har været gjort forsøg på at karakterisere kolonmekanikken ved divertikulose23,30,40,41. Der er imidlertid ikke blevet fastlagt nogen materialeegenskaber specifikt for divertikelvævet (dvs. karakterisering af den isolerede divertikelpose). Egenskaber svarende til det normale væv blev derfor overvejet for posesektionen.
Grænsebetingelser
Grænsebetingelserne blev pålagt på baggrund af vores tidligere undersøgelse af svinekolon15. Vi observerede, at luminalt tryk (tryk på vævets indervæg med et ydre tryk lig nul) værdier over 1,5 kPa førte til permanent beskadigelse af vævet under passive forhold. I de beregningsmæssige simuleringer blev denne værdi skubbet yderligere for at visualisere tendensen i udviklingen af stressværdierne, og der blev pålagt et luminalt tryk på op til 5 kPa på den indre overflade af tyktarmen, herunder posen. Det ydre vægtryk blev holdt lig med 0. Da vi arbejder med et tyndvægget væv (forholdet mellem indre diameter og tykkelse er tæt på 20), kan det luminale tryk således sidestilles med deltaet mellem det indre tryk (f.eks. det tryk, der pålægges inde i vævets kanal af afføring eller tarmgas) og det ydre tryk (f.eks. abdominalt tryk). Typisk vil der være behov for et lille tryk for at åbne tyktarmsvævet, som ellers ville kollapse. Dette tryk (~0,2 kPa) er imidlertid meget lille i forhold til det område, der er undersøgt her, så der ses bort fra det. Trykket blev påført i kvasi-statiske trin. En aksial strækning på ca. 10 % blev typisk set in vivo i tyktarmen i vores tidligere undersøgelse. Der blev således pålagt en 10 % aksial strækning langs \(Z\)-retningen før trykpåvirkning. Vi har fundet meget små åbningsvinkelværdier i svinekolonvævet, og de resterende cirkumferentielle spændinger blev derfor negligeret i modellen. For både normale og syge modeller anvendte vi symmetri med hensyn til \(YZ\)-planet og \(XY\)-planet og løste kun en fjerdedel af modellen for at opnå højere effektivitet, før vi rekonstruerede hele modellen til visualiseringsformål.
Meshing
Der blev foretaget passende geometriudsnit i hvert enkelt tilfælde for at sikre, at meshing er mulig ved hjælp af lineære hexaedriske elementer. Der blev valgt en hybrid elementformulering for at håndhæve materialets inkompressibilitet (dvs. element C3D8H fra Abaqus-biblioteket). Der blev foretaget en undersøgelse af netkonvergens for den normale model ved at reducere referencestørrelsen af nettet fra den standardstørrelse, der automatisk foreslås af Abaqus (1,3 mm), indtil variationen i den maksimale værdi af de maksimale hovedspændinger mellem to på hinanden følgende net var mindre end 10 %. Den konvergerende maskestørrelse blev derefter også anvendt i de syge modeller med nogle yderligere finjusteringer omkring posen. Der blev anvendt i alt 69608 elementer og 88350 knuder til den normale model ved den fineste maskestørrelse, der blev bibeholdt til analyse (en fjerdedel af geometrien). Der blev anvendt mellem 69764 elementer og 80656 elementer (henholdsvis 88605 og 102440 knuder) for de ni syge modeller (en fjerdedel af geometrien). En lokal materialeorientering blev specificeret ved hjælp af den diskrete formuleringsmulighed, der er tilgængelig i Abaqus, for at definere fiberorienteringen i hvert element korrekt, som det kræves i strain-energifunktionen.
Resultater analyse
I alt 10 modeller (en normal og ni syge med forskellige posestørrelser) blev overvejet. Hver model blev simuleret med 5 forskellige sæt materialeparametre, hvilket resulterede i i alt 50 simulerede tilfælde. Et mål for belastningsniveauet er nødvendigt for at kunne sammenligne på tværs af modellerne. Da vi ikke fastlægger et fejlkriterium, men blot ser på variationer under forskellige betingelser, vil ethvert mål som f.eks. maksimal hovedspænding eller von Mises-spænding være tilstrækkeligt. Resultaterne er her kun vist med hensyn til den maksimale hovedspænding for at undgå overflødig information, da de samlede observationer med von Mises-spænding viste sig at være ens. Analyser og resultater med von Mises-spænding er inkluderet i de data, der er knyttet til denne undersøgelse, hvis læseren ønsker at konsultere dem. For hvert simuleret tilfælde blev den maksimale hovedspænding (ved centroiden af et netelement), benævnt \({\sigma }_{MPS}\), observeret i et trykområde fra 0 til 5 kPa med en stigning på 0,25 kPa. Der blev udviklet et Python-script til automatisk at søge efter den maksimale værdi af \({\sigma }_{MPS}\), der betegnes \({\sigma }_{MPS,max}\) ved de forskellige tryk for hvert simuleret tilfælde. For hver model blev den gennemsnitlige værdi, der blev opnået for hvert simuleret tilfælde (dvs. for alle fem sæt materialeparametre), beregnet og anvendt til efterfølgende analyse. Disse værdier betegnes \({\sigma }_{MPS,\,max}^{avg}\) og \({\sigma }_{MPS,\,max}^{avg,\,n}\) for henholdsvis syge og normale modeller.
Værdierne \({\sigma }_{MPS,\,max}^{avg}\) blev for tryk mellem 1 og 5 kPa (i intervaller på 0,25 kPa) korreleret med forskellige parametre for posens geometri (fig. 1c) for at identificere sammenhængen mellem posens geometri og høj stress: posens halsbredde langs \(x\)-aksen \({D}_{x}\), posens halsbredde langs \(z\)-aksen \({D}_{z}\), arealet ved posens halsbasis \({A}_{b}\) og posens overflade på lumen-siden \({S}_{p}\). Pearson-korrelationskoefficienten \({r}^{2}\) blev beregnet for hver parameter. Den dobbeltsidede p-værdi blev også beregnet for at vurdere korrelationens signifikans.
Dertil kommer, at vi beregnede trykket \({P}_{d}\) i de syge modeller således, at for et givet tryk \({P}_{n}\) i den normale tyktarm har vi:
Den følgende trykfaktor blev efterfølgende beregnet som:
Denne faktor angiver den %, hvormed trykket \({P}_{n}}\)skal reduceres i den syge tyktarm, således at \({\sigma }_{MPS,max}}\) svarer til \({{\sigma }_{MPS,\,max}^{n}\), dvs.dvs. at den maksimale værdi af den maksimale hovedspænding i den syge model ikke overstiger den tilsvarende værdi i den normale model. Det reducerede tryk betegnes derefter som \({P}_{d}\). Konkret blev \({P}_{n}\) for hvert simuleret tilfælde reduceret med 0,25 kPa, indtil relationen fra Eq. 6 var opfyldt, hvilket gav den tilsvarende værdi af \({P}_{d}\), og trykfaktoren \(f\) blev beregnet. Dette blev gjort for \({P}_{n}\) i intervallet 2 til 5 kPa (i intervaller på 1 kPa). Den gennemsnitlige værdi, der blev opnået for modellerne over alle de simulerede tilfælde, blev beregnet og bibeholdt til den efterfølgende analyse.
Slutteligt fastlagde vi indflydelseszonen for en pose, dvs. afstanden omkring en pose, hvor spændingen er forhøjet. For at opnå dette etablerede vi for hvert sygdomstilfælde et kort over den maksimale hovedspænding på lumen-sidefladen af tyktarmen (da spændingen er størst på denne overflade) langs meshing-knuder placeret på de langsgående og cirkumferentielle baner, der går gennem centrum af posen. Udviklingen af den maksimale hovedspænding langs disse baner blev observeret som en funktion af afstandene \({d}_{l}\) og \({d}_{c}\) målt i forhold til posens centrum. Derefter beregnede vi de længste afstande \({d}_{l}^{inf}\) langs den langsgående bane og \({d}_{c}^{inf}\) langs den cirkumferentielle bane, hvorefter den maksimale hovedspænding falder inden for 10 % af den tilsvarende værdi i den normale model langs de samme baner. Disse afstande udgjorde vores mål for indflydelseszonen. Da de når en plateauværdi efter et vist tryk, blev deres gennemsnitsværdier over alle tilfælde for hver model ved 5 kPa korreleret med forskellige parametre for posens geometri på samme måde som \({\sigma }_{MPS,\,max}^{avg}\).
Alle beregninger blev udført i programmeringssproget Python ved at udvikle Jupyter notebooks42. Data blev gemt og manipuleret med Pandas-biblioteket i Python43. Grundlæggende matematiske operationer blev udført med NumPy-biblioteket44. Korrelationskoefficienten \({r}^{2}\) og den tosidede p-værdi blev beregnet med “stats”-funktionen i Scipy-biblioteket45. Alle plots blev realiseret med Matplotlib-biblioteket46.