불연속을 가지는 심혈관 모델 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}$

• 박출량 ${\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$에서 저계 도함수(Low-order Derivative)에 불연속을 가집니다.

상수 해 내역은 다음과 같은 물리적 파라미터로 정의됩니다.

$${\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는 파라미터 값을 전달하는 데 사용되는 네 번째 선택적 입력값입니다.

처음 세 개의 입력값은 솔버에 의해 자동으로 함수에 전달되며 변수 이름에 따라 방정식 코딩 방식이 결정됩니다. 파라미터 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 

참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.

해 내역 코딩하기

다음으로 세 개의 성분 ${\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를 풉니다. 파라미터 구조체 p에서 전달하는 익명 함수를 사용하여 ddefun을 지정합니다.

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을 사용하여 특정 지점의 해를 구할 수 있습니다.)

시간에 대해 세 번째 해 성분(심박수)을 플로팅합니다.

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

로컬 함수(Local Function)

여기 나열된 함수는 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

관련 항목

• 지연 미분 방정식 풀기
• 중립(Neutral) 유형의 DDE