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

Открыть сценарий в реальном времени

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

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

$\overset{Pa˙}{_{}} (t) \frac{=}{_{-}}_{} 1caRPa \frac{(}{t_{)}} +_{} \frac{}{{1caRPv}_{}} (_{t)} +_{}^{} 1caVstr (P$ α (t)) H (t)

${\overset{}{}}_{P˙v} (t) \frac{}{=_{}}_{} 1cvRPa (\frac{t}{)_{} -} \frac{(}{_{}} 1cvR_{+}$ 1cvr) Pv (t)

$\overset{H˙}{} (t) \frac{_{=}_{}}{_{αHTs1}_{+}}_{}_{}$ γ HTp-βHTp.

Члены для $_{Ts}$ и $_{Tp}$ являются вариациями одного и того же уравнения с временной задержкой и без нее. $_{Pa}^{startи}$ $_{Pa}$ представляют среднее артериальное давление с задержкой по времени и без нее, соответственно.

$_{Ts} \frac{=}{{11 \frac{_{+}^{(}}{_{Pa}}}^{_{}}}$ αs) βs

$_{Tp} \frac{=}{{11 \frac{_{+}}{_{(}}}^{_{Paαp}}})$ -βp.

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

Артериальное соответствие $_{ca} = 1,55$ мл/мм рт.ст.
Венозное соответствие $_{cv} = 519$ мл/мм рт.ст.
Периферическое сопротивление $R = 1,05 (0,84) мм рт.ст./мл$
Сопротивление венозному оттоку $r = 0,068 мм рт.ст./мл$
Ударный объем $_{Vstr} = 67,9 (77,9)$ мл
Типичное среднее артериальное давление $_{P0} = 93$ мм рт.ст.
$_{α0} =_{} {αs}_{} = αp = 93 ($ 121) мм рт.ст.
$_{αH} = {0,84}^{}$ с-2
$_{β0} =_{} {βs}_{} =$ βp = 7
$_{βH} =$ 1,17
$_{γ H} =$ 0

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

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

$_{Pa} =_{} P0, Pv (t)_{} = 11 \frac{}{+ \frac{}{}}_{RrP0}, H (t) = \frac{}{{1RVstr11}_{+}} \frac{}{\frac{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;

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

Далее создайте переменную tau , чтобы представить постоянную временную задержку $в$ уравнениях для членов $_{Pa (}^{t)}_{=} Pa ($ t-start).

tau = 4;

Уравнения кода

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

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

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

Z(:,1) $_{→Pa} (t-start$ )

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

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

Журнал кодовых решений

Затем создайте вектор для определения постоянной истории решения для трех компонентов $_{Pa}$ , $_{Pv}$ и H. История решения является решением для времени $_{t≤t0}$ .

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$ ] и решите 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 с задержкой по времени». J. Math. Biol. Том 36, номер 1, 1997, стр. 41-63.

См. также

dde23 | ddensd | ddesd | deval

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