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

See Also

Related Examples

Was this topic helpful?