具有不连续性的心血管模型 dde -凯发k8网页登录

此示例说明如何使用 dde23 对具有不连续导数的心血管模型求解。此示例最初由 ottesen [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)$ 由 dely(t,y) 的分量 j 给出。
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}]$ 并使用 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.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 object. the axes object with title heart rate for baroreflex-feedback mechanism., xlabel time t, ylabel h(t) 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.

另请参阅

| | |