ヨーの構成モデルを使った層状血管

はじめに比較
解析解近似用語集
3層モデル参考文献
均質化モデル

はじめに

このモデルは生理的条件下での冠動脈の応力状態を分析する．モデルは実験手順から抽出したモデルパラメータを持つヨーの超弾性材料モデルを利用する．

冠動脈は心臓に血液を供給する大きい血管である．冠動脈に関連する疾病の中には，ある生命体を傷つける可能性のあるものもある．一般的な心臓疾患は冠動脈壁にたまったプラークによって生じる．また，冠動脈瘤は危険な病状を呈する．これは冠動脈拡張と定義され，壁が正常な直径を超えている状態である．これらの疾病は動脈壁に大きく関連している．動脈壁は人間が生きている間，常に変化するアクティブな構造である．これは生物的・化学的・力学的な環境を感知し，恒常性と心細胞への最適な血液供給を保証するものである[Mishani, 2021]．

1.gif

この図は，冠血管の解剖学的構造から，関心のある力学特性を維持する簡約化モデルまでの複雑性低減を示している．2D構造は3層と対称条件（端に平行な方向にのみスライドする）でできており，これで円環全体を表している．

動脈壁の構造は複雑である．位相的な見地からすると，動脈壁は3つの異なる層からできている．内側の内膜は上皮組織と弾性繊維を持つ結合組織でできている．真ん中の層は中膜と言い，主に滑らかな筋細胞でできており，局所的および中心的な刺激による動脈直径を規制することができる．外側の層は外膜と言い，コラーゲンと弾性繊維を持つ結合組織でできている．これは動脈を周囲の構造と繋ぐ．

ここで3層血管と均質の血管を比較してみよう．3層は内膜，中膜，外膜を表す．それぞれ文献から，異なる力学的特性を持つ．形状は主となる左の冠動脈の解剖学的構造を表す直径を持つ3つの同心円筒に簡約される[Gholipour, 2018]．さらに血管特性を一つの均質層構造にする2つの異なる手法も示す．

このモデルでは，次の4つのモデリングアプローチを示す．1つ目は簡単な分析モデルである．2つ目のモデルでは，内膜，中膜，外膜を表す各層に異なる材料パラメータを持つ3層を使う．後の2つのモデルは簡約された均質化を表す．その一つは内膜の材料パラメータを使うというものであり，もう一つは3つの材料層すべての面積平均を使うというものである．

多層構造モデルと単層モデルの差を示すために，一軸引張試験の分析を行う．1つの方向に力を加えて実行すると，，λ2=λ3となる．

非圧縮性を強制するために，λ1λ2λ3=1を設定すると，最初のひずみの不変量 I1は以下のように一軸引張で表すことができる：

第一の主方向上のコーシー応力は以下になる：

圧力：

以下になる：

それぞれの層に対して別々の材料パラメータを指定する：

材料構成モデルを二次元のチャート（ひずみvs応力）としてプロットすることができる．生物組織では，kPaを使うと便利なことが多く，数量を再スケールするために倍率を導入する．

数量を再スケールするのに便利な倍率を導入する：
一軸応力をkPaで計算するヘルパー関数を定義する：
3層の一軸応力を可視化する：

解析解近似

この力学的問題は，円筒または球が内部圧力pの対象となるという特徴的な工学問題に還元することができる．これは圧力容器と呼ばれ，多様な分野で応用されている．

7.gif

この図は長さ l，厚さ t，内部圧力 の理想管を示している．周方向応力 および半径方向応力 を知ることは可能である．管の両端は開いているので，軸応力は除外することができる．

ここで，膨張して厚くなった壁の管の解析解を上限として提示する．これは，応力状態とFEM解を比較するために使うことができる．生体組織では，ひずみ仮定は失敗することが多いが，この還元された解析解と数値解のいくつかの違いを導く．

最初と2番目の主応力方向を示す．最初の主応力 は円周方向である．2番目の主応力 は半径方向である．

膨張した管では，内部圧力のため半径方向応力は少なくとも内側でこの境界圧力に関連している．また，外部の力または圧力がないため，半径方向応力は外側では0となる傾向がある．

さらにこの円筒はという比率を示しており，ラメ方程式を使って計算することのできる厚い壁の管となる．半径方向応力と接線応力を以下のように表すことができる：

ここで と は

および

である．は内部圧力，は外部圧力である．

軸力または軸制約条件がないため，軸応力はゼロに近い．

壁の厚い膨張した円筒の解析解：
形状の数量を導入する：
解析解を可視化する：

のプロットから分かるように，半径方向応力 は境界条件に関連している．半径方向応力 は血管の内側では内部圧力 p=140 mmHg 18.6 kPa ，外側では p=0と変化している．周方向応力 は血管組織との関係が強い．これは非線形性が強いため，次のFEM解析に関して大きい差を見ることができる．

3層モデル

次のFEM解析は3層モデルを表す．これは適切に簡素化した形状を使い，計算コストを減少させながらも正確な結果が得られるようにモデルを構築している．この比較的簡約化された解析は動脈の力学的環境を調べ，さらなる調査を進め，最終的には患者特有の形状を使って動脈疾患のリスクや進行を評価する出発点となり得る[Holzapfel, 2018]．

有限要素パッケージをロードし，$HistoryLengthを0に設定する：
モデル変数を設定する：

各層は，生体組織等の軟組織の力学応答に存在する非線形性を再生することができる，ヨーの超弾性モデルで記述される均質材料として表される．各層はそれぞれ，文献からの異なる材料パラメータを使う[Gervaso, 2018]．

ヨーの超弾性構造モデルは，ゴムや生体組織等の等方性で，ほぼ非圧縮性の非線形弾性材料の変形を表すための現象論的モデルである．このモデルはひずみ不変量のベキ級数として表されるひずみエネルギー密度関数によって記述される．

実装されたヨーモデルはほぼ非圧縮性モデルであり，パラメータ pars で"MaterialMoldel"->"Yeoh"を指定して使うことができる．ヨーモデルは"c1"，"c2"，"c3"の３つのパラメータを取り，ひずみエネルギー密度関数 の体積部分を表す多項式関数で記述される：

また，は第一不変量である．

それぞれの層について材料パラメータを指定する必要がある．このプロセスを簡単にするために，特性を各領域にマップするヘルパー関数を作成する．

材料特性を各層にマップする関数を定義する：
異なる層を指定するマーカーを定義する：
モデルの材料パラメータを指定する：

計算的に効率の良いPDE係数を設定する方法は，こちらをご覧いただきたい．

メッシュ

計算コストを下げ，層と層の関係をよりよく調べるために，解析を血管の2D断面に還元する．さらに，これは回転対称であるため，対応する対称条件を先端に課して，円環の四分の一だけを調べる．

形状の境界を生成する：
異なる領域を指定するマーカーを定義する：
メッシュを生成する：
メッシュを可視化する：

層の厚さ全体に十分な数の要素があることを保証することが大切である．導関数のサンプリングが不足せずにひずみと応力を推定するために，各領域に対して少なくとも3つの層が必要である．要素が多いほど正確さは向上するが，要素数は計算コストとバランスが取れている必要がある．

それぞれの部分領域についての要素数を調べる：

モデルの設定

このモデルでは2種類の境界条件が必要となる．縮小されたシミュレーション領域をモデル化するための対称条件と，内部圧力をモデル化するための境界条件である．対称条件についての詳細は，固体力学モノグラフの変位制約のセクションをご参照いただきたい．圧力をモデル化することについては，そのゴールは，課せられる外部荷重であるパラメトリック変数 でパラメトリックソルバを設定することである．パラメータ を使うことで，適用される力を徐々に増大していくことができる．これが必要なのは，この問題の非線形性が強すぎるため，1ステップでは解けないからである．非線形性の強い超弾性材料モデルを解くことに関する情報は，超弾性モノグラフのネオ・フックモデルのセクションに記載されている．

冠動脈の生理学的状態に関する内部圧力を指定する必要がある．140 mmHgの圧力を指定する[Ramanathan, 2005]．

合計の圧力と最大値 を指定する：

次のステップはこの圧力が適用される位置を指定することである．

圧力境界条件の配置：

対称条件は と でそれぞれ適用される．変位はそれぞれの方向で0に設定されるので，もう一方の方向には自由に動く．境界における不正確性を避けるために許容誤差を設定する．

小さい境界許容誤差を設定する：
境界条件を持つPDEモデルを作成する：
パラメトリックソルバを生成する：

解法

力学的問題に強い非線形性があるため，外部荷重をゆっくりと増やしていく．

荷重 をゆっくりと増やすことによって変位について解く：

シミュレーションが長く実行すると，再びシミュレーションを実行しなくても利用できるように，得られた解は保存される．

以下のコメントを無効にして解をExportする：
以下のコメントを無効にして，前にエキスポートされた解をImportする：

後処理

解を得た後，データを処理してサンプルの領域における応力場とひずみ場を調べることができる．

応力とひずみを抽出して，それらをもとのメッシュおよび変形したメッシュの両方にプロットすることができる．

もとの領域の境界を持つ変形したメッシュを作成し可視化する：

応力場とひずみ場を解析するために，これらの場を解から抽出する必要がある．

解からひずみ場を計算する：
ひずみ場をプロットする：

方向と 方向に沿ったひずみ場は，45°における二等分線 で対称になっている．このことは，この円環が暗示的に示している回転対称性に関係する．この理由で，別のオプションは等価ひずみをプロットすることである．等価ひずみはすべてのひずみ成分を1つのスカラーにまとめ理解しやすくする．等価ひずみは均質領域と不均質領域の以下の比較で使われる．

解から等価ひずみ場を計算する：
等価ひずみ場をプロットする：

厚さ全体での応力分布を調べるのもおもしろい． 方向と 方向に対する対称プロットも示す．外側の層よりも内側の層で高い値を示すコーシー応力をプロットすることができる．

解からコーシー応力場を計算する：

応力場を表すためには，kPa単位を使うと便利なので，変換係数を掛ける．

コーシー応力をkPa単位でプロットする：

別の方法として，応力のさまざまな成分全体の平均を表すミーゼス応力をプロットするというものもある．まさにそのプロットは内側の層で応力がいかに高く，内側から外側に行くにつれて減少することを示している．

ミーゼス応力を計算する：
変形したメッシュ上でミーゼス応力を計算する：
ミーゼス応力の最小値と最大値を計算する：
変形した血管全体のミーゼス応力をプロットする：

また主応力方向を調べると，応力分布の理解に役に立つ．

主応力の最初と2番目の方向を計算する：
主応力の最初と2番目の方向を可視化する：

主応力は応力行列の固有値であり，せん断応力 が0である斜め応力値 に対応する．主応力についての詳細は，固体力学モノグラフの主方向に記載されている．主応力方向をプロットすることによって，最初の方向が円環に対して接線方向であり同心円周全体への伸張にどのように関連しているかを示すことができる．さらに，2つ目の主方向は半径方向に現れ，内部境界圧力に主に関連しており外側に向かうにつれゼロになっていく．以下のプロットは主方向を示す一方，応力値は以下の比較で示され，前の解析解によって説明される．

均質化モデル

前述のモデルは，血管壁を記述するために3層のモデルを使った．ここで問題となるのは，1層モデルと3層モデルはどこが違うかということである．次では異なる2つの均質化アプローチを提示する．最初のアプローチは，シミュレーション領域のすべてに対して内側の層，内膜の材料特性を使っている．このアプローチは，特に内膜肥厚がある状態では妥当である．その場合内膜は主たる荷重負担層となるためである．2つ目のアプローチは3つすべての層の平均として計算される均質化を使う．

内膜値を使った均質化

材料特性を均質化するために，内側の層，内膜の特性だけを使うことができる．このことは，主たる荷重負担層となる内膜における大きい応力を示す，前の結果に従うと妥当と言える．

内膜の係数を使うパラメータ：
内膜均質化のためのパラメトリックソルバを生成する：
内膜値均質化で解く：
内膜均質化領域に対する後処理数量を計算する：

平均値を使った均質化

2つ目の均質化アプローチは3層に渡る特性を1つに均質化するために，数学平均を使うというものである．

材料係数の平均値を求めるヘルパー関数：

上のヘルパー関数でRoundを使うのは，第2ピオラ応力をコーシー応力に変換するときの数値的不安定性を避けるためである．平均値は近似に過ぎないので，数は最近の整数に変換される．整数自体は厳密数なので，コーシー応力への変換で問題に遭遇することは少ない．

係数の平均値を計算してそれを新しい材料パラメータに割り当てる：
平均均質化のためのパラメトリックソルバを作成する：
平均値の均質化で解く：
平均値で均質化された領域に対する後処理の数量を計算する：
最小・最大ひずみ測定を計算する：

比較

このセクションでは，3つの異なるFEMモデルを比較する．

層構造は，層の接触面全体でより正確な結果を示す．しかしこれには，均質化の計算時間の20倍以上の時間を必要とする．シミュレーションの関心部分によっては，一方のモデルより他方を使った方がよい場合がある．最適なモデルを選ぶには，応力とひずみの結果の差分を評価する必要がある．

各モデルにおける最小・最大の等価ひずみを計算する：
それぞれのモデルに対する等価ひずみをプロットする：
最小・最大のミーゼス応力を計算する：
ミーゼス応力の大域的最大値・最小値を計算してカラーマップをスケールする：
それぞれのモデルに対するミーゼス応力をプロットする：

均質，および均質化されたものは層のアプローチより密集している．層モデルと均質化モデルの違いは，層モデルには内部材料境界があるのに対し，均質化モデルにはそのような層が存在したいために生じる．

次に厚み全体の応力測定を抽出し，それを単軸の表にプロットすることができる．対称であるため，これらのすうりょはどの半径方向についても同じである．

以下のグラフの異なる領域をハイライトする補助関数を生成する：
におけるデータを抽出するヘルパー関数を作成する：
結果を抽出する点の数の定義：
厚み全体でのミーゼス応力をプロットする：
3層モデルの最初と2番目の主応力を抽出する：
内膜の係数を持つ均質モデルの最初と2番目の主応力を抽出する：
平均された係数を持つ均質モデルの最初と2番目の主応力を抽出する：
動径座標上に厚み全体の主 応力をプロットする：

最初の主応力は接線方向への伸長に関係している．実際，モデルが異なるとひずみ分布も異なる．図を見ると分かるように，これらの差はそれぞれのモデルでのミーゼス応力の差を引き起こしているものである．

動径座標上に厚み全体の主 応力をプロットする：

2番目の主応力はそれぞれのモデルで類似している．まさにこれは動径方向の応力であり，内側では境界値問題（平衡の法則）の内部圧力と対比するものである．外部圧力がないため，厚み全体でゼロになるはずである．

内膜と中膜の間の接触面全体で，応力分布は3つのモデルすべてて一致しているのも興味深い．さらに，3層モデルは内膜が主たる荷重負担層で他の層は荷重負担が低い．均質モデルと均質化モデルでは，応力場は継続的に内側から外側へと小さくなっている．

用語集

参考文献

1.  A. Gholipour, M. H. Ghayesh, A. Zander, R. Mahajan (2018). Three-dimensional biomechanics of coronary arteries. International Journal of Engineering Science, vol 130, pp. 93-114, https://doi.org/10.1016/j.ijengsci.2018.03.002

2.  S. Mishani, H. Belhoul-Fakir, C. Lagat, S. Jansen, B. Evans, M. Lawrence-Brown, (2021). Stress distribution in the walls of major arteries: implications for atherogenesis. Quantitative Imaging in Medicine and Surgery, vol 11, https://qims.amegroups.com/article/view/68560

3.  G. A. Holzapfel, R. W. Ogden (2018). Biomechanical relevance of the microstructure in artery walls with a focus on passive and active components. American Journal of Physiology-Heart and Circulatory Physiology, vol 315, https://doi.org/10.1152/ajpheart.00117.2018

4.  F. Gervaso, C. Capelli, L. Petrini, S. Lattanzio, L. Di Virgilio, F. Migliavacca (2016). On the effects of different strategies in modelling balloon-expandable stenting by means of finite element method. Journal of Biomechanics, vol 41, pp. 1206-1212, https://doi.org/10.1016/j.jbiomech.2008.01.027

5.  T. Ramanathan, H. Skinner, (2005). Coronary blood flow. Continuing Education in Anaesthesia Critical Care & Pain, vol 5, pp 1743-1816, https://doi.org/10.1093/bjaceaccp/mki012