Сердечно-сосудистый образцовый DDE с разрывами

Этот пример показывает, как использовать 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 {секунда}^{- 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)')

Localfunctions

Перечисленный здесь локальные функции помощника, которые решатель 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] Оттесен, J. T. “Моделирование Механизма Baroreflex-обратной-связи с Задержкой”. J. Математика. Biol. Издание 36, Номер 1, 1997, стр 41–63.

Смотрите также

dde23 | ddensd | ddesd | deval

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