Cardiovascular Model DDE with Discontinuities

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

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.

The terms for Ts and Tp are variations of the same equation with and without time delay. Paτ and Pa represent the mean arterial pressure with and without time delay, respectively.

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

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

This problem has a number of physical parameters:

  • Arterial compliance ca=1.55ml/mmHg

  • Venous compliance cv=519ml/mmHg

  • Peripheral resistance R=1.05(0.84)mmHgs/ml

  • Venous outflow resistance r=0.068mmHgs/ml

  • Stroke volume Vstr=67.9(77.9)ml

  • Typical mean arterial pressure P0=93mmHg

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

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

Pa=P0,        Pv(t)=11+RrP0,         H(t)=1RVstr11+rRP0.

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 Paτ(t)=Pa(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 yn(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)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 

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 Pa, Pv, and H. The solution history is the solution for times 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];

Solve Equation

Use ddeset to specify the presence of the discontinuity at t=600. Finally, define the interval of integration [t0  tf] 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)')

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.

See Also

| | |

Related Topics