Colon model
有限要素ソフトウェアAbaqus (version 6.13, SIMULIA, Providence, RI) はモデルの開発と関連シミュレーションを行うために使用されました。 すべてのシミュレーションは、2 つの Intel Xeon プロセッサと 32 GBytes のメモリを搭載した Dell Workstation T7810 で実行されました。 豚の組織挙動に関する我々の以前の研究では、豚の結腸の下行および螺旋領域について、膨張と伸展の同時試験に基づいて、異方性材料構成モデルを確立した15。 このモデルは、結腸の異方性挙動を適切に捉えることができることが、これまでのいくつかの研究で示されています19,34,35,36。 大腸の2大変形モードである膨張と伸展をモデル化した研究は、我々の知る限りこの研究が唯一であるため、前回の研究で確立した材料パラメータを維持した。 さらに、豚の大腸は憩室症の研究に適していることが立証されている33,37. このように、選択した材料モデルとパラメータは、本研究にとって適切な選択である。 らせん状結腸の解剖学的構造はヒトには当てはまらないため、ここでは下行結腸の試料から得られた結果を結腸組織材料の特性評価に使用しました。 我々は、豚の下行結腸の長さ方向に直径の局所的な変化はほとんどないことを以前の研究で報告しているため、正常結腸を円柱状のセグメントとしてモデル化した（図1a）。 内径と外径は、これまでの研究で得られた平均値から、それぞれ \({R}_{i}=11.5) mm と \({R}_{o}=12.7) mm とした。 また、軸方向の長さは半径の数倍となるように(L=10㎠)cmとした。 具体的には、大腸組織を均質・異方性・非圧縮性の超弾性材料とし、その力学的挙動を以下の異方性歪みエネルギー関数で記述することとした。\)15,35
第一項はNeo_BAR_Wを表し、第二項はKey_BAR_Wを表します。フック応答は、大腸組織の非コラーゲン成分の挙動を特徴づける。 このとき、量({I}_{1}={{lambda }_{z}^{2}+{{plambda }_{theta }^{2}+δfrac{1}{{lambda }_{z}^{2}{lambda }_{theta }^{2}})は第1項を表し、δfrac{1}は第2項を表す。 変形テンソルの不変量で、 \({ Θlambda }_{z}}) と \({ Θlambda }_{theta }} は軸方向と周方向の伸びを表します。 それぞれ 後続の項はコラーゲン繊維の寄与を特徴づけるもので、 \({k}_{1}) と \({k}_{2}) は繊維の硬さの尺度である。 具体的には、第2項（上付き記号 “li”）は長手方向（軸方向）に沿ったコラーゲン繊維の挙動を特徴付け、第3項（上付き記号 “s”）は円周方向に関して優先的に対称な2方向に分散したコラーゲン繊維の応答を特徴付ける。 この量は擬似不変量と呼ばれ、繊維の力学的応答を特徴づける。 第2項は、長手方向に沿った繊維を考慮するため、 \({I}_{4}^{l}={{lambda }_{z}^{2}}) とする。 第3項について、 \({I}_{4}^{s}} は次のように表される。
ここで、。 \を示し、円周方向に対して対称な2つの繊維配向角度を示しています。 材料パラメーターは、豚の下行結腸5サンプルについて過去に導出した5セットを使用し、材料パラメーター \({C}_{10},\,{k}_{1}^{l},{k}_{2}^{l},{k}_{1}^{s},\,{k}_{2}^{s}, \) と \({hatgamma }^{s}) の2つのパラメーターを使用することとしました。 このひずみエネルギー関数はAbaqusでは直接利用することができません。 そこで、FortranベースのユーザーサブルーチンであるUANISOHYPER_INVを実装することで、初めて統合され、不変式を使って異方性超弾性材料挙動を定義することが可能となりました。 このサブルーチンでは，ひずみエネルギー関数（式1）の1次導関数と2次導関数を，スカラー不変量({I}_{1})と({I}_{4})に関して定義する必要があった． 本サブルーチンが適切に実装されていることを確認し、大腸組織の計算モデルを検証するために、大腸外表面における等時性周方向応力度⊖({⊖bar{sigma }}_{theta }^{m})、等時性軸方向応力度⊖({⊖bar{sigma }}_{z}^{m}) のモデル推定値と歪エネルギー関数に基づいて解析的に得られた値との比較により下記の計算を実施しました。
計算値と解析値の偏差を定量化するために、決定係数({R}^{2})を次のように計算した。
where \(A={bar{sigma }}_{z},\,{bar{sigma }}_{theta }}であり、添え字の “avg}”は全データポイントの解析値の平均値を示しています。 3532>
憩室モデル
疾患大腸モデルでは、正常大腸の円筒形状に、表面から高さ(H)で直径(D)の球が交わる憩室様構造をモデル化した。 本研究では、複数のパウチから取り込まれる追加的な変動要因（パウチの数、結腸に沿った位置、互いの相対位置、個体サイズの変動など）を排除して、単一のパウチの効果に焦点を当てるために、1つのパウチのみを考慮した。 我々の知る限りでは、憩室の形状を系統的に報告した研究はない。 一般的に直径5〜10mmと言われているが、その直径がどこで測定されたものなのか、適切な指標はない38,39。 様々なサイズの袋をカバーするため、シミュレーションでは、初期袋の高さ(H)と初期袋の直径(D)を3つの値で検討した。 \また、袋の高さ(H)と直径(D)については、以下の3つの値を考慮してシミュレーションを行った。 これらのDとHの組み合わせにより、袋の大きさが異なる9種類の疾患モデルを作成した。 また，正常結腸を表す円柱とパウチを表す球の交点に半径 Δ({r}_{f}=2}mm の円形フィレットを入れ，実際の組織にはない鋭いエッジが応力値を 変える可能性があることを示した． また、憩室組織の厚さは大腸と同じとした（すなわち、 \({R}_{o}-{R}_{i}) ）。 これまでにも憩室症における大腸の力学的特性を明らかにする試みがなされてきた23,30,40,41. しかし、憩室組織に特化した材料特性は確立されていない（すなわち、孤立した憩室袋の特性）。 3532>
境界条件
境界条件は我々が以前に行った豚の大腸の研究15に基づいて設定した。 我々は、内腔圧（外圧をゼロとしたときの組織内壁の圧力）が1.5kPaを超えると、受動条件下で組織の永久的な損傷につながることを観察した。 計算機シミュレーションでは、応力値の変化の傾向を可視化するためにこの値をさらに押し上げ、5kPaまでの内腔圧を袋を含む大腸の内表面に課した。 外壁の圧力は0に保った。薄壁の組織（内径と厚さの比が20に近い）を扱っているので、内腔圧は内圧（例えば、糞便や腸内ガスによって組織の導管内に与えられた圧力）と外圧（例えば、腹圧）の差として同化される可能性がある。 一般に、そうでなければ崩壊してしまう結腸の組織を開くには、小さな圧力が必要であろう。 しかし、その圧力（〜0.2kPa）は、ここで検討した範囲に比べれば非常に小さいので、無視することにする。 圧力は準静的な増分で適用した。 私たちの以前の研究では、約10％の軸方向伸縮が結腸のin vivoで典型的に見られた。 そこで、加圧前に軸方向の伸縮を10%付与した。 また、豚の大腸組織では開口角の値が非常に小さく、残留する周方向応力を無視したモデルとなっています。 正常モデル、疾患モデルともに、healthy planeとhealthy planeに関して対称性を持たせ、効率を上げるためにモデルの1/4のみを解き、可視化のためにモデル全体を再構築した
Meshing
線形六面体要素でメッシングできるよう、各ケースで適切な形状切断を行った。 材料の非圧縮性を強制するためにハイブリッド要素定式化を選択しました（すなわち、Abaqusライブラリの要素C3D8H）。 通常モデルでは，Abaqus が自動的に提案するデフォルトのメッシュサイズ（1.3mm）か ら，連続する 2 つのメッシュ間での最大主応力の最大値の変動が 10%未満になるまで参照メッシュサイズを小さくして，メッシュ収束の調査を実施しました． その後，収束したメッシュサイズを疾患モデルにも使用し，袋の周囲をさらに細分化した． 正常なモデルでは，解析に使用した最も細かいメッシュサイズ（形状の4分の1）で，合計69608要素および88350節点が使用されました． 9つの疾患モデルには，69764要素と80656要素（それぞれ88605と102440ノード）が使用されました（ジオメトリの4分の1）． ひずみエネルギー関数で必要とされる各要素の繊維配向を適切に定義するため、Abaqusで利用可能な離散定式化オプションを使用して、局所的な材料配向を指定しました。 各モデルは5つの異なる材料パラメータのセットでシミュレーションされ、合計50のシミュレーションケースとなった。 各モデルを比較するためには、応力レベルの測定が必要である。 ここでは破壊基準を設定するのではなく，異なる条件下でのばらつきを調べるだけなので，最大主応力やフォンミーゼス応力など，どのような指標でも適切に用いることができます． フォンミーゼス応力による全体的な観察結果は同様であったため、ここでは情報の重複を避けるために、最大主応力の観点からのみ結果を示しています。 なお，von Mises 応力を用いた解析とその結果は，本研究の関連資料に含まれているので，そちらを参照されたい． 各シミュレーションケースにおいて，0～5kPa の圧力範囲で 0.25kPa 刻みで最大主応力（メッシュ要素のセントロイドにおける）を観測した（ \({sigma}_{MPS}）． また、Pythonスクリプトを開発し、各シミュレーションケースで圧力が変化したときの最大値( \({Θsigma}_{MPS,max})を自動的に検索するようにしました。 また，各モデルについて，各ケース（5 つの材料パ ラメータセットすべて）の平均値を算出し，以降の解析 に使用した． これらの値を、疾患モデル、正常モデルそれぞれについて、 \({ Θsigma }_{MPS,Θ,max}^{avg,Θ,n}}) と呼びます。
また、1～5kPa（0.25kPa刻み）の圧力に対して、その値であるⒶ({Θsigma }_{MPS,\,max}^{avg}} )と袋体形状の異なるパラメータとの相関を調べた（Fig.5）。 1c）では、パウチ形状と高応力の関係を明らかにするため、パウチネックの \(x) 軸の幅、 \(z) 軸の幅、パウチネックの基部の面積、ルーメン側のパウチの表面形状をパラメータとし、相関を求めた（図1）。 各パラメーターについてピアソン相関係数({r}^{2}}を算出した。 相関の有意性を評価するために、両側p値も計算した。
さらに、正常大腸の圧({P}_{n}}が与えられた場合、疾患モデルにおける圧({P}_{d}})が以下のようになるように計算した。
続いて以下の圧力係数が計算されました。
この係数は、病変大腸において、｛P｝_{n}の圧力を何％低減すれば、｛P}_{sigma }_{MPS.P}｝が得られるかを示しています。max}) が、 \({ {hatsigma }_{MPS,\,max}^{n}}) と同じような値であること。のように、疾患モデルにおける最大主応力の最大値が正常モデルにおける対応する値を超えないようにする。 また、減圧したときの圧力は、Ⓐ({P}_{d}Ⓑ)とする。 具体的には、各シミュレーションケースについて、式（6）の関係を満たすまで0.25kPaずつ減圧していき、対応する “P}_{d}”の値と圧力係数 “f}”を計算した。 これは、2～5kPaの範囲（1kPa刻み）で行われました。 3532>
最後に、パウチの影響領域、つまり応力が上昇するパウチの周囲の距離を設定しました。 そのために、各症例について、結腸の内腔側表面における最大主応力（この表面が最も応力が高いため）のマップを、袋の中心を通る縦方向および周方向の経路上にあるメッシュ状のノードに沿って作成した。 これらの経路に沿った最大主応力の変化を、袋の中心を基準として測定した距離(δ({d}_{l})とδ({d}_{c})の関数として観測した。 次に、縦方向の経路において、最大主応力が通常モデルにおける対応値の10%以内に低下する最も遠い距離(縦方向の経路)、および周方向の経路における距離(周方向の経路)を算出した。 これらの距離を影響領域とした。 ある圧力を超えるとプラトー値になることを利用して、5kPaにおける各モデルの全ケースでの平均値を、 \({sigma }_{MPS,\,max}^{avg}}) と同様に袋の形状の異なるパラメータと相関させた。
計算はすべてプログラミング言語Pythonで、Jupyter notebooksを開発した42。 データの保存と操作はPythonのPandasライブラリで行った43。 基本的な数学的演算はNumPyライブラリで行った44。 相関係数({r}^{2}}と両側p値は、Scipyライブラリのstats関数で計算した45。 プロットはすべて Matplotlib ライブラリで実現した46.
．