# Documentation

## Cardiovascular Model with Discontinuities

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``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.

1. 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.

2. 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;```
3. 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]; ```
4. Define the delay, `tau`. Add this code to the main function.

`tau = 4;`
5. 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.

6. 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.

7. 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)')```
8. 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.

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

## References

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