Modelo colónico
El software de elementos finitos Abaqus (versión 6.13, SIMULIA, Providence, RI) se utilizó para desarrollar los modelos y realizar las simulaciones asociadas. Todas las simulaciones se realizaron en una estación de trabajo Dell T7810 con dos procesadores Intel Xeon y 32 GBytes de memoria. En nuestro estudio anterior sobre el comportamiento del tejido porcino, establecimos un modelo constitutivo de material anisotrópico para las regiones descendente y espiral del colon porcino, basado en ensayos simultáneos de inflación y extensión15. Este modelo ha demostrado en varios estudios anteriores que captura adecuadamente el comportamiento anisotrópico del colon19,34,35,36. Los parámetros del material establecidos en nuestro estudio anterior se mantuvieron ya que, hasta donde sabemos, es el único estudio que modela el colon de un animal grande bajo inflación y extensión, los dos principales modos de deformación del colon. Además, se ha establecido que el colon porcino es relevante para el estudio de la diverticulosis33,37. Por lo tanto, el modelo de material y los parámetros seleccionados son opciones adecuadas para este estudio. Dado que la anatomía del colon espiral no es relevante para los humanos, los resultados de las muestras de colon descendente se utilizaron aquí para caracterizar el material del tejido del colon. En nuestro trabajo anterior hemos informado de que hay una variación insignificante en el diámetro localmente a lo largo de la longitud de un segmento de colon descendente porcino, por lo que el colon normal se modeló como un segmento cilíndrico (Fig. 1a). Los radios interior y exterior se fijaron en \({R}_{i}=11,5\) mm y \({R}_{o}=12,7\) mm, respectivamente, basándose en los valores medios medidos en nuestro trabajo anterior. La longitud axial se fijó en \(L=10\) cm de forma que es varios pliegues mayor que el radio. En concreto, se asumió que el tejido del colon es un material hiperelástico homogéneo, anisotrópico e incompresible cuyo comportamiento mecánico se describe mediante la siguiente función de energía de deformación anisotrópica \(\bar{W}:\)15,35
El primer término representa una respuesta neo-Hookean que caracteriza el comportamiento de los componentes no colágenos del tejido del colon. La cantidad \({I}_{1}={lambda }_{z}^{2}+{lambda }_{theta }^{2}+{frac{1}{{lambda }_{z}^{2}{lambda }_{theta }^{2}) representa la primera invariante del tensor de deformación, con \ {\lambda }_{z}\) y \ {\lambda }_{theta }\) referidos a los tramos axiales y circunferenciales respectivamente. Los términos siguientes caracterizan la contribución de las fibras de colágeno, siendo \({k}_{1}}) y \({k}_{2}}) una medida de la rigidez de las fibras. En concreto, el segundo término (con superíndice \(l\)) caracteriza el comportamiento de las fibras de colágeno alineadas a lo largo de la dirección longitudinal (axial), mientras que el tercer término (con superíndice \(s\)) caracteriza la respuesta de las fibras de colágeno dispersas en dos direcciones simétricas preferentes con respecto a la dirección circunferencial. La cantidad \({I}_{4}\) se denomina pseudoinvariante y caracteriza la respuesta mecánica de las fibras a lo largo de las direcciones preferenciales. Para el segundo término, \({I}_4}^{l}={lambda }_{z}^{2}\) ya que da cuenta de las fibras a lo largo de la dirección longitudinal. Para el tercer término, \ {I}_{4}^{s}} se expresa como:
Aquí, \(\pm {\gamma }^{s}\) indica el ángulo de las dos orientaciones simétricas de las fibras con respecto a la dirección circunferencial. En este trabajo se utilizaron cinco conjuntos de parámetros de material, derivados previamente para cinco muestras de colon descendente porcino, para los parámetros de material \({C}_{10},,{k}_{1}^{l},,{k}_{2}^{l},,{k}_{1}^{s},,{k}_{2}^{s},\) y \({\gamma }^{s}). Esta función de energía de deformación no está disponible directamente en Abaqus. Por lo tanto, se integró por primera vez mediante la implementación de una subrutina de usuario basada en Fortran llamada UANISOHYPER_INV, que permite definir el comportamiento del material hiperelástico anisotrópico utilizando la formulación invariante. Esta subrutina requería definir la primera y segunda derivadas de la función tensión-energía (Ecuación 1) con respecto a los invariantes escalares \({I}_{1}}) y \({I}_{4}}). Para verificar que la subrutina se implementó correctamente y validar el modelo computacional del tejido del colon, comparamos los valores de tensión estimados por el modelo de la tensión circunferencial isocórica \({bar{sigma}_{theta }^{m}) y la tensión axial isocórica \({bar{sigma}_{z}^{m}) en la superficie externa del colon contra los valores obtenidos analíticamente sobre la base de la función de energía de deformación y calculados como sigue:
El coeficiente de determinación \({R}^{2}\) se calculó de la siguiente manera para cuantificar la desviación entre los valores computacionales y los analíticos:
donde \ {A={bar{{{sigma}}_{z},\bar{{sigma}}_{theta}), y el subíndice «\(avg\)» indica la media de los valores analíticos sobre todos los puntos de datos \(n\). Un valor de \({R}^{2}\) cercano a 1 indica que se obtiene globalmente una buena correlación.
Modelo de divertículo
Para los modelos de colon enfermo, se modeló una estructura similar a un divertículo como la intersección de una esfera de diámetro \(D\) con la forma cilíndrica del colon normal a una altura \(H\) sobre la superficie. En este estudio sólo se consideró una bolsa para centrarse en el efecto de una sola bolsa sin las variabilidades adicionales que incorporarían las bolsas múltiples (número de bolsas, posiciones a lo largo del colon, posición relativa entre ellas, variabilidad del tamaño individual, etc.). Hasta donde sabemos, no hay ningún estudio que informe sistemáticamente sobre la geometría de los divertículos. Se suele mencionar que suelen tener entre 5 y 10 mm de diámetro, sin que se indique adecuadamente dónde se mide el diámetro38,39. Para cubrir una variedad de tamaños de bolsas, se consideraron tres valores en las simulaciones para la altura inicial de la bolsa \(H\)y el diámetro inicial de la bolsa \(D\): \(H=2\) mm, 4 mm y 6 mm, y \(D\) = 8 mm, 10 mm y 12 mm. Las combinaciones de estos valores de D y H dieron lugar a 9 modelos enfermos con diferentes tamaños de bolsa. Se incluyó un filete redondo de radio \({r}_{f}=2\) mm en la intersección entre el cilindro que representa el colon normal y la esfera que representa la bolsa para eliminar los bordes afilados que pueden no encontrarse en el tejido real y pueden alterar los valores de tensión. Se supuso que el grosor del tejido del divertículo era el mismo que el del colon (es decir, \({R}_{o}-{R}_{i}\)). Se ha intentado caracterizar la mecánica del colon en la diverticulosis23,30,40,41. Sin embargo, no se han establecido las propiedades del material específicamente para el tejido del divertículo (es decir, la caracterización de la bolsa del divertículo aislada). Por lo tanto, se consideraron propiedades similares a las del tejido normal para la sección de la bolsa.
Condiciones de contorno
Las condiciones de contorno se impusieron basándose en nuestro estudio anterior del colon porcino15. Observamos que los valores de presión luminal (presión en la pared interior del tejido con una presión exterior igual a cero) superiores a 1,5 kPa conducían a un daño permanente del tejido en condiciones pasivas. En las simulaciones computacionales, este valor se amplió para visualizar la tendencia de la evolución de los valores de tensión y se impuso una presión luminal de hasta 5 kPa en la superficie interna del colon, incluida la bolsa. La presión de la pared exterior se mantuvo igual a 0. Dado que estamos trabajando con un tejido de paredes finas (relación entre el diámetro interior y el grosor cercana a 20), la presión luminal podría asimilarse así como el delta entre la presión interior (por ejemplo, la presión impuesta dentro del conducto del tejido por la materia fecal o el gas intestinal) y la presión exterior (por ejemplo, la presión abdominal). Normalmente, se necesitaría una pequeña presión para abrir el tejido del colon que, de lo contrario, se colapsaría. Esa presión (~0,2 kPa) es, sin embargo, muy pequeña en comparación con el rango considerado aquí, por lo que se ignora. La presión se aplicó en incrementos cuasiestáticos. En nuestro estudio anterior se observó típicamente un estiramiento axial de alrededor del 10% in vivo en el colon. Por lo tanto, se impuso un estiramiento axial del 10% a lo largo de la dirección \ (Z\) antes de la presurización. Hemos encontrado valores de ángulo de apertura muy pequeños en el tejido del colon porcino, por lo que se despreciaron las tensiones circunferenciales residuales en el modelo. Tanto para los modelos normales como para los enfermos, utilizamos la simetría con respecto al plano \(YZ\) y al plano \(XY\) y resolvimos sólo una cuarta parte del modelo para lograr una mayor eficiencia, antes de reconstruir el modelo completo con fines de visualización.
Mezhing
Se realizó un seccionamiento geométrico apropiado en cada caso para asegurar que el mallado es posible utilizando elementos hexaédricos lineales. Se seleccionó la formulación de elementos híbridos para reforzar la incompresibilidad del material (es decir, el elemento C3D8H de la biblioteca Abaqus). Se realizó un estudio de convergencia de la malla para el modelo normal reduciendo el tamaño de la malla de referencia a partir del tamaño por defecto sugerido automáticamente por Abaqus (1,3 mm) hasta que la variación del valor máximo de las tensiones principales máximas entre dos mallas sucesivas fuera inferior al 10%. El tamaño de malla convergente se utilizó entonces también en los modelos enfermos, con algunos refinamientos adicionales alrededor de la bolsa. Se utilizó un total de 69608 elementos y 88350 nodos para el modelo normal en el tamaño de malla más fino retenido para el análisis (un cuarto de la geometría). Para los nueve modelos enfermos se utilizaron entre 69764 elementos y 80656 elementos (88605 y 102440 nodos respectivamente) (una cuarta parte de la geometría). Se especificó una orientación local del material utilizando la opción de formulación discreta disponible en Abaqus para definir adecuadamente la orientación de las fibras en cada elemento tal y como se requiere en la función de energía de deformación.
Análisis de resultados
Se consideraron un total de 10 modelos (uno normal y nueve enfermos con diferentes tamaños de bolsa). Cada modelo se simuló con 5 conjuntos diferentes de parámetros de material, lo que dio lugar a un total de 50 casos simulados. Se necesita una medida del nivel de tensión para comparar los distintos modelos. Dado que no estamos estableciendo un criterio de fallo, sino simplemente observando las variaciones bajo diferentes condiciones, cualquier medida como la tensión principal máxima o la tensión de von Mises sería adecuada. Los resultados se muestran aquí sólo en términos de la tensión principal máxima para evitar información redundante, ya que las observaciones generales con la tensión de von Mises resultaron ser similares. El análisis y los resultados con la tensión de von Mises se incluyen en los datos asociados a este estudio si el lector desea consultarlos. Para cada caso simulado, se observó la tensión principal máxima (en el centroide de un elemento de la malla), denominada \({\sigma }_{MPS}\), en un rango de presión de 0 a 5 kPa con un incremento de 0,25 kPa. Se desarrolló un script en Python para buscar automáticamente el valor máximo de \ {{sigma }_{MPS}}, designado como \ {{sigma }_{MPS,max}} en las distintas presiones para cada caso simulado. Para cada modelo, se calculó el valor medio obtenido para cada caso simulado (es decir, para los cinco conjuntos de parámetros del material), y se retuvo para el análisis posterior. Estos valores se denominan \ {{sigma}_{MPS,\}^{avg}\} y \ {{sigma}_{MPS,\}^{avg,\}\} para los modelos enfermos y normales, respectivamente.
Los valores \({\sigma }_{MPS,\}^{avg}\), se correlacionaron, para una presión entre 1 y 5 kPa (en incremento de 0,25 kPa), con diferentes parámetros de la geometría de la bolsa (Fig. 1c) para identificar la relación entre la geometría de la bolsa y la tensión elevada: la anchura del cuello de la bolsa a lo largo del eje \(x) \({D}_{x}\), la anchura del cuello de la bolsa a lo largo del eje \(z) \({D}_{z}\), el área en la base del cuello de la bolsa \({A}_{b}\), y la superficie de la bolsa en el lado del lumen \({S}_{p}\). Se calculó el coeficiente de correlación de Pearson \ {r}^{2}} para cada parámetro. También se calculó el valor p de dos colas para evaluar la importancia de la correlación.
Además, se calculó la presión \({P}_{d}\} en los modelos enfermos de forma que para una presión dada \({P}_{n}\}) en el colon normal tenemos:
El siguiente factor de presión se calculó como:
Este factor indica el % en el que debe reducirse la presión \({P}_{n}} en el colon enfermo para que \({sigma }_MPS,max}\) sea similar a \ {{sigma }_{MPS,\}^{n}\}, es decire valor máximo de la tensión principal máxima en el modelo enfermo no supera el valor correspondiente en el modelo normal. La presión reducida se designa entonces como \({P}_{d}\}). Específicamente, para cada caso simulado, \({P}_{n}) se redujo en 0,25 kPa hasta que se satisfizo la relación de la Ec. 6, proporcionando el valor correspondiente de \({P}_{d}\) y se calculó el factor de presión \(f\). Esto se hizo para \({P}_{n}} en el rango de 2 a 5 kPa (en incremento de 1 kPa). Se calculó el valor medio obtenido para los modelos en todos los casos simulados y se retuvo para el análisis posterior.
Por último, establecimos la zona de influencia de una bolsa, es decir, la distancia alrededor de una bolsa en la que se eleva la tensión. Para ello, establecimos para cada caso enfermo un mapa de la tensión principal máxima en la superficie lateral del lumen del colon (ya que la tensión es mayor en esta superficie) a lo largo de los nodos de malla situados en las trayectorias longitudinal y circunferencial que pasan por el centro de la bolsa. La evolución de la tensión principal máxima a lo largo de estas trayectorias se observó como una función de las distancias \ ({d}_{l}\) y \ ({d}_{c}\) medidas con respecto al centro de la bolsa. A continuación, calculamos las distancias más lejanas ({d}_{l}^{inf}} a lo largo de la trayectoria longitudinal y \({d}_{c}^{inf}}) a lo largo de la trayectoria circunferencial después de las cuales la tensión principal máxima cae dentro del 10% del valor correspondiente en el modelo normal a lo largo de las mismas trayectorias. Estas distancias constituyen nuestra medida de la zona de influencia. Observando que alcanzan un valor de meseta a partir de una determinada presión, sus valores medios sobre todos los casos para cada modelo a 5 kPa se correlacionaron con diferentes parámetros de la geometría de la bolsa, de forma similar a \({{sigma }_{MPS,\,max}^{avg}\).
Todos los cálculos se realizaron en el lenguaje de programación Python mediante el desarrollo de Jupyter notebooks42. Los datos fueron almacenados y manipulados con la librería Pandas de Python43. Las operaciones matemáticas básicas se realizaron con la librería NumPy44. El coeficiente de correlación \({r}^{2}\) y el valor p de dos colas se calcularon con la función «stats» de la biblioteca Scipy45. Todos los gráficos se realizaron con la biblioteca Matplotlib46.