%% 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
% \dot{P_a } \left(t\right)=-\frac{1}{c_a R}P_a \left(t\right)+\frac{1}{c_a R}P_v \left(t\right)+\frac{1}{c_a }V_{\textrm{str}} \left({P\;}_a^{\tau } \left(t\right)\right)H\left(t\right)
% {\dot{P} }_v \left(t\right)=\frac{1}{c_v R}P_a \left(t\right)-\left(\frac{1}{c_v R}+\frac{1}{c_v r}\right)P_v \left(t\right)
% \dot{H} \left(t\right)=\frac{\alpha_H T_s }{1+\gamma_H T_p }-\beta_H T_p \ldotp
% The terms for T_s and T_p are variations of the same equation with and without time delay. P_a^{\;\tau } and P_a represent the mean arterial pressure with and without time delay, respectively.
% T_s =\frac{1}{1+{\left(\frac{P_a^{\;\tau } }{\alpha_s }\right)}^{\beta_s } }
% T_p =\frac{1}{1+{\left(\frac{P_a }{\alpha_p }\right)}^{-\beta_p } }\ldotp
% This problem has a number of physical parameters:
% Arterial compliance c_a =1\ldotp 55\;\textrm{ml}/\textrm{mmHg}
% Venous compliance c_v =519\;\textrm{ml}/\textrm{mmHg}
% Peripheral resistance R=1\ldotp 05\left(0\ldotp 84\right)\textrm{mmHg}\;s/\textrm{ml}
% Venous outflow resistance r=0\ldotp 068\;\textrm{mmHg}\;s/\textrm{ml}
% Stroke volume V_{\textrm{str}} =67\ldotp 9\left(77\ldotp 9\right)\;\textrm{ml}
% Typical mean arterial pressure P_0 =93\;\textrm{mmHg}
% \alpha_0 =\alpha_s =\alpha_p =93\left(121\right)\;\textrm{mmHg}
% \alpha_H =0\ldotp 84\;\sec^{-2}
% \beta_0 =\beta_s =\beta_p =7
% \beta_H =1\ldotp 17
% \gamma_H =0
% The system is heavily influenced by peripheral pressure, which decreases exponentially from R=1\ldotp 05 to R=0\ldotp 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
% P_a =P_0 ,\;\;\;\;\;\;\;\;P_v \left(t\right)=\frac{1}{1+\frac{R}{r}}P_0 ,\;\;\;\;\;\;\;\;\;H\left(t\right)=\frac{1}{{\textrm{RV}}_{\textrm{str}} }\frac{1}{1+\frac{r}{R}}P_0 \ldotp
% 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 \tau in the equations for the terms P_a^{\;\tau } \left(t\right)=P_a \left(t-\tau \right).

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 y_n \left(d\left(j\right)\right), where the delay d\left(j\right) 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)\;\to P_a \left(t-\tau \right)
% 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 P_a, P_v, and H. The solution history is the solution for times t\le t_0.

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 \left\lbrack t_0 \;\;t_f \right\rbrack 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.
% Copyright 2018 The MathWorks, Inc.
%
%