Сердечно-сосудистая модель DDE с разрывами

Открыть Live Script

В этом примере показано, как использовать dde23 решить сердечно-сосудистую модель, которая имеет прерывистое производное. Пример был первоначально представлен Оттесеном [1].

Система уравнений

$\dot{P_{a}} (t) = - \frac{1}{c_{a} R} P_{a} (t) + \frac{1}{c_{a} R} P_{v} (t) + \frac{1}{c_{a}} V_{str} (P_{a}^{τ} (t)) H (t)$

${\dot{P}}_{v} (t) = \frac{1}{c_{v} R} P_{a} (t) - (\frac{1}{c_{v} R} + \frac{1}{c_{v} r}) P_{v} (t)$

$\dot{H} (t) = \frac{α_{H} T_{s}}{1 + γ_{H} T_{p}} - β_{H} T_{p} .$

Условия для $T_{s}$ и $T_{p}$ являются изменениями того же уравнения с временной задержкой и без нее. $P_{a}^{τ}$ и $P_{a}$ представляют среднее артериальное давление с задержкой и без нее, соответственно.

$T_{s} = \frac{1}{1 + {(\frac{P_{a}^{τ}}{α_{s}})}^{β_{s}}}$

$T_{p} = \frac{1}{1 + {(\frac{P_{a}}{α_{p}})}^{- β_{p}}} .$

Эта задача имеет ряд физических параметров:

Артериальная податливость $c_{a} = 1.55 ml / mmHg$
Венозная податливость $c_{v} = 519 ml / mmHg$
Периферийное сопротивление $R = 1.05 (0.84) mmHg s / ml$
Сопротивление венозному оттоку $r = 0.068 mmHg s / ml$
Объем штриха $V_{str} = 67.9 (77.9) ml$
Типичное среднее артериальное давление $P_{0} = 93 mmHg$
$α_{0} = α_{s} = α_{p} = 93 (121) mmHg$
$α_{H} = 0.84 \sec^{- 2}$
$β_{0} = β_{s} = β_{p} = 7$
$β_{H} = 1.17$
$γ_{H} = 0$

На систему сильно влияет периферийное давление, которое уменьшается экспоненциально $R = 1.05$ кому $R = 0.84$ начиная с $t = 600$ . В результате система имеет разрыв в производной низкого порядка $t = 600$ .

История постоянных решений определяется в терминах физических параметров

$P_{a} = P_{0}, P_{v} (t) = \frac{1}{1 + \frac{R}{r}} P_{0}, H (t) = \frac{1}{{RV}_{str}} \frac{1}{1 + \frac{r}{R}} 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;

Задержка кода

Затем создайте переменную tau для представления постоянной временной задержки $τ$ в уравнениях для членов $P_{a}^{τ} (t) = P_{a} (t - τ)$ .

tau = 4;

Кодовые уравнения

Теперь создайте функцию, чтобы кодировать уравнения. Эта функция должна иметь подпись dydt = ddefun(t,y,Z,p), где:

t является временем (независимая переменная).
y - решение (зависимая переменная).
Z(n,j) аппроксимирует задержки $y_{n} (d (j))$ , где задержка $d (j)$ задается как компонент j от dely(t,y).
p является необязательным четвертым входом, используемым для прохождения в значениях параметров.

Первые три входов автоматически передаются функции решателем, и имена переменных определяют, как вы кодируете уравнения. Структура параметров p передается в функцию, когда вы вызываете решатель. В этом случае задержки представлены:

Z(:,1) $\to P_{a} (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

Примечание.Все функции включены в качестве локальных функций в конце примера.

История решения кода

Затем создайте вектор, чтобы задать историю постоянных решений для трех компонентов $P_{a}$ , $P_{v}$ , и $H$ . История решений является решением для времени $t \leq 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 для определения наличия разрыва в $t = 600$ . Наконец, задайте интервал интегрирования $[t_{0} t_{f}]$ и решить DDE используя dde23 решатель. Задайте 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 для оценки решения в конкретных точках.)

Постройте график компонента третьего решения (частота сердечных сокращений) относительно времени.

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

Figure contains an axes. The axes with title Heart Rate for Baroreflex-Feedback Mechanism. 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] Оттесен, Дж. Т. «Моделирование механизма Baroreflex-Feedback с задержкой по времени». J. Math. Biol. vol. 36, Number 1, 1997, pp. 41-63.

См. также

dde23 | ddensd | ddesd | deval

Документация