微分方程式で数学モデルを作ろう　2章

人工腎臓器、いわゆる透析装置のモデルについてです。
前回の薬のモデルは濃度によって排出量が決まりましたが、
今回は流速によって除去率が決まるモデルです。
老廃物のある血液濃度と透析液の濃度をを変数 $u,v$ で表し、 $x$ を透析開始点からの位置とすると、
$u,v$ は $u=u(x),v=v(x)$ で表されると仮定します。

ここでは以下のフィックの法則を用いて、膜を通過する物質の総量を考えます。

フィックの法則：単位時間に単位面積の幕を通過する物質の総量は、
　　　　　　　　　　膜のその位置における濃度差に比例する

位置 $x$ から $x+\delta x$ の微小区間について
単位時間当たり、血液から透析液まで移行する量は比例定数 $k$ を用いて、
$k\{u(x)-v(x)\}\delta x$
と与えられます。
さらに、(位置 $x$ に入ってきた血液の量)=(膜を通して移行した量)+( $x+\delta x$ で出ていく血液の量)が成立していますので、
微小区間内での血液の量を $Q_{b}$ とすると、
$Q_{b}u(x)=k\{u(x)-v(x)\}\delta x+Q_{b}u(x+\delta x)$
と書き換えられます。
$\delta x \rightarrow 0$ とすると、以下の微分方程式が成り立ちます。
$Q_{b}\frac{du}{dx}=-k(u-v)$
逆に、微小区間内での透析液の量を $Q_{d}$ とすると、
$Q_{d}v(x)=k\{u(x)-v(x)\}\delta x+Q_{d}v(x+\delta x)$
から、以下の微分方程式が成り立ちます。
$Q_{d}\frac{dv}{dx}=k(u-v)$
ここで、 $z=u-v , \alpha=\frac{k}{Q_{b}}-\frac{k}{Q_{d}}$ とすると、
$\frac{du}{dx}-\frac{dv}{dx}=-(\frac{k}{Q_{b}}-\frac{k}{Q_{d}})(u-v)$
より、
$\frac{dz}{dx}=-\alpha z$
となります。
よって、 $u,v$ は任意定数 $C_{1},C_{2}$ を用いて
$u=C_{1}+\frac{kC_{2}}{\alpha Q_{b}}e^{-\alpha x}$
$v=C_{1}+\frac{kC_{2}}{\alpha Q_{d}}e^{-\alpha x}$
と表されます。
ここで、初期条件について
$x=0$ で $u=u_{0}$
$x=L$ で $v=0$ 　( $L$ は透析の開始位置)
と定めると、
$u=u_{0}[\frac{(e^{-\alpha L}/Q_{d})-(e^{-\alpha x}/Q_{b})}{(e^{-\alpha L}/Q_{d})-(1/Q_{b})}]$
$v=\frac{u_{0}}{Q_{d}}[\frac{e^{-\alpha L}-e^{-\alpha x}}{(e^{-\alpha L}/Q_{d})-(1/Q_{d})}]$
となります。

透析装置の設計者は浄化値（クリアランス） $Cl$ を考えるそうです。
浄化値とは
$Cl=\frac{\int^{L}_{0}k\{u(x)-v(x)\}dx}{u_{0}}=\frac{Q_{b}}{u_{0}}\{u_{0}-u(L)\}$
で定義される数値のことで、具体的には
$Cl=Q_{b}[\frac{1-e^{-\alpha L}}{1-(Q_{b}/Q_{d})e^{-\alpha L}}]$
と表されます。

このモデルのパラメーターは $Q_{b}/Q_{d}$ と $kL/Q_{b}$ であるそうです。
なぜこの2つなのかと考えると、可変なのが $Q_{d}$ と $L$ であり、
これらを含み、モデルを考えるうえでわかりやすい独立変数がこの二つになるからだと思われます。

for i=1:1:800;
    k=1;
    L=200;
    Qd=400;
    a=k*(50/i-1/Qd);
    Cl=(i/50)*((1-exp(-a*L))/(1-(i/50*Qd)*exp(-a*L)));
    hold on;
    plot(i/50,Cl,'b-');
    axis([0 17 0 20]);
    title('浄化値の理論曲線')
    xlabel('Qb ml/分');
    ylabel('Cl ml/分');
    set(gca,'XTick',2:2:16);
    set(gca,'YTick',2:2:20);
    hold off;
end;
saveas(figure(1),'Dialysis','png');

グラフは本とはかなり異なった形になりました。
数式的には $Q_{b}/Q_{d}$ が小さい時は $Cl$ は $Q_{b}$ の一次近似になるので、
本のグラフとは $Q_{d}$ の値が異なるのかもしれません。