CircleHeartシミュレーションを始めるログイン

Standard 74のしくみ

循環の接続、心筋の構成則、圧・流量の計算と解析方法。

baseline・プリセットの設定と検証

横長の数式・表は左右にスクロールできます。

保存
文書一式(HTML)

しくみ・設定・全検証記録をまとめて、オフラインで読めます。

設定・検証記録(JSON)数値表(CSV)

全体像と回路

心筋が張力を生み、圧差が血液を動かし、変わった体積が再び筋長と圧に影響します。心筋・心臓の形状・弁・血管の相互作用から、圧波形とPV loopを計算する0Dモデルです。

四心腔と11の主循環血管区画に、16の冠血管区画を接続します。各区画は血液を蓄える場所で、接続部が流れにくさを表します。名前は代表する領域を示すもので、カテーテル先端の正確な位置を指定するものではありません。中隔は心筋壁であり、血液を蓄えるnodeではありません。

体循環

  1. 左室 (LV)
  2. 近位大動脈 (Ao)
  3. 体動脈 (SA)
  4. 体抵抗血管側 (Art)
  5. 体毛細管 (Cap)
  6. 体静脈 (SV)
  7. 大静脈 (VC)
  8. 右房 (RA)

肺循環

  1. 右室 (RV)
  2. 近位肺動脈 (PA)
  3. 肺抵抗血管側 (PArt)
  4. 肺毛細管 (PCap)
  5. 肺細静脈側 (PVen)
  6. 肺静脈・左房流入口側 (PVein)
  7. 左房 (LA)
左房 —MV→ 左室、右房 —TV→ 右室で閉じた回路になります。矢印は流量の正方向です。逆流を許す接続では、負の流量も体積収支に含めます。冠循環はAoから分岐し、RAへ戻ります。

Ao・SA・Artを分けることで、近位の血液貯留と末梢へ向かう圧低下を別々に扱います。ただし、区画の間を脈波が伝わる時間や反射は表していません。SVは主な静脈貯留部、VCは胸腔圧の影響を受ける右房直前の区画です。肺側では、PCapに肺胞圧、PVen・PVeinに胸腔圧を作用させ、左房までの貯留と流入抵抗を分けています。この分割は、各血管を解剖学的に一本ずつ再現するためではありません。

主循環の全区画:役割と圧の基準
区画圧を決める関係外圧表示との対応
左室 (LV)心筋・形状の釣り合いPth+PperiP_{\mathrm{th}}+P_{\mathrm{peri}}LVP
左房 (LA)心筋・形状の釣り合いPth+PperiP_{\mathrm{th}}+P_{\mathrm{peri}}LAP
右室 (RV)心筋・形状の釣り合いPth+PperiP_{\mathrm{th}}+P_{\mathrm{peri}}RVP
右房 (RA)心筋・形状の釣り合いPth+PperiP_{\mathrm{th}}+P_{\mathrm{peri}}RAP / CVP
近位大動脈 (Ao)動脈の指数則0(基準圧)AoP
体動脈 (SA)動脈の指数則0(基準圧)ABP
体抵抗血管側 (Art)動脈の指数則0(基準圧)
体毛細管 (Cap)線形容量則0(基準圧)
体静脈 (SV)静脈型の非線形容量則0(基準圧)
大静脈 (VC)静脈型の非線形容量則PthP_{\mathrm{th}}
近位肺動脈 (PA)動脈の指数則PthP_{\mathrm{th}}PAP
肺抵抗血管側 (PArt)動脈の指数則PthP_{\mathrm{th}}
肺毛細管 (PCap)静脈型の非線形容量則PalvP_{\mathrm{alv}}
肺細静脈側 (PVen)静脈型の非線形容量則PthP_{\mathrm{th}}
肺静脈・左房流入口側 (PVein)静脈型の非線形容量則PthP_{\mathrm{th}}

AoP・PAPはAo・PA区画の内圧で、表示の段階でZcQを足しません。ABPはSAの代表圧で、上腕カフ圧の再現ではありません。CVPは平均右房圧、PCWPは平均左房圧を代用する表示であり、肺毛細管楔入の手技は計算していません。PV loopは心室経壁圧を使います。

V˙i=jNijQj,Nij={+1j enters i1j leaves i0otherwise,i=131Vi=TBV\dot V_i=\sum_j N_{ij}Q_j,\qquad N_{ij}=\begin{cases}+1&j\text{ enters }i\\-1&j\text{ leaves }i\\0&\text{otherwise}\end{cases},\qquad\sum_{i=1}^{31}V_i=TBV

接続行列Nは、プリセットの接続表から組み立てられます。どの接続も一方から出た量が他方へ入るため、閉回路全体では血液量が保存されます。心筋体積・心嚢液量はこの合計とは別です。TBVを変更しない限り、圧や拍出量を合わせるための血液の追加・削除は行いません。

構成と数式

以下はこのモデルに共通する構成則です。係数の採用値・形状・初期状態は、baselineまたは各プリセットの「設定」にまとめています。両方を合わせて、選んだ作動点のモデルを組み立てられます。

係数・形状・初期状態を見る

興奮とCa:収縮のきっかけ

心房・心室が興奮すると、細胞内のCa濃度が一過性に上昇します。このCaの変化を入力として、筋原線維の張力を計算します。

x˙r=xr/τr,x˙d=xd/τd,[Ca]=Ca0+g(xdxr)\dot x_r=-x_r/\tau_r,\quad \dot x_d=-x_d/\tau_d,\qquad [Ca]=\mathrm{Ca}_0+g(x_d-x_r)

時定数が異なる場合の式です。興奮時に二つの状態変数へ同じ量を加えます。速く減る成分と遅く減る成分の差によって、Caが立ち上がり、その後低下します。

数式・仮定を詳しく

baselineは規則的な洞調律です。心房の興奮と、遅れて起こる心室の興奮を扱います。心室では二つの指数関数を組み合わせてCa波形を定め、心房には別の設定を使います。細胞内のCa輸送や活動電位の全過程は解きません。前の拍の成分が残るため、Ca₀と実際の1拍の最低Ca濃度は異なります。

xr,  xdx_r,\;x_d
Ca波形を作る無次元の状態変数
τr,  τd\tau_r,\;\tau_d
各成分が減衰する時定数(s)
Ca0,  g\mathrm{Ca}_0,\;g
興奮が長く途絶えたときのCa濃度と、波形の振幅係数(µM)

興奮時刻・状態の更新

時間の単位はsです。規則的洞調律では、心房興奮の間隔T=60/HR、心室興奮はその120 ms後(房室伝導80 ms+遠位伝導40 ms)です。各心房・心室のCaへの入力は、対応する興奮の12 ms後に加えます。各壁は独立した二つのCa状態を持ちます。

qは興奮ごとの入力強度です。baselineでは両心房・三つの心室壁ともq=1です。イベント間は指数関数で厳密に減衰し、イベント直後に両状態へqを加えます。下式の周期解は規則的なq=1の場合で、心腔容積やLand状態の初期値を定めるものではありません。

tA,k=tA,0+kT,tV,k=tA,k+0.120,tCa,w,k=tA or V,k+0.012t_{A,k}=t_{A,0}+kT,\quad t_{V,k}=t_{A,k}+0.120,\quad t_{Ca,w,k}=t_{A\text{ or }V,k}+0.012
xj(t+Δt)=xj(t)eΔt/τj,xj(tk+)=xj(tk)+qk(j=r,d)x_j(t+\Delta t)=x_j(t)e^{-\Delta t/\tau_j},\quad x_j(t_k^+)=x_j(t_k^-)+q_k\quad(j=r,d)
xj(tk+)=11eT/τj(periodic, qk=1)x_j(t_k^+)=\frac{1}{1-e^{-T/\tau_j}}\quad\text{(periodic, }q_k=1\text{)}

心拍間隔とCa入力強度

心室には拍間隔依存の離散的な負荷状態Lを持たせます。Iは直前の心室興奮からの間隔、aは回復割合、qは次のCa入力強度です。LはSR負荷を規格化した量で、実測のSR内Ca濃度ではありません。係数と規則的洞調律での基準状態は下の表に示します。基準心拍数に対応するL・aとq=1が固定点を作ります。

ここで記載する起動条件は規則的洞調律・補助循環なしです。異所性刺激・ペーシング・補助循環装置を有効にした場合の全モードは、このbaselineの方程式系には含めません。

ak=1eIk/τrec,qk=akβLk,Lk+1=Lk(1r)qk+γ(1hak)a_k=1-e^{-I_k/\tau_{\mathrm{rec}}},\quad q_k=a_k\beta L_k,\qquad L_{k+1}=L_k-(1-r)q_k+\gamma(1-ha_k)

筋原線維:Caから張力へ

アクチンとミオシンが結びつくクロスブリッジの割合と、その歪みから能動張力を計算します。同じCa濃度でも、心筋の長さや短縮速度、過去に受けた負荷によって力が変わります。

Ta=h(λ)Trefrs[S(1+ζs)+Wζw]T_a=\frac{h(\lambda)T_{\mathrm{ref}}}{r_s}\left[S(1+\zeta_s)+W\zeta_w\right]

Land由来の能動線維応力の出力式(Pa)。W/Sは弱・強結合の割合、ζは歪み、hは長さ依存係数。心腔圧そのものではありません。

数式・仮定を詳しく

Landら(2017)の収縮モデルを基礎に、心室の張力係数・Ca感受性・結合速度を調整しています。さらに、Ca結合が低下したときに強結合の余剰分を非結合状態へ戻す項を加えました。各状態の割合の合計は保ちますが、追加した機構の生理的妥当性には別途検証が必要です。

Ta,  TrefT_a,\;T_{\mathrm{ref}}
能動線維応力と、その大きさを決める係数(Pa)
λ,  h(λ)\lambda,\;h(\lambda)
線維の長さ/基準長と、長さ依存の補正係数(無次元)
W,  S,  rsW,\;S,\;r_s
弱結合・強結合の割合と、基準の強結合割合(無次元)
ζw,  ζs\zeta_w,\;\zeta_s
結合したクロスブリッジの歪み(無次元)。短縮速度と過去の変形に依存
Jexit=kmax(θnθn+cn)pmax(SrW,0),S˙exit=Jexit,U˙exit=JexitJ_{\mathrm{exit}}=k_{\mathrm{max}}\left(\frac{\theta^n}{\theta^n+c^n}\right)^p\max(S-rW,0),\quad \dot S|_{\mathrm{exit}}=-J_{\mathrm{exit}},\quad \dot U|_{\mathrm{exit}}=J_{\mathrm{exit}}

追加した強結合離脱の式。cはCaが結合したトロポニンの割合、θはその基準値、nとpは感受性を決める指数です。kmaxは最大離脱速度(s⁻¹)、r=kws/ksuは歪みゼロでの強/弱結合比、Uは非結合割合です。SがrWを超えた分だけを移します。

六つの状態と保存条件

各壁でc=Ca結合トロポニンの割合、b=結合が阻害された割合、W=弱結合、S=強結合、ζw・ζs=架橋の歪みを持ちます。結合可能な非結合割合Uは独立状態ではなく、U=1−b−W−Sです。0<c≤1、b・U・W・S≥0を満たす解を使います。bは血管・弁の二次抵抗係数とは別の記号です。

Caは自由Ca濃度(µM)、λはLandに入力する伸長比です。gwu・gsuは歪みに伴う離脱速度(s⁻¹)です。Jexitは上で示した追加離脱項で、心室壁に適用します。心房ではJexit=0です。離脱した割合はUへ戻ります。

c˙=kTRPN{(Ca/Ca50)nTRPN(1c)c},b˙=kbmin(cnTm/2,100)UkucnTm/2b,W˙=kuwU(kwu+kws+gwu)W,S˙=kwsW(ksu+gsu)SJexit,ζ˙w=Awλ˙cwζw,ζ˙s=Asλ˙csζs.\begin{aligned}\dot c&=k_{\mathrm{TRPN}}\{(\mathrm{Ca}/\mathrm{Ca}_{50})^{n_{\mathrm{TRPN}}}(1-c)-c\},\\\dot b&=k_b\min(c^{-n_{\mathrm{Tm}}/2},100)U-k_u c^{n_{\mathrm{Tm}}/2}b,\\\dot W&=k_{\mathrm{uw}}U-(k_{\mathrm{wu}}+k_{\mathrm{ws}}+g_{\mathrm{wu}})W,\\\dot S&=k_{\mathrm{ws}}W-(k_{\mathrm{su}}+g_{\mathrm{su}})S-J_{\mathrm{exit}},\\\dot\zeta_w&=A_w\dot\lambda-c_w\zeta_w,\qquad\dot\zeta_s=A_s\dot\lambda-c_s\zeta_s.\end{aligned}
λc=min(λ,1.2),Ca50=Ca50,ref+β1(λc1),h(λ)=max{0,1+β0[λc+min(λc,0.87)1.87]},gwu=γwζw,gsu=γsmax(ζs1,ζs,0).\begin{aligned}\lambda_c&=\min(\lambda,1.2),\\\mathrm{Ca}_{50}&=\mathrm{Ca}_{50,\mathrm{ref}}+\beta_1(\lambda_c-1),\\h(\lambda)&=\max\{0,1+\beta_0[\lambda_c+\min(\lambda_c,0.87)-1.87]\},\\g_{\mathrm{wu}}&=\gamma_w|\zeta_w|,\quad g_{\mathrm{su}}=\gamma_s\max(-\zeta_s-1,\zeta_s,0).\end{aligned}

派生係数と心臓壁への結合

表に示す独立係数から、以下の係数を計算します。rₛ・r𝓌は基準結合割合、θ=TRPN50です。k・cの単位はs⁻¹、Aは無次元です。Ca50が正でない場合は有効な材料条件ではありません。長さ依存をλ=1.2で飽和させる処理と、トロポニン項の上限100も採用した式の一部です。

形状の対数歪みをeとするとλ=s₀ exp(e)。s₀は筋長基準の倍率で、baselineは全壁で1です。短縮速度には(λn+1−λn)/Δtを使います。Landの出力Tₐをλ・線維配向割合fo・活動割合faで変換し、Kirchhoff応力として受動・粘弾性応力に加えます。fo=fa=1です。心室の収縮力倍率はTrefだけに、受動倍率は受動応力と粘弾性弾性率に掛かります。

kb=kuθnTm(1rs)(1rw),kwu=kuw(1/rw1)kws,ksu=kwsrw(1/rs1),Aw=As=Aeffrs(1rs)rw+rs,cw=ϕkuw(1rw)rw,cs=ϕkws(1rs)rwrs,τf=λfofaTa+τpass+τvis.\begin{aligned}k_b&=\frac{k_u\theta^{n_{\mathrm{Tm}}}}{(1-r_s)(1-r_w)},\quad k_{\mathrm{wu}}=k_{\mathrm{uw}}(1/r_w-1)-k_{\mathrm{ws}},\\k_{\mathrm{su}}&=k_{\mathrm{ws}}r_w(1/r_s-1),\quad A_w=A_s=\frac{A_{\mathrm{eff}}r_s}{(1-r_s)r_w+r_s},\\c_w&=\frac{\phi k_{\mathrm{uw}}(1-r_w)}{r_w},\quad c_s=\frac{\phi k_{\mathrm{ws}}(1-r_s)r_w}{r_s},\\\tau_f&=\lambda f_o f_aT_a+\tau_{\mathrm{pass}}+\tau_{\mathrm{vis}}.\end{aligned}
Land et al. 2017

受動特性:伸びにくさと変形の履歴

拡張期の圧には、心筋が伸ばされたときの抵抗と、時間をかけて緩む粘弾性が関わります。これらを能動張力に加えて、心臓壁に働く応力を求めます。

e=lnλg,τpass=dΨde,τvis=Ev(eα),α˙=eατve=\ln\lambda_g,\quad \tau_{\mathrm{pass}}=\frac{d\Psi}{de},\quad \tau_{\mathrm{vis}}=E_v(e-\alpha),\quad \dot\alpha=\frac{e-\alpha}{\tau_v}

受動・粘弾性部分の関係式です。能動部分はLandの応力に線維伸長比と配向・活動割合を掛け、同じ応力基準に変換して加えます。

数式・仮定を詳しく

心室の受動応力は、伸びに対して増加する弾性エネルギーから求めます。心房には別の受動則を使います。粘弾性は、変形の履歴を一つの変数で持つMaxwell要素で表します。心室の受動則はKlotzの心室全体のEDPVRを参考にした設計であり、個々の心筋の物性を直接測定した値ではありません。

λg,  e,  α\lambda_g,\;e,\;\alpha
形状から求めた伸長比、その対数歪み、粘性部分の歪み(無次元)
Ψ,  τpass,  τvis\Psi,\;\tau_{\mathrm{pass}},\;\tau_{\mathrm{vis}}
弾性エネルギー密度(J/m³)と、受動・粘性のKirchhoff応力(Pa)
Ev,  τvE_v,\;\tau_v
粘弾性の弾性率(Pa)と材料の緩和時間(s)。LVPから測る弛緩時定数τとは別

心室の弾性エネルギー

e=ln λg、p=Hδ(e)、q=Hδ(−e)とします。K0は原点付近の弾性率、Kt・aは伸張時の指数的な硬化、Kcは圧縮側の追加弾性率です。Hδは正の歪みだけを滑らかに取り出す関数で、u=z/δです。心室三壁には同じ式を使い、エネルギーとその微分に各壁の受動倍率sₚ,wを掛けます。倍率はbaselineの設定表に示します。

Hδ(z)={0z0δ(u312u4)0<z<δzδ/2zδH_\delta(z)=\begin{cases}0&z\le0\\\delta(u^3-\tfrac12u^4)&0<z<\delta\\z-\delta/2&z\ge\delta\end{cases}
Ψ(e)=12K0e2+Kta2(eap1ap)+12Kcq2,τpass=sp,wdΨde\Psi(e)=\tfrac12K_0e^2+\frac{K_t}{a^2}(e^{ap}-1-ap)+\tfrac12K_cq^2,\qquad \tau_{\mathrm{pass}}=s_{p,w}\,\frac{d\Psi}{de}

心房の受動則と粘弾性

心房は等二軸変形に縮約した別の材料則を使います。下式のsₚ,wは各壁の受動倍率です。λg=exp(e)、C1・C2・C3はPa、C4は無次元です。線維項は伸張側のみ働きます。対応するエネルギーはΨ(e)=∫₀ᵉτpass(s)dsで定まります。支持する歪み範囲は−0.5≤e≤0.5です。左房由来の材料を右房にも適用する点は外挿です。

全壁の粘弾性は1状態のMaxwell枝です。αは粘性歪み、Evは弾性率、τvは緩和時間です。下式は採用している後退Euler更新で、各力学ステップにつき1回適用します。Evの表は受動倍率を適用済みです。

τpass=sp,w[4C1(λg2λg4)+4C2(λg4λg2)+{C3(eC4(λg1)1)e>00e0]\tau_{\mathrm{pass}}=s_{p,w}\left[4C_1(\lambda_g^2-\lambda_g^{-4})+4C_2(\lambda_g^4-\lambda_g^{-2})+\begin{cases}C_3(e^{C_4(\lambda_g-1)}-1)&e>0\\0&e\le0\end{cases}\right]
αn+1=αn+(Δt/τv)en+11+Δt/τv,τvis,n+1=Ev(en+1αn+1)\alpha_{n+1}=\frac{\alpha_n+(\Delta t/\tau_v)e_{n+1}}{1+\Delta t/\tau_v},\qquad \tau_{vis,n+1}=E_v(e_{n+1}-\alpha_{n+1})

Klotz et al. 2006心室全体のEDPVRの参照。採用した心筋の受動則そのものの出典ではありません。

Moyer et al. 2015心房受動材料の出発点。ここでは等二軸変形に縮約して使います。

心房・心室:張力から圧へ

左室自由壁・中隔・右室自由壁を結合し、左右の心室が中隔を介して影響し合います。これに左右心房を加えた五壁構成です。圧は決め打ちの時間波形ではなく、壁応力と形状の釣り合いから求まります。

δW=wVm,wτf,wδlnλw,Pcavity=Ptm+Pext\delta W=\sum_w V_{\mathrm{m},w}\tau_{\mathrm{f},w}\,\delta\ln\lambda_w,\qquad P_{\mathrm{cavity}}=P_{\mathrm{tm}}+P_{\mathrm{ext}}

壁をわずかに変形させたときの仕事と、内外の圧の関係です。壁の仕事と心腔の圧–体積仕事を対応させて圧を求め、中隔位置と接合円半径についても力の釣り合いを解きます。

数式・仮定を詳しく

左室自由壁・中隔・右室自由壁を、共通の円で接する三つの球冠として表します。これがTriSegの形状近似です。患者の3D画像を再現するものではなく、局所的な壁応力やねじれは扱いません。PV loopの縦軸は経壁圧、圧波形の縦軸は心腔内圧です。心膜圧や胸腔圧が変わると、両者の差も変わります。

δW,  Vm,w\delta W,\;V_{\mathrm{m},w}
仮想的な微小変形に伴う仕事(J)と、壁wの心筋体積(m³)
τf,w,  λw\tau_{\mathrm{f},w},\;\lambda_w
各壁のKirchhoff線維応力(Pa)と、形状から決まる伸長比。ここでの応力を表すτは、LVPから測る弛緩時定数とは別の量
Pcavity,  Ptm,  PextP_{\mathrm{cavity}},\;P_{\mathrm{tm}},\;P_{\mathrm{ext}}
心腔内圧、壁を内側から広げる経壁圧、壁の外側にかかる圧。同じ単位で加算
LVRVSEPy接合円の断面右室方向 +
三つの壁が同じ円で接する球冠近似の模式図です。断面形状や壁厚の実測図ではありません。hは接合面から各球冠頂点までの符号付き高さです。

三つの球冠と心筋長

以下の力学式はm・m²・m³・PaのSI単位で書きます。VL・VRは心腔血液量、Mwは各壁の心筋体積、vSは中隔の符号付き球冠体積、y>0は三壁が接する円の半径です。心筋体積MはTBVには含めません。球冠高さhは右室方向を正とし、左室自由壁では通常負になります。

各壁wについて球冠体積vwからhを求め、面積Aw・曲率κw・厚さ補正zwを計算します。Aref,wは基準中壁面積です。ここでのzwはLandの架橋歪みζとは別物です。形状から得たewを材料則へ入力します。

vL=VL12(ML+MS)+vS,vR=VR+12(MR+MS)+vSv_L=-V_L-\tfrac12(M_L+M_S)+v_S,\qquad v_R=V_R+\tfrac12(M_R+M_S)+v_S
vw=πhw(hw2+3y2)6,Aw=π(hw2+y2),κw=2hwhw2+y2v_w=\frac{\pi h_w(h_w^2+3y^2)}6,\quad A_w=\pi(h_w^2+y^2),\quad\kappa_w=\frac{2h_w}{h_w^2+y^2}
zw=3κwMw2Aw,ew=12lnAwAref,wzw2120.019zw4z_w=\frac{3\kappa_wM_w}{2A_w},\qquad e_w=\frac12\ln\frac{A_w}{A_{\mathrm{ref},w}}-\frac{z_w^2}{12}-0.019z_w^4

内部の釣り合いと圧の決定

Fは球冠体積方向、Gは接合円半径方向の一般化力です。各壁の応力は材料状態と候補形状から求め、その時点の応力を掛けて歪みの幾何学微分を計算します。これを『能動応力がポテンシャルから生じる』と仮定して微分するわけではありません。

心腔容積を与え、三壁の力が釣り合うvS・yを連立して求めます。その解からLV・RV経壁圧が定まります。両心房には球状の一線維近似を使い、Vrefは表の基準心腔血液量です。圧に任意の時間波形や追加の壁別圧倍率を掛けません。

Fw=Mwτf,wewvwy,Gw=Mwτf,wewyvwF_w=M_w\tau_{\mathrm{f},w}\left.\frac{\partial e_w}{\partial v_w}\right|_y,\qquad G_w=M_w\tau_{\mathrm{f},w}\left.\frac{\partial e_w}{\partial y}\right|_{v_w}
FL+FS+FR=0,GL+GS+GR=0,Ptm,L=FL,Ptm,R=FRF_L+F_S+F_R=0,\quad G_L+G_S+G_R=0,\qquad P_{\mathrm{tm},L}=-F_L,\quad P_{\mathrm{tm},R}=F_R
eA=13lnVA+MA/2Vref,A+MA/2,Ptm,A=MAτf,A3(VA+MA/2)e_A=\frac13\ln\frac{V_A+M_A/2}{V_{\mathrm{ref},A}+M_A/2},\qquad P_{\mathrm{tm},A}=\frac{M_A\tau_{\mathrm{f},A}}{3(V_A+M_A/2)}

胸腔圧・心膜圧の加算

心膜内を占める容積VHは、四心腔の血液量+五壁の心筋体積+処方した心嚢液量です。冠血液量をここへもう一度加える処理はありません。四心腔に同じ心膜外圧を加えます。V0は基準容量、P*は圧係数、kは硬化係数です。

Hは心膜が張り始める領域を滑らかにした関数です。δ=0.001、u=(x+δ)/(2δ)とします。H′はxでの微分です。baselineでは圧オフセット0、心嚢液0 mL、胸腔・肺胞の呼吸振幅0です。心膜の張りがなければ、心腔内圧と経壁圧は一致します。

VH=c=LA,LV,RA,RVVc+wMw+Vfluid,x=(VHV0)/V0V_H=\sum_{c=LA,LV,RA,RV}V_c+\sum_w M_w+V_{\mathrm{fluid}},\quad x=(V_H-V_0)/V_0
H(x)={0xδδ(2u3u4)x<δxxδH(x)=\begin{cases}0&x\le-\delta\\\delta(2u^3-u^4)&|x|<\delta\\x&x\ge\delta\end{cases}
Pperi=Poffset+P[ekH(x)1]H(x),Pc=Ptm,c+Pth+PperiP_{peri}=P_{\mathrm{offset}}+P_*[e^{kH(x)}-1]H'(x),\qquad P_c=P_{\mathrm{tm},c}+P_{\mathrm{th}}+P_{peri}

Lumens et al. 2009 · TriSeg三つの壁の形状近似の基礎です。心筋の構成則や圧への変換は本モデルの実装に従い、原著モデル全体の再現とは区別します。

四弁:圧差から流れへ

僧帽弁・大動脈弁・三尖弁・肺動脈弁は、共通の開口・圧損失の構造を使います。面積や逆流口の設定は弁ごとに異なります。

ΔP=RQ+B(A)QQ,B(A)ρ/A2\Delta P=R Q+B(A)Q|Q|,\qquad B(A)\propto \rho/A^2

弁を通る流れがある場合の式です。Qと面積Aの単位換算は係数に含めます。完全閉鎖時は、圧差が残っていてもQ=0となる閉鎖条件を別に適用します。

数式・仮定を詳しく

各時刻の流量は圧差から代数的に求め、流量そのものの慣性は持たせません。弁尖の開口割合は過去の状態を引き継ぎ、0〜1の範囲で変わります。有効開口面積(EOA)には縮流の影響を含めるため、流出係数Cdを別に掛けません。弁通過後の圧回復は組み込んでいません。この圧差を、Dopplerのジェット速度から求める勾配や、異なる時刻の圧ピーク同士の差と読み替えることはできません。

ΔP,  Q\Delta P,\;Q
弁前後の同時圧差(mmHg)と流量(mL/s)。順行を正とする
A,  ρA,\;\rho
その時点の有効開口面積と血液密度。面積は最大面積と開口割合から求める
R,  BR,\;B
線形抵抗(mmHg·s/mL)と二次損失係数(mmHg·s²/mL²)

開口状態の時間発展

ξは開口割合(0〜1)、ΔPは上流圧−下流圧です。dは不感帯、pₒは圧オフセット(baselineは四弁とも0 mmHg)、kₒは開口感度です。Fεは正の開口駆動を滑らかにした関数で、ε=0.1 mmHgです。ξ∞が直前のξより大きければ開口時定数、それ以外は閉鎖時定数を使います。

Fϵ(z)={0z0z2/(2ϵ)0<z<ϵzϵ/2zϵ,ξ=1ekoFϵ(ΔPdpo)F_\epsilon(z)=\begin{cases}0&z\le0\\z^2/(2\epsilon)&0<z<\epsilon\\z-\epsilon/2&z\ge\epsilon\end{cases},\quad \xi_\infty=1-e^{-k_oF_\epsilon(\Delta P-d-p_o)}
ξ˙=ξξτξ,ξn+1=ξn+(Δt/τξ)ξ(ΔPn+1)1+Δt/τξ\dot\xi=\frac{\xi_\infty-\xi}{\tau_\xi},\qquad \xi_{n+1}=\frac{\xi_n+(\Delta t/\tau_\xi)\xi_\infty(\Delta P_{n+1})}{1+\Delta t/\tau_\xi}

面積・方向・完全閉鎖

Amaxは最大順行EOA、Arは閉鎖時の逆流EOAです。順行では残存開口Arにξで開く部分を加え、逆行ではξによらずArだけを使います。従ってArを増やすと、逆行だけでなく閉鎖近くの順行流も変わります。baselineは全弁Ar=0です。

QはmL/s、Aはcm²、圧はmmHg、ρ=1060 kg/m³、cP=133.322387415 Pa/mmHgです。面積が0ならQ=0とし、残る圧差を閉鎖支持反力として扱います。小さな面積の下限や流量の平滑化は置きません。Rは面積に応じて再スケールしません。

A={Ar+ξ(AmaxAr)ΔP0ArΔP<0,B(A)=ρ2cP(106104A)2A=\begin{cases}A_r+\xi(A_{\mathrm{max}}-A_r)&\Delta P\ge0\\A_r&\Delta P<0\end{cases},\quad B(A)=\frac{\rho}{2c_P}\left(\frac{10^{-6}}{10^{-4}A}\right)^2
Q={0A=0 or ΔP=0sgn(ΔP)2ΔPR+R2+4B(A)ΔPotherwiseQ=\begin{cases}0&A=0\text{ or }\Delta P=0\\\operatorname{sgn}(\Delta P)\frac{2|\Delta P|}{R+\sqrt{R^2+4B(A)|\Delta P|}}&\text{otherwise}\end{cases}

血管:蓄える・流す

血管のコンプライアンスは血液を蓄える性質、抵抗は流れにくさを表します。体循環と肺循環を閉じた回路として結び、各部の体積と圧を同時に求めます。

dVidt=Qin,iQout,i,Ci(Ptm)=dVidPtm,i\frac{dV_i}{dt}=\sum Q_{\mathrm{in},i}-\sum Q_{\mathrm{out},i},\qquad C_i(P_{\mathrm{tm}})=\frac{dV_i}{dP_{\mathrm{tm},i}}

体積保存とコンプライアンスの定義。静脈を含む圧–体積関係は必ずしも一定Cの直線ではありません。

数式・仮定を詳しく

大動脈と肺動脈の近位部では、流体の慣性係数Lを0としています。AoPとPAPには、それぞれの血管区画の圧を表示します。特性インピーダンスと流量の積(ZcQ)を表示時に足す処理はありません。SAは下流の体動脈をまとめた区画です。圧波の伝播や反射を解かないため、カフや特定の動脈ライン位置の波形とは区別します。

Vi,  Qin,  QoutV_i,\;Q_{\mathrm{in}},\;Q_{\mathrm{out}}
区画iの血液量(mL)と、流入・流出量(mL/s)
Ci,  Ptm,iC_i,\;P_{\mathrm{tm},i}
圧に応じたコンプライアンス(mL/mmHg)と、血管の経壁圧(mmHg)

血管区画の圧–容量関係

圧はmmHg、容積はmL、流量はmL/sです。Pは内圧、p=P−Pextは経壁圧です。Vuはp=0での容量で、V−Vuをstressed volumeと呼びます。動脈のVsは表では硬さ・コンプライアンス倍率を適用した値です。Ao・SA・ArtはVs,元×0.65/1.42、PA・PArtはVs,元/1.42です。

静脈型区画は、虚脱時Cc・開通時Co・高圧伸展時Cdを滑らかに接続します。S(z)=ln(1+exp z)、σ(z)=1/(1+exp(−z))、ΔS(p;a,d)=S((p−a)/d)−S(−a/d)とします。静脈トーンuはVu=Vu,元−G uとして作用し、baselineのu=0.15を適用したVuを表に示します。

静脈の逆計算は−20〜45 mmHgに制限し、範囲外では端の圧を返します。動脈は(V−Vu)/Vsの下限をln(0.05)とします。これは正常範囲ではなく計算上の境界で、境界に達した結果を生理的な予測として外挿できません。コンプライアンス評価には下限10⁻⁴ mL/mmHgがありますが、V(p)自体を別の曲線に置き換える処理ではありません。

p=P0(emax[(VVu)/Vs,ln0.05]1)(arterial),p=(VVu)/C(linear)p=P_0\left(e^{\max[(V-V_u)/V_s,\ln0.05]}-1\right)\quad\text{(arterial)},\qquad p=(V-V_u)/C\quad\text{(linear)}
V(p)=Vu+Ccp+(CoCc)doΔS(p;po,do)(CoCd)dsΔS(p;ps,ds)V(p)=V_u+C_cp+(C_o-C_c)d_o\Delta S(p;p_o,d_o)-(C_o-C_d)d_s\Delta S(p;p_s,d_s)
dVdp=Cc+(CoCc)σ(ppodo)(CoCd)σ(ppsds)\frac{dV}{dp}=C_c+(C_o-C_c)\sigma\left(\frac{p-p_o}{d_o}\right)-(C_o-C_d)\sigma\left(\frac{p-p_s}{d_s}\right)

冠循環・外圧:心臓と循環の結合

冠動脈への流入と静脈への還流も全体の体積収支に含みます。心筋収縮による血管の圧迫が冠血流に影響し、心膜・胸腔圧は心腔や血管を外側から負荷します。

数式・仮定を詳しく

冠循環は、左前下行枝・回旋枝・右冠動脈の領域を、それぞれ心外膜側と心内膜側に分けます。血管の抵抗と容量に加え、心筋内圧による圧迫・虚脱と自己調節を扱います。一部の係数には、正常成人を想定した暫定値を使っています。baselineでは補助循環・弁逆流・呼吸性変動を加えていません。

AoArtC1C1C2C2CVRAR1RmR2心外膜側心内膜側
LAD・LCx・RCAそれぞれにこの分岐があり、CVだけを共有します。Art・CVには共通心外圧、C1・C2には各領域・層の心筋内圧が作用します。

冠血管の容量と抵抗

各領域(LAD・LCx・RCA)は1つの大動脈側容量Artを持ち、そこから心外膜側・心内膜側の2経路へ分かれます。各経路は近位容量C1と遠位容量C2を持ち、共通冠静脈容量CVへ合流して右房へ戻ります。C1・C2という名称は区画名で、コンプライアンスそのものの数値ではありません。合計16容量を主循環と重複なくTBVへ含めます。

全冠区画に以下の圧–容量関係を使います。v=V/Vref、m=4、n=2。CrefはV=Vrefでの接線コンプライアンスです。虚脱による抵抗倍率fには別の基準容量Vhを使うため、Vrefと混同しないでください。a=0.67、x=min(1,max(0,V/Vh))とします。

各枝の流量は両端圧差/有効抵抗です。R1には自己調節トーンθとC1の虚脱倍率、Rmには両区画の虚脱倍率の幾何平均、R2にはC2の虚脱倍率を掛けます。構造的な抵抗倍率はbaselineで1、局所狭窄による追加損失は0です。

p=P(vmvn),P=VrefCref(m+n),V>0p=P_*(v^m-v^{-n}),\quad P_*=\frac{V_{\mathrm{ref}}}{C_{\mathrm{ref}}(m+n)},\quad V>0
f(V)={a+(1a)x2(32x)}2f(V)=\{a+(1-a)x^2(3-2x)\}^{-2}
R1,eff=R1θf(VC1),Rm,eff=Rmf(VC1)f(VC2),R2,eff=R2f(VC2)R_{1,eff}=R_1\theta f(V_{C1}),\quad R_{\mathrm{m},eff}=R_m\sqrt{f(V_{C1})f(V_{C2})},\quad R_{2,eff}=R_2f(V_{C2})

心筋内圧と拍ごとの自己調節

ArtとCVの外圧は共通心外圧Pe=Pth+Pperi、C1とC2の外圧は心筋内圧PIMです。壁重みwL・wS・wR、層深さd、短縮圧係数Kは表を参照します。中隔の深さsはLAD/LCxでd、RCAで1−dです。eMVC,wは直前の僧帽弁閉鎖時に更新する各壁の歪み。Fは弁の節で定義した正部分関数で幅0.005です。

短縮の基準更新には、受理されたステップでMV順行流が1 mL/s超から1 mL/s以下へ変わったときの歪みを使います。更新した基準は次のステップから使います。これは画像上の弁尖接合時刻や、指標計算で補間した弁閉鎖時刻と同じ定義ではありません。

自己調節は各心周期を終えたときに更新し、その周期内ではトーンを保持します。Q̄mはRmを通る符号付き流量の周期平均、Qtargetは領域・層の安静目標です。ℓ=ln θ、Tは周期長、τa=25 s、θmin=4/45、θmax=2です。baselineでは需要倍率1・充血駆動0で、下式を使います。これは酸素需要から自律的に収縮力を変える機構ではありません。

PIM=Pe+(wLd+wSs)Ptm,L+[wRd+wS(1s)]Ptm,R+KwwwF0.005(1eeweMVC,w)P_{\mathrm{IM}}=P_e+(w_Ld+w_Ss)P_{\mathrm{tm},L}+[w_Rd+w_S(1-s)]P_{\mathrm{tm},R}+K\sum_w w_wF_{0.005}(1-e^{e_w-e_{\mathrm{MVC},w}})
Qˉm,k=1Tn in cycle kΔtnQm,n+1\bar Q_{\mathrm{m},k}=\frac{1}{T}\sum_{n\text{ in cycle }k}\Delta t_n Q_{m,n+1}
k+1=cliplnθmin,lnθmax[k+Tτalnmax(Qˉm,kQtarget,0.05)]\ell_{k+1}=\operatorname{clip}_{\ln\theta_{\mathrm{min}},\ln\theta_{\mathrm{max}}}\left[\ell_k+\frac{T}{\tau_a}\ln\max\left(\frac{\bar Q_{\mathrm{m},k}}{Q_{\mathrm{target}}},0.05\right)\right]

酸素輸送:血流から供給と消費を計算

1拍平均の血流と、Hb・吸入酸素濃度・酸素消費量の設定から、全身の酸素需給を計算します。心筋の収縮モデルとは別の計算です。

D˙O2=10COCaO2,V˙O2=10CO(CaO2CvO2)\dot D_{\mathrm{O_2}}=10\,\mathrm{CO}\,C_{\mathrm{aO}_2},\qquad \dot V_{\mathrm{O_2}}=10\,\mathrm{CO}(C_{\mathrm{aO}_2}-C_{\mathrm{vO}_2})

数式・仮定を詳しく

肺胞気式、酸素解離曲線、シャント血の混合を使って動脈血酸素含量を求めます。Fickの関係から必要な混合静脈血酸素含量を逆算し、負の含量を必要とする条件は計算上成立しないと判定します。局所組織の酸素拡散や代謝適応の全過程は再現しません。

CO\mathrm{\mathrm{CO}}
1拍平均の血流量(L/min)
CaO2,  CvO2C_{\mathrm{aO}_2},\;C_{\mathrm{vO}_2}
動脈血と混合静脈血の酸素含量(mL O₂/dL)
D˙O2,  V˙O2\dot D_{\mathrm{O}_2},\;\dot V_{\mathrm{O}_2}
酸素供給量と設定した酸素消費量(mL O₂/min)。係数10はLからdLへの換算

拍平均の酸素収支

酸素輸送は循環の解から計算する拍平均の代数評価で、血液量・張力へはフィードバックしません。Pは酸素分圧(mmHg)、Sは飽和度、Hbはg/dL、CはmL O₂/dLです。肺胞気と肺毛細管終末の酸素分圧を等しいと仮定します。PBは気圧、Rは呼吸交換比です。

COはL/min、VO₂はmL/min、sは酸素混合上のシャント率です。sは循環回路に追加する血流接続ではありません。肺毛細管終末含量Cc、静脈含量Cv、動脈含量CaをFick収支と混合から求めます。CO≤0や負の必要含量なら評価不能とし、飽和度を見栄えのよい値に置換しません。

PAO2=FIO2(PB47)PaCO2/R>0,S(P)=P2.726.82.7+P2.7P_{\mathrm{AO}_2}=F_{\mathrm{IO}_2}(P_B-47)-P_{\mathrm{aCO}_2}/R>0,\quad S(P)=\frac{P^{2.7}}{26.8^{2.7}+P^{2.7}}
C(P)=1.34HbS(P)+0.0031P,D=VO210CO,Ca=Ccs1sD,Cv=CaDC(P)=1.34\,Hb\,S(P)+0.0031P,\quad D=\frac{VO_2}{10\mathrm{CO}},\quad C_a=C_c-\frac{s}{1-s}D,\quad C_v=C_a-D

解析:負荷を変えて心臓の応答を調べる

1拍のPV loopはシミュレーションが作る軌道です。ESPVR・EDPVR・Starling/Guyton曲線やPVA/PEは、別の負荷条件を計算する解析から得ます。

数式・仮定を詳しく

解析用に状態を複製するため、操作中のシミュレーションは変わりません。収縮末期圧–容積関係(ESPVR)は低容量側からbaselineまでの複数条件で求めます。拡張末期圧–容積関係(EDPVR)とStarling曲線には高容量側も含めます。推定値は負荷範囲や各条件の定常化に依存し、ESPVRが直線になるとは限りません。曲線や面積の推定が成立しない場合は未評価と表示します。

結合と時間積分

独立した血液量は31区画に分布し、その総和をTBVに固定します(独立自由度は30)。各壁にLandの6状態・粘性歪み1状態・Ca源2状態、四弁に開口状態を持ちます。さらに冠循環の6トーンと、拍時刻・Ca入力強度・直前の僧帽弁閉鎖時歪みを引き継ぎます。心腔圧・血管圧・流量・中隔位置・接合円半径は、その時点の連立条件で求める量です。

Vn+1Vn=ΔtNQn+1,zn+1zn=Δtf(zn+1,Can+1,λn+1,λ˙n+1),0=g(Vn+1,zn+1,Pn+1,Qn+1,vS,n+1,yn+1).\begin{aligned}V_{n+1}-V_n&=\Delta t\,NQ_{n+1},\\z_{n+1}-z_n&=\Delta t\,f(z_{n+1},Ca_{n+1},\lambda_{n+1},\dot\lambda_{n+1}),\\0&=g(V_{n+1},z_{n+1},P_{n+1},Q_{n+1},v_{S,n+1},y_{n+1}).\end{aligned}

zはLand・粘弾性・開口の状態、fは各節の時間発展式、gは圧–容量関係・流量則・力の釣り合いです。全体を後退Eulerで連立し、Ca源はイベント間の指数解を使います。候補容積→形状・筋長→材料応力→圧→流量→容積収支が同時に整合するまで解きます。血管だけを先に更新して古い心腔圧を使い続ける手順ではありません。

通常刻みは2 msです。興奮・Ca入力・制御周期の境界ではステップを分けます。冠トーンは完了周期の流量積分を用いて更新します。状態の非負性・有限性・体積保存・非線形残差を満たさない試行は採用しません。別の数値積分法でも同じ連続モデルを組めますが、有限刻みのピーク・弁イベント・保存されたbaselineとの一致は別途検証が必要です。

適用範囲と限界

各区画は空間的な分布を平均した代表値です。局所ジェット、波の伝播・反射、心室内の圧分布、局所虚血や組織の不均一性をそのまま再現しません。測定位置・解析方法が異なる実測値を、同じ名称だけで直接比較することはできません。

数式を定めること、ある設定で計算が安定すること、その症例の生理的妥当性を示すことは別です。採用条件・文献との比較・操作後の検証範囲は、各プリセットに記録しています。教育・研究のためのモデルで、患者の診断や治療判断を目的としません。