This example shows how to use `dde23`

to
solve a cardiovascular model that has a discontinuous derivative as
presented by Ottesen [1].

Click `ddex2.m`

or
type `edit ddex2.m`

in a command window to view the
code for this example in an editor.

This is a problem with 1 delay, constant history, and 3 differential
equations with 14 physical parameters. The system is heavily influenced
by peripheral pressure, *R*, which decreases exponentially
from 1.05 to 0.84, beginning at *t* = 600. As a result,
the system has a discontinuity in a low-order derivative at *t * =
600.

Create a new program file in the editor. This file will contain a main function and a nested function. The main function accepts no inputs and returns no outputs.

Define the physical parameters. Add this code to the main function.

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;

Define the solution history. Add this code to the main function.

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];

Define the delay,

`tau`

. Add this code to the main function.tau = 4;

Define the location of discontinuity, which occurs at

*t*= 600. Add this code to the main function.`options = ddeset('Jumps',600);`

When your DDE has discontinuities in low-order derivatives, and you know the locations in advance, it is better to use

`ddeset`

with the`Jumps`

property.Solve the DDE over the interval [0, 1000]. Add this code to the main function.

sol = dde23(@ddex2de,tau,history,[0,1000],options);

The function,

`@ddex2de`

, which defines the system of DDEs, is the first input argument. You define this function in 8.Plot the solution. Add this code to the main function.

figure plot(sol.x,sol.y(3,:)) title('Heart Rate for Baroflex-Feedback Mechanism.') xlabel('time t') ylabel('H(t)')

Define the system of DDEs as a nested function inside the main function.

function dydt = ddex2de(t,y,Z) 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

This function is nested so that the main function can access the 14 parameters defined in 2.

Run your program to calculate the solution and display the plot.

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

Was this topic helpful?