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 .
The system of equations is
The terms for and are variations of the same equation with and without time delay. and represent the mean arterial pressure with and without time delay, respectively.
This problem has a number of physical parameters:
Venous outflow resistance
Typical mean arterial pressure
The system is heavily influenced by peripheral pressure, which decreases exponentially from to beginning at . As a result, the system has a discontinuity in a low-order derivative at .
The constant solution history is defined in terms of the physical parameters
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;
Next, create a variable
tau to represent the constant time delay in the equations for the terms .
tau = 4;
Now, create a function to code the equations. This function should have the signature
dydt = ddefun(t,y,Z,p), where:
tis time (independent variable).
yis the solution (dependent variable).
Z(n,j)approximates the delays , where the delay is given by component
pis 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:
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 , , and . The solution history is the solution for times .
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 to specify the presence of the discontinuity at . Finally, define the interval of integration and solve the DDE using the
dde23 solver. Specify
ddefun using an anonymous function to pass in the structure of parameters,
options = ddeset('Jumps',600); tspan = [0 1000]; sol = dde23(@(t,y,Z) ddefun(t,y,Z,p), tau, history, tspan, options);
The solution structure
sol has the fields
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)')
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
 Ottesen, J. T. “Modelling of the Baroreflex-Feedback Mechanism with Time-Delay.” J. Math. Biol. Vol. 36, Number 1, 1997, pp. 41–63.