Analyse computationnelle de la contrainte mécanique dans la diverticulose colique

Modèle de côlon

Le logiciel d’éléments finis Abaqus (version 6.13, SIMULIA, Providence, RI) a été utilisé pour développer les modèles et effectuer les simulations associées. Toutes les simulations ont été effectuées sur une station de travail Dell T7810 avec deux processeurs Intel Xeon et 32 Go de mémoire. Dans notre étude précédente du comportement des tissus porcins, nous avons établi un modèle constitutif de matériau anisotrope pour les régions descendante et spirale du côlon porcin, basé sur des tests simultanés d’inflation et d’extension15. Ce modèle a été démontré dans plusieurs études précédentes pour capturer de manière adéquate le comportement anisotrope du côlon19,34,35,36. Les paramètres matériels établis dans notre étude précédente ont été conservés car, à notre connaissance, il s’agit de la seule étude modélisant un côlon de grand animal sous inflation et extension, les deux principaux modes de déformation du côlon. De plus, il a été établi que le côlon porcin est pertinent pour l’étude de la diverticulose33,37. Le modèle de matériau et les paramètres sélectionnés sont donc des choix adéquats pour cette étude. Comme l’anatomie du côlon spiralé n’est pas pertinente pour les humains, les résultats des échantillons du côlon descendant ont été utilisés ici pour caractériser le matériau du tissu du côlon. Nous avons rapporté dans notre travail précédent qu’il y a une variation négligeable du diamètre localement le long de la longueur d’un segment de côlon descendant de porc, ainsi le côlon normal a été modélisé comme un segment cylindrique (Fig. 1a). Les rayons interne et externe ont été fixés respectivement à \({R}_{i}=11,5\) mm et \({R}_{o}=12,7\) mm, sur la base des valeurs moyennes mesurées lors de nos précédents travaux. La longueur axiale a été fixée à \(L=10\) cm de sorte qu’elle soit plusieurs fois supérieure au rayon. Plus précisément, le tissu du côlon a été supposé être un matériau hyperélastique homogène, anisotrope et incompressible dont le comportement mécanique est décrit par la fonction d’énergie de déformation anisotrope suivante : \(\bar{W} :\)15,35

$$\bar{W}={C}_{10}({I}_{1}-3)+\frac{{k}_{1}^{l}}{k}_{2}^{l}}+\frac{{k}_{1}^{s}{k}_{2}^{s}}$
(1)

Le premier terme représente une réponse néo-Hookean caractérisant le comportement des constituants non collagéniques du tissu du côlon. La quantité \({I}_{1}={\lambda }_{z}^{2}+{\lambda }_{\theta }^{2}+\frac{1}{{\lambda }_{z}^{2}{\lambda }_{\theta }^{2}}\) représente le premier invariant du tenseur de déformation avec \({\lambda }_{z}\) et \({\lambda }_{\theta }\) faisant référence aux étirements axial et circonférentiel, respectivement. Les termes suivants caractérisent la contribution des fibres de collagène, \({k}_{1}\) et \({k}_{2}\) étant une mesure de la rigidité des fibres. Plus précisément, le deuxième terme (avec l’exposant \(l\)) caractérise le comportement des fibres de collagène alignées dans la direction longitudinale (axiale), tandis que le troisième terme (avec l’exposant \(s\)) caractérise la réponse des fibres de collagène dispersées dans deux directions symétriques préférentielles par rapport à la direction circonférentielle. La quantité \({I}_{4}\) est appelée pseudo-invariante et caractérise la réponse mécanique des fibres le long des directions préférentielles. Pour le deuxième terme, \({I}_{4}^{l}={\lambda }_{z}^{2}\) car il tient compte des fibres le long de la direction longitudinale. Pour le troisième terme, \({I}_{4}^{s}\) s’exprime comme suit :

${I}_{4}^{s}={\lambda }_{\theta }^{2}{\cos }^{2}{\gamma }^{s}+{\lambda }_{z}^{2}{\sin }^{2}{\gamma }$
(2)

Ici, \(\N{2}{{gamma }^{s}}) indique l’angle des deux orientations symétriques des fibres par rapport à la direction circonférentielle. Cinq ensembles de paramètres de matériau, précédemment dérivés pour cinq échantillons de côlon descendant de porc, ont été utilisés dans ce travail pour les paramètres de matériau \({C}_{10},\,{k}_{1}^{l},\,{k}_{2}^{l},\,{k}_{1}^{s},\,{k}_{2}^{s},\) et \({\gamma }^{s}\). Cette fonction d’énergie de déformation n’est pas directement disponible dans Abaqus. Elle a donc été intégrée pour la première fois en implémentant une sous-routine utilisateur basée sur Fortran appelée UANISOHYPER_INV, qui permet de définir le comportement hyperélastique anisotrope du matériau en utilisant la formulation invariante. Ce sous-programme a nécessité la définition des dérivées premières et secondes de la fonction déformation-énergie (Eq. 1) par rapport aux invariants scalaires \({I}_{1}\) et \({I}_{4}\). Pour vérifier que le sous-programme a été mis en œuvre correctement et valider le modèle de calcul du tissu du côlon, nous avons comparé les valeurs de contrainte estimées par le modèle de la contrainte circonférentielle isochore \({\bar{\sigma }}_{\theta }^{m}\) et de la contrainte axiale isochore \({\bar{\sigma }}_{z}^{m}\) sur la surface externe du côlon aux valeurs obtenues analytiquement sur la base de la fonction d’énergie de déformation et calculées comme suit :

${\bar{\sigma }}_{\theta }^{a}={\lambda }_{\theta }\frac{\partial \bar{W}{\partial {\lambda }_{\theta }}. }_{\theta }}$$
(3)

${bar{\sigma }}_{z}^{a}={\lambda }_{z}\frac{\partial \bar{W}}{\partial {\lambda }_{z}}$
(4)

Le coefficient de détermination \({R}^{2}\) a été calculé comme suit pour quantifier l’écart entre les valeurs calculées et les valeurs analytiques :

$${R}^{2}(A)=1-\frac{{\sum }_{q=1}^{n}{({A}_{q}^{m}-{A}_{q}^{a})}^{2}}{{\sum }_{q=1}^{n}{({A}_{avg}^{a}-{A}_{q}^{a})}^{2}}$
(5)

où \(A={bar{\sigma }}_{z},\,{\bar{\sigma }}_{\theta }\), et l’indice « \(avg\) » indique la moyenne des valeurs analytiques sur tous les points de données \(n\). Une valeur de \({R}^{2}\) proche de 1 indique qu’une bonne corrélation est globalement obtenue.

Modèle de diverticule

Pour les modèles de colon malade, une structure de type diverticule a été modélisée comme l’intersection d’une sphère de diamètre \(D\) avec la forme cylindrique du colon normal à une hauteur \(H\) au-dessus de la surface. Une seule poche a été prise en compte dans cette étude afin de se concentrer sur l’effet d’une seule poche sans les variabilités supplémentaires qui seraient incorporées par des poches multiples (nombre de poches, positions le long du côlon, position relative les unes par rapport aux autres, variabilité de la taille des individus, etc.) A notre connaissance, il n’existe pas d’étude rapportant systématiquement la géométrie des diverticules. Il est couramment mentionné qu’ils ont typiquement un diamètre de 5 à 10 mm sans indication précise de l’endroit où le diamètre est mesuré38,39. Pour couvrir une variété de tailles de poches, trois valeurs ont été considérées dans les simulations pour la hauteur initiale de la poche \(H\) et le diamètre initial de la poche \(D\) : \(H=2\) mm, 4 mm et 6 mm, et \(D\) = 8 mm, 10 mm et 12 mm. Les combinaisons de ces valeurs de D et H ont permis d’obtenir 9 modèles malades avec différentes tailles de poches. Un congé rond de rayon \({r}_{f}=2\) mm a été inclus à l’intersection entre le cylindre représentant le côlon normal et la sphère représentant la poche afin d’éliminer les arêtes vives qui ne sont pas forcément rencontrées dans le tissu réel et qui peuvent modifier les valeurs de contrainte. L’épaisseur du tissu du diverticule a été supposée être la même que celle du côlon (c’est-à-dire \({R}_{o}-{R}_{i}\)). Des tentatives ont été faites pour caractériser la mécanique du côlon dans la diverticulose23,30,40,41. Aucune propriété matérielle, cependant, n’a été établie spécifiquement pour le tissu du diverticule (c’est-à-dire la caractérisation de la poche diverticulaire isolée). Des propriétés similaires au tissu normal ont donc été considérées pour la section de la poche.

Conditions limites

Les conditions limites ont été imposées sur la base de notre étude précédente du côlon porcin15. Nous avons observé que les valeurs de pression luminale (pression sur la paroi interne du tissu avec une pression externe égale à zéro) supérieures à 1,5 kPa entraînaient des dommages permanents au tissu dans des conditions passives. Dans les simulations informatiques, cette valeur a été poussée plus loin pour visualiser la tendance de l’évolution des valeurs de contrainte et une pression luminale allant jusqu’à 5 kPa a été imposée sur la surface interne du côlon, y compris la poche. La pression de la paroi extérieure a été maintenue égale à 0. Comme nous travaillons avec un tissu à paroi mince (rapport diamètre intérieur/épaisseur proche de 20), la pression luminale pourrait donc être assimilée au delta entre la pression intérieure (par exemple, la pression imposée à l’intérieur du conduit du tissu par les matières fécales ou les gaz intestinaux) et la pression extérieure (par exemple, la pression abdominale). Typiquement, une petite pression est nécessaire pour ouvrir le tissu du côlon qui, sinon, s’effondrerait. Cette pression (~0,2 kPa) est cependant très faible par rapport à la gamme considérée ici, elle est donc ignorée. La pression a été appliquée par incréments quasi-statiques. Un étirement axial d’environ 10% a été typiquement observé in vivo dans le côlon lors de notre étude précédente. Ainsi, un étirement axial de 10 % a été imposé dans la direction \(Z\) avant la pressurisation. Nous avons trouvé de très petites valeurs d’angle d’ouverture dans le tissu du côlon porcin, et les contraintes circonférentielles résiduelles ont donc été négligées dans le modèle. Pour les modèles normaux et malades, nous avons utilisé la symétrie par rapport au plan \(YZ\) et au plan \(XY\) et nous n’avons résolu qu’un quart du modèle pour obtenir une plus grande efficacité, avant de reconstruire le modèle entier à des fins de visualisation.

Maillage

Un sectionnement approprié de la géométrie a été effectué dans chaque cas pour s’assurer que le maillage est possible en utilisant des éléments hexaédriques linéaires. Une formulation d’éléments hybrides a été sélectionnée pour renforcer l’incompressibilité des matériaux (c’est-à-dire l’élément C3D8H de la bibliothèque Abaqus). Une étude de convergence du maillage a été menée pour le modèle normal en réduisant la taille du maillage de référence à partir de la taille par défaut proposée automatiquement par Abaqus (1,3 mm) jusqu’à ce que la variation de la valeur maximale des contraintes principales maximales entre deux maillages successifs soit inférieure à 10 %. Le maillage convergent a ensuite été également utilisé dans les modèles malades, avec quelques raffinements supplémentaires autour de la poche. Un total de 69608 éléments et 88350 nœuds ont été utilisés pour le modèle normal au maillage le plus fin retenu pour l’analyse (quart de la géométrie). Entre 69764 éléments et 80656 éléments (88605 et 102440 noeuds respectivement) ont été utilisés pour les neuf modèles malades (quart de la géométrie). Une orientation locale du matériau a été spécifiée en utilisant l’option de formulation discrète disponible dans Abaqus pour définir correctement l’orientation des fibres dans chaque élément comme requis dans la fonction d’énergie de déformation.

Analyse des résultats

Un total de 10 modèles (un normal et neuf malades avec différentes tailles de poches) ont été considérés. Chaque modèle a été simulé avec 5 ensembles différents de paramètres de matériaux, ce qui a donné un total de 50 cas simulés. Une mesure du niveau de contrainte est nécessaire pour comparer les différents modèles. Puisque nous n’établissons pas de critère de défaillance mais que nous examinons simplement les variations dans différentes conditions, toute mesure telle que la contrainte principale maximale ou la contrainte de von Mises serait adéquate. Les résultats ne sont présentés ici qu’en termes de contrainte principale maximale afin d’éviter les informations redondantes puisque les observations globales avec la contrainte de von Mises se sont avérées similaires. L’analyse et les résultats avec la contrainte de von Mises sont inclus dans les données associées à cette étude si le lecteur souhaite les consulter. Pour chaque cas simulé, la contrainte principale maximale (au centroïde d’un élément de maillage), appelée \({\sigma }_{MPS}\), a été observée dans une plage de pression de 0 à 5 kPa avec un incrément de 0,25 kPa. Un script Python a été développé pour rechercher automatiquement la valeur maximale de \({\sigma }_{MPS}\), désignée par \({\sigma }_{MPS,max}\) aux différentes pressions pour chaque cas simulé. Pour chaque modèle, la valeur moyenne obtenue pour chaque cas simulé (c’est-à-dire pour les cinq ensembles de paramètres de matériaux), a été calculée et retenue pour l’analyse ultérieure. Ces valeurs sont appelées \({\sigma }_{MPS,\,max}^{avg}\) et \({\sigma }_{MPS,\,max}^{avg,\,n}\) pour les modèles malades et normaux, respectivement.

Les valeurs \({\sigma }_{MPS,\,max}^{avg}\), ont été corrélées, pour une pression comprise entre 1 et 5 kPa (par incrément de 0,25 kPa), à différents paramètres de la géométrie de la poche (Fig. 1c) afin d’identifier la relation entre la géométrie de la poche et les contraintes élevées : la largeur du col de la poche le long de l’axe \(x\) \({D}_{x}\), la largeur du col de la poche le long de l’axe \(z\) \({D}_{z}\), la surface à la base du col de la poche \({A}_{b}\), et la surface de la poche du côté de la lumière \({S}_{p}\). Le coefficient de corrélation de Pearson \({r}^{2}\) a été calculé pour chaque paramètre. La valeur p bilatérale a également été calculée pour évaluer la signification de la corrélation.

De plus, nous avons calculé la pression \({P}_{d}\) dans les modèles malades de telle sorte que pour une pression donnée \({P}_{n}\) dans le colon normal, nous avons :

$$\max ({\sigma }_{MPS,max}({P}_{d}))\le \,\max ({\sigma }_{MPS,\,max}^{n}({P}_{n}))$$
(6)

Le facteur de pression suivant a ensuite été calculé comme :

$$f=\left(\frac{{P}_{n}-{P}_{d}}{P}_{n}}\right)\times 100$$
(7)

Ce facteur indique le % par lequel la pression \({P}_{n}\)doit être réduite dans le colon malade pour que \({\sigma }_{MPS,max}\) soit similaire à \({\sigma }_{MPS,\,max}^{n}\), c’est-à-dire que la valeur maximale de la contrainte principale max.e la valeur maximale de la contrainte principale maximale dans le modèle malade ne dépasse pas la valeur correspondante dans le modèle normal. La pression réduite est alors désignée par \({P}_{d}\\). Plus précisément, pour chaque cas simulé, \({P}_{n}\) a été réduite de 0,25 kPa jusqu’à ce que la relation de l’équation 6 soit satisfaite, fournissant la valeur correspondante de \({P}_{d}\) et le facteur de pression \(f\) a été calculé. Ceci a été fait pour \({P}_{n}\) dans la gamme de 2 à 5 kPa (par incrément de 1 kPa). La valeur moyenne obtenue pour les modèles sur l’ensemble des cas simulés a été calculée et retenue pour l’analyse ultérieure.

Enfin, nous avons établi la zone d’influence d’une poche, c’est-à-dire la distance autour d’une poche où la contrainte est élevée. Pour ce faire, nous avons établi pour chaque cas malade une carte de la contrainte principale maximale sur la surface latérale de la lumière du côlon (puisque la contrainte est la plus élevée sur cette surface) le long des nœuds de maillage situés sur les chemins longitudinal et circonférentiel passant par le centre de la poche. L’évolution de la contrainte principale maximale le long de ces chemins a été observée en fonction des distances \({d}_{l}\) et \({d}_{c}\) mesurées par rapport au centre de la poche. Ensuite, nous avons calculé les distances les plus éloignées \({d}_{l}^{inf}\) le long du trajet longitudinal et \({d}_{c}^{inf}\) le long du trajet circonférentiel après lesquelles la contrainte principale maximale chute à moins de 10% de la valeur correspondante dans le modèle normal le long des mêmes trajets. Ces distances constituent notre mesure de la zone d’influence. Observant qu’ils atteignent une valeur plateau après une certaine pression, leurs valeurs moyennes sur l’ensemble des cas pour chaque modèle à 5 kPa ont été corrélées à différents paramètres de la géométrie de la poche, de manière similaire à \({\sigma }_{MPS,\,max}^{avg}\).

Tous les calculs ont été effectués en langage de programmation Python en développant des carnets Jupyter42. Les données ont été stockées et manipulées avec la bibliothèque Pandas de Python43. Les opérations mathématiques de base ont été effectuées avec la bibliothèque NumPy44. Le coefficient de corrélation \({r}^{2}\) et la valeur p bilatérale ont été calculés avec la fonction « stats » de la bibliothèque Scipy45. Tous les graphiques ont été réalisés avec la bibliothèque Matplotlib46.

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.