不連続をもつ心臓血管モデル DDE

この例では、dde23 を使用して、不連続導関数をもつ心臓血管モデルを解く方法を説明します。この例は、元々 Ottesen [1] により提示されました。

方程式系は次のとおりです。

$$\dot{{\mathit{P}}_{\mathit{a}}}\left(\mathit{t}\right)=-\frac{1}{{\mathit{c}}_{\mathit{a}}\mathit{R}}{\mathit{P}}_{\mathit{a}}\left(\mathit{t}\right)+\frac{1}{{\mathit{c}}_{\mathit{a}}\mathit{R}}{\mathit{P}}_{\mathit{v}}\left(\mathit{t}\right)+\frac{1}{{\mathit{c}}_{\mathit{a}}}{\mathit{V}}_{\mathrm{str}}\left({\mathit{P}}_{\mathit{a}}^{\tau }\left(\mathit{t}\right)\right)\mathit{H}\left(\mathit{t}\right)$$

$${\dot{\mathit{P}}}_{\mathit{v}}\left(\mathit{t}\right)=\frac{1}{{\mathit{c}}_{\mathit{v}}\mathit{R}}{\mathit{P}}_{\mathit{a}}\left(\mathit{t}\right)-\left(\frac{1}{{\mathit{c}}_{\mathit{v}}\mathit{R}}+\frac{1}{{\mathit{c}}_{\mathit{v}}\mathit{r}}\right){\mathit{P}}_{\mathit{v}}\left(\mathit{t}\right)$$

$$\dot{\mathit{H}}\left(\mathit{t}\right)=\frac{{\alpha }_{\mathit{H}}{\mathit{T}}_{\mathit{s}}}{1+{\gamma }_{\mathit{H}}{\mathit{T}}_{\mathit{p}}}-{\beta }_{\mathit{H}}{\mathit{T}}_{\mathit{p}}.$$

${\mathit{T}}_{\mathit{s}}$ および ${\mathit{T}}_{\mathit{p}}$ の項は、時間遅延の有無が異なる同じ方程式です。${\mathit{P}}_{\mathit{a}}^{\tau }$ および ${\mathit{P}}_{\mathit{a}}$ は、それぞれ遅延を伴う場合と伴わない場合の平均動脈圧です。

$${\mathit{T}}_{\mathit{s}}=\frac{1}{1+{\left(\frac{{\mathit{P}}_{\mathit{a}}^{\tau }}{{\alpha }_{\mathit{s}}}\right)}^{{\beta }_{\mathit{s}}}}$$

$${\mathit{T}}_{\mathit{p}}=\frac{1}{1+{\left(\frac{{\mathit{P}}_{\mathit{a}}}{{\alpha }_{\mathit{p}}}\right)}^{-{\beta }_{\mathit{p}}}}.$$

この問題には、次のような多くの物理パラメーターがあります。

• 動脈コンプライアンス ${\mathit{c}}_{\mathit{a}}=1.55\mathrm{ml}/\mathrm{mmHg}$

• 静脈コンプライアンス ${\mathit{c}}_{\mathit{v}}=519\mathrm{ml}/\mathrm{mmHg}$

• 末梢抵抗 $\mathit{R}=1.05\left(0.84\right)\mathrm{mmHg}\mathit{s}/\mathrm{ml}$

• 静脈流出量抵抗 $\mathit{r}=0.068\mathrm{mmHg}\mathit{s}/\mathrm{ml}$

• 1 回拍出量 ${\mathit{V}}_{\mathrm{str}}=67.9\left(77.9\right)\mathrm{ml}$

• 一般的な平均動脈圧 ${\mathit{P}}_0=93\mathrm{mmHg}$

• $${\alpha }_0={\alpha }_{\mathit{s}}={\alpha }_{\mathit{p}}=93\left(121\right)\mathrm{mmHg}$$

• $${\alpha }_{\mathit{H}}=0.84{\sec }^{-2}$$

• $${\beta }_0={\beta }_{\mathit{s}}={\beta }_{\mathit{p}}=7$$

• $${\beta }_{\mathit{H}}=1.17$$

• $${\gamma }_{\mathit{H}}=0$$

この方程式系は末梢血圧によって多大な影響を受けますが、この血圧は、$\mathit{t}=600$ を始点として $\mathit{R}=1.05$ から $\mathit{R}=0.84$ に指数関数的に減少します。その結果、この方程式系には $\mathit{t}=600$ において低次導関数に不連続性があります。

定数解履歴は物理パラメーターで定義されます。

$${\mathit{P}}_{\mathit{a}}={\mathit{P}}_0,{\mathit{P}}_{\mathit{v}}\left(\mathit{t}\right)=\frac{1}{1+\frac{\mathit{R}}{\mathit{r}}}{\mathit{P}}_0,\mathit{H}\left(\mathit{t}\right)=\frac{1}{{\mathrm{RV}}_{\mathrm{str}}}\frac{1}{1+\frac{\mathit{r}}{\mathit{R}}}{\mathit{P}}_0.$$

MATLAB® でこの方程式系を解くには、定数の時間遅延をもつ方程式系に適した遅延微分方程式ソルバー dde23 を呼び出す前に、方程式、パラメーター、遅延、および履歴をコード化する必要があります。(ここで行ったように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。

物理パラメーターの定義

最初に、問題の物理パラメーターを構造体のフィールドとして定義します。

p.ca     = 1.55;
p.cv     = 519;
p.R      = 1.05;
p.r      = 0.068;
p.Vstr   = 67.9;
p.alpha0 = 93;
p.alphas = 93;
p.alphap = 93;
p.alphaH = 0.84;
p.beta0  = 7;
p.betas  = 7;
p.betap  = 7;
p.betaH  = 1.17;
p.gammaH = 0;

遅延のコード化

次に、項 ${\mathit{P}}_{\mathit{a}}^{\tau }\left(\mathit{t}\right)={\mathit{P}}_{\mathit{a}}\left(\mathit{t}-\tau \right)$ の方程式で、定数の時間遅延 $\tau$ を表す変数 tau を作成します。

tau = 4;

方程式のコード作成

ここで、方程式をコード化する関数を作成します。この関数にはシグネチャ dydt = ddefun(t,y,Z,p) がなければなりません。ここでは以下のようになります。

• t は時間 (独立変数) です。

• y は解 (従属変数) です。

• Z(n,j) は遅延 ${\mathit{y}}_{\mathit{n}}\left(\mathit{d}\left(\mathit{j}\right)\right)$ を近似します。ここで、遅延 $\mathit{d}\left(\mathit{j}\right)$ は dely(t,y) の要素 j によって与えられます。

• p は、パラメーターの値を渡すために必要に応じて使用される 4 番目の入力です。

最初の 3 つの入力はソルバーによって関数に自動的に渡され、変数名によって方程式のコード化方法が決まります。パラメーター p の構造体は、ソルバーを呼び出すときに関数に渡されます。この場合、遅延は次で表されます。

• Z(:,1)$\to {\mathit{P}}_{\mathit{a}}\left(\mathit{t}-\tau \right)$

function dydt = ddefun(t,y,Z,p)
    if t <= 600
      p.R = 1.05;
    else
      p.R = 0.21 * exp(600-t) + 0.84;
    end    
    ylag = Z(:,1);
    Patau = ylag(1);
    Paoft = y(1);
    Pvoft = y(2);
    Hoft  = y(3);

    dPadt = - (1 / (p.ca * p.R)) * Paoft ...
            + (1/(p.ca * p.R)) * Pvoft ...
            + (1/p.ca) * p.Vstr * Hoft;

    dPvdt = (1 / (p.cv * p.R)) * Paoft...
            - ( 1 / (p.cv * p.R)...
            + 1 / (p.cv * p.r) ) * Pvoft;

    Ts = 1 / ( 1 + (Patau / p.alphas)^p.betas );
    Tp = 1 / ( 1 + (p.alphap / Paoft)^p.betap );

    dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp) ...
           - p.betaH * Tp;

    dydt = [dPadt; dPvdt; dHdt];
end 

メモ: 関数はすべて例の最後にローカル関数として含まれます。

解の履歴のコード化

次に、ベクトルを作成して、3 つの要素 ${\mathit{P}}_{\mathit{a}}$、${\mathit{P}}_{\mathit{v}}$、および $\mathit{H}$ の定数解履歴を定義します。解の履歴は時間 $\mathit{t}\le {\mathit{t}}_0$ での解です。

P0 = 93;
Paval = P0;
Pvval = (1 / (1 + p.R/p.r)) * P0;
Hval = (1 / (p.R * p.Vstr)) * (1 / (1 + p.r/p.R)) * P0;
history = [Paval; Pvval; Hval];

方程式の求解

ddeset を使用して、$\mathit{t}=600$ での不連続性の存在を指定します。最後に、積分区間 $\left[{\mathit{t}}_0{\mathit{t}}_{\mathit{f}}\right]$ を定義し、dde23 ソルバーを使用して DDE を解きます。無名関数を使用して ddefun を指定し、パラメーターの構造体 p を渡します。

options = ddeset('Jumps',600);
tspan = [0 1000];
sol = dde23(@(t,y,Z) ddefun(t,y,Z,p), tau, history, tspan, options);

解のプロット

解の構造体 sol には sol.x フィールドと sol.y フィールドがあり、ソルバーが受け入れる内部タイム ステップと、その時間に対応する解が含まれます (特定の点での解が必要な場合、deval を使用して特定の点での解を評価することができます)。

3 番目の解要素 (心拍数) を時間に対してプロットします。

plot(sol.x,sol.y(3,:))
title('Heart Rate for Baroreflex-Feedback Mechanism.')
xlabel('Time t')
ylabel('H(t)')

ローカル関数

ここでは、DDE ソルバー dde23 が解を計算するために呼び出すローカル補助関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。

function dydt = ddefun(t,y,Z,p) % equation being solved
    if t <= 600
      p.R = 1.05;
    else
      p.R = 0.21 * exp(600-t) + 0.84;
    end    
    ylag = Z(:,1);
    Patau = ylag(1);
    Paoft = y(1);
    Pvoft = y(2);
    Hoft  = y(3);

    dPadt = - (1 / (p.ca * p.R)) * Paoft ...
            + (1/(p.ca * p.R)) * Pvoft ...
            + (1/p.ca) * p.Vstr * Hoft;

    dPvdt = (1 / (p.cv * p.R)) * Paoft...
            - ( 1 / (p.cv * p.R)...
            + 1 / (p.cv * p.r) ) * Pvoft;

    Ts = 1 / ( 1 + (Patau / p.alphas)^p.betas );
    Tp = 1 / ( 1 + (p.alphap / Paoft)^p.betap );

    dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp) ...
           - p.betaH * Tp;

    dydt = [dPadt; dPvdt; dHdt];
end 

参照

[1] Ottesen, J. T."Modelling of the Baroreflex-Feedback Mechanism with Time-Delay."J. Math.Biol.Vol. 36, Number 1, 1997, pp. 41–63.

参考

ddensd | ddesd | dde23 | deval

関連するトピック

• 遅延微分方程式の求解
• 中立型 DDE