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

В этом примере показано, как использовать dde23 решить сердечно-сосудистую модель, которая имеет прерывистую производную. Пример был первоначально представлен Оттесеном [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.55ml/mmHg

  • Венозная податливость cv=519ml/mmHg

  • Периферийное сопротивление R=1.05(0.84)mmHgs/ml

  • Венозное сопротивление оттока r=0.068mmHgs/ml

  • Ударный объем Vstr=67.9(77.9)ml

  • Типичное среднее артериальное давление P0=93mmHg

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

  • αH=0.84sec-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+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;

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

Затем создайте переменную tau представлять постоянную задержку τ в уравнениях для условий Paτ(t)=Pa(t-τ).

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-τ)

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. История решения является решением в течение многих времен 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] и решите 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)')

Локальные функции

Перечисленный здесь локальные функции помощника что решатель 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.

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

| | |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте