Haskell による血圧のシミュレーション

Haskell で左室系の血圧と心収縮力の関係をシミュレーションするプログラムを作ってみた。

動脈を流れる血圧と血流量の関係は、大まかには、空気室モデルでモデル化できる。

左心室から収縮によって大動脈へ流れ込んだ血液は、いったん大動脈の弾性によって動脈内に貯蔵され、心室が拡張して血液が大動脈に流れ込まなくなったら、大動脈の弾性力によって、末梢へ血液が供給される。また、末梢では一定の末梢血管抵抗を受けながら、血液は動脈内から静脈系へ流れ出していく。

したがって、大動脈の血圧と血流の関係は電気回路のRC回路でモデル化できる。このとき、血流と血圧の間のインパルス応答(インピーダンスのインパルス応答)を計算することができるので、後で述べる心室の収縮力とのインターフェースを考えることができる。RC回路のインピーダンスのインパルス応答関数は、P = Po * exp (-1/RC * t) という指数関数になる。

一方、左心室の計量的な性質は P = V * E(t) という簡潔な関係式であらわすことができる。心室内の圧力 P と心室内容積の比 P / V は心筋の張力が一定であれば常に一定である。この関係式は心筋の聴力や収縮力を数量的に表わすことができる。これを大動脈のインピーダンスにつなげてやると、心筋の収縮力や、大動脈のコンプライアンスや末梢血管抵抗を変化させた時に、血圧や血流や心室内の容積がどのように変わるかをシミュレーションできる。

動脈のインピーダンスは線形だが、心室の収縮力と圧力の関係は非線形なので、プログラムしないと両者をつなげて考えることができない。Octave は線形システムの計算は簡単にできるが、プログラム環境が貧弱なので使えなかった。Haskell のリスト操作機能が強力なのでこのような問題に使えるのではないかと考えて今朝からやってみたら半日でできた。

下がそのプログラムだ。左心室と大動脈のインターフェースは、大動脈弁があったり、左心室圧と大動脈圧の二つの圧力があるため難しいが、今回は単純に大動脈入口部に抵抗成分があると考え、左心室からの血液の流入は左心室圧と動脈圧の差をこの抵抗で割ったものが流れると仮定した(乱暴な仮定だ)。また、大動脈弁があるので血液は大動脈から左心室へ向かうことはできないので、左室圧が大動脈圧より低い時は血流が0になるようにプログラムした。

このように、左心室と大動脈の結合は複雑なので単純な線形システムでは記述できない。

さて、プログラムの方は下のようになるが、あまり分かりやすくは説明できないので、結果を、Excel でグラフにしたものを表示した。Haskell だとこういうことがわりに簡単にできる。いまのところ、思考補助手段としてのプログラム言語は Haskell が最右翼だ。

1.左心系の血圧シミュレーションプログラム

main = writeFile "bpdata.csv" $ init $ tail $ show $ map fst result

-- Blood pressure of the first 3 contractions
result = first ++ second ++ third

first = bp registance inductance contraction [0.0] 1.0
second = bp registance inductance contraction (reverse (map snd first)) 1.0
third = bp registance inductance contraction (reverse (map snd (first ++ second))) 1.0

-- The contractiliti - pressure interaction model
bp :: Float -> [Float] -> [Float] -> [Float] -> Float -> [(Float,Float)]
bp r hs (e:es) fs v =
  (pa, f):bp r hs es (f:fs) v'
    where v' = v - head fs
          pv = v * e
          pa = pressure hs fs
          f | pv > pa = (pv - pa) / r
            | otherwise = 0

bp _ _ _ _ _ = []

-- Parameters

registance :: Float
registance = 4.0

pressure :: [Float] -> [Float] -> Float
pressure hs ps = sum (zipWith (*) hs ps)

inductance :: [Float]
inductance = map (\x -> exp (-0.2*x) :: Float) [0..20]

contraction :: [Float]
contraction = replicate 4 1.0 ++ replicate 6 0.0

2.上のプログラムで計算した血圧のシミュレーション

d0038298_1515665.jpg

[PR]
by tnomura9 | 2009-09-27 15:19 | Haskell | Comments(0)
<< Simple, simple ... 費用対効果比 >>