This example shows how to use dde23 to solve a cardiovascular model that has a discontinuous derivative as presented by Ottesen .
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.
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)')
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.
 Ottesen, J. T. "Modelling of the Baroflex-Feedback Mechanism with Time-Delay." J. Math. Biol. Vol. 36, Number 1, 1997, pp. 41–63.