Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

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

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

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

Pa˙(t)=-1caRPa(t)+1caRPv(t)+1caVstr(Paτ(t))H(t)

P˙v(t)=1cvRPa(t)-(1cvR+1cvr)Pv(t)

H˙(t)=αHTs1+γHTp-βHTp.

Ts および Tp の項は、時間遅延の有無が異なる同じ方程式です。Pa τ および Pa は、それぞれ遅延を伴う場合と伴わない場合の平均動脈圧です。

Ts=11+(Paταs)βs

Tp=11+(Paαp)-βp.

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

  • 動脈コンプライアンス ca=1.55 ml/mmHg

  • 静脈コンプライアンス cv=519 ml/mmHg

  • 末梢抵抗 R=1.05(0.84)mmHg s/ml

  • 静脈流出量抵抗 r=0.068 mmHg s/ml

  • 1 回拍出量 Vstr=67.9(77.9) ml

  • 一般的な平均動脈圧 P0=93 mmHg

  • α0=αs=αp=93(121) mmHg

  • αH=0.84 sec-2

  • β0=βs=βp=7

  • βH=1.17

  • γH=0

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

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

Pa=P0,        Pv(t)=11+RrP0,         H(t)=1RVstr11+rRP0.

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;

遅延のコード化

次に、項 Pa τ(t)=Pa(t-τ) の方程式で、定数の時間遅延 τ を表す変数 tau を作成します。

tau = 4;

方程式のコード作成

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

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

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

  • Z(n,j) は遅延 yn(d(j)) を近似します。ここで、遅延 d(j)dely(t,y) の要素 j によって与えられます。

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

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

  • Z(:,1) Pa(t-τ)

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 つの要素 PaPv、および H の定数解履歴を定義します。解の履歴は時間 tt0 での解です。

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 を使用して、t=600 での不連続性の存在を指定します。最後に、積分区間 [t0 tf] を定義し、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)')

Figure contains an axes object. The axes object with title Heart Rate for Baroreflex-Feedback Mechanism., xlabel Time t, ylabel H(t) contains an object of type line.

ローカル関数

ここでは、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.

参考

| | |

関連するトピック