Cardiovascular Model DDE with Discontinuities

Open Live Script

This example shows how to use dde23 to solve a cardiovascular model that has a discontinuous derivative. The example was originally presented by Ottesen [1].

The system of equations is

$\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} .$

The terms for $T_{s}$ and $T_{p}$ are variations of the same equation with and without time delay. $P_{a}^{τ}$ and $P_{a}$ represent the mean arterial pressure with and without time delay, respectively.

$T_{s} = \frac{1}{1 + {(\frac{P_{a}^{τ}}{α_{s}})}^{β_{s}}}$

$T_{p} = \frac{1}{1 + {(\frac{P_{a}}{α_{p}})}^{- β_{p}}} .$

This problem has a number of physical parameters:

Arterial compliance $c_{a} = 1.55 ml / mmHg$
Venous compliance $c_{v} = 519 ml / mmHg$
Peripheral resistance $R = 1.05 (0.84) mmHg s / ml$
Venous outflow resistance $r = 0.068 mmHg s / ml$
Stroke volume $V_{str} = 67.9 (77.9) ml$
Typical mean arterial pressure $P_{0} = 93 mmHg$
$α_{0} = α_{s} = α_{p} = 93 (121) mmHg$
$α_{H} = 0.84 \sec^{- 2}$
$β_{0} = β_{s} = β_{p} = 7$
$β_{H} = 1.17$
$γ_{H} = 0$

The system is heavily influenced by peripheral pressure, which decreases exponentially from $R = 1.05$ to $R = 0.84$ beginning at $t = 600$ . As a result, the system has a discontinuity in a low-order derivative at $t = 600$ .

The constant solution history is defined in terms of the physical parameters

$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} .$

To solve this system of equations in MATLAB®, you need to code the equations, parameters, delays, and history before calling the delay differential equation solver dde23, which is meant for systems with constant time delays. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Define Physical Parameters

First, define the physical parameters of the problem as fields in a structure.

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;

Code Delay

Next, create a variable tau to represent the constant time delay $τ$ in the equations for the terms $P_{a}^{τ} (t) = P_{a} (t - τ)$ .

tau = 4;

Code Equations

Now, create a function to code the equations. This function should have the signature dydt = ddefun(t,y,Z,p), where:

t is time (independent variable).
y is the solution (dependent variable).
Z(n,j) approximates the delays $y_{n} (d (j))$ , where the delay $d (j)$ is given by component j of dely(t,y).
p is an optional fourth input used to pass in the values of the parameters.

The first three inputs are automatically passed to the function by the solver, and the variable names determine how you code the equations. The structure of parameters p is passed to the function when you call the solver. In this case the delays are represented with:

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

Note: All functions are included as local functions at the end of the example.

Code Solution History

Next, create a vector to define the constant solution history for the three components $P_{a}$ , $P_{v}$ , and $H$ . The solution history is the solution for times $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];

Solve Equation

Use ddeset to specify the presence of the discontinuity at $t = 600$ . Finally, define the interval of integration $[t_{0} t_{f}]$ and solve the DDE using the dde23 solver. Specify ddefun using an anonymous function to pass in the structure of parameters, p.

options = ddeset('Jumps',600);
tspan = [0 1000];
sol = dde23(@(t,y,Z) ddefun(t,y,Z,p), tau, history, tspan, options);

Plot Solution

The solution structure sol has the fields sol.x and sol.y that contain the internal time steps taken by the solver and corresponding solutions at those times. (If you need the solution at specific points, you can use deval to evaluate the solution at the specific points.)

Plot the third solution component (heart rate) against time.

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.

Local Functions

Listed here are the local helper functions that the DDE solver dde23 calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

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

References

[1] Ottesen, J. T. “Modelling of the Baroreflex-Feedback Mechanism with Time-Delay.” J. Math. Biol. Vol. 36, Number 1, 1997, pp. 41–63.