Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

具有不连续性的心血管模型 DDE

此示例说明如何使用 dde23 对具有不连续导数的心血管模型求解。此示例最初由 Ottesen [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.

TsTp 的项分别是同一方程在有时滞和没有时滞状态下的变体。Pa τPa 分别代表在有时滞和没有时滞状态下的平均动脉压。

Ts=11+(Paταs)βs

Tp=11+(Paαp)-βp.

此问题有许多物理参数:

  • 动脉顺应性 ca=1.55 ml/mmHg

  • 静脉顺应性 cv=519 ml/mmHg

  • 外周阻力 R=1.05(0.84)mmHg s/ml

  • 静脉流出阻力 r=0.068 mmHg s/ml

  • 心博量 Vstr=67.9(77.9) ml

  • 典型平均动脉压 P0=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 处的低阶导数具有不连续性。

常历史解由以下物理参数定义

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)dely(t,y) 的分量 j 给出。

  • 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 

注意:所有函数都作为局部函数包含在示例的末尾。

编写历史解代码

接下来,创建一个向量来定义三个分量 PaPvH 的常历史解。历史解是时间 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] 并使用 dde23 求解器对 DDE 求解。使用匿名函数指定 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.xsol.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] Ottesen, J. T.“Modelling of the Baroreflex-Feedback Mechanism with Time-Delay.”J. Math.Biol.Vol. 36, Number 1, 1997, pp. 41–63.

另请参阅

| | |

相关主题