Changing numeric variables based on symbolic variables in a set of differential equations.

I am implementing a mathematical model of the cardiovascular system based on a paper by Burkhoff and Tyberg.

Below is my code, with screenshots of the equations in the paper that I am attempting to represent. Any variable that is not symbolic is a constant parameter that is set at the start of the code (not shown).

%create symbols

syms V_lv(t) V_rv(t) V_C_as(t) V_C_vs(t) V_C_ap(t) V_C_vp(t) e(t)

tAll  = [];
yAll  = [];

tPool = [0 2*T_es T]; % time periods in the cardiac cycle 

yInit = [12 16 32 23 40 50]; % initial conditions (volumes)

%initialise ejection/filling coefficients
alpha_rv = 0;
alpha_lv = 1;
beta_rv = 0;
beta_lv = 1;

for k = 1:numel(tPool) - 1

    %ventricles (time varying elastance)
    
    P_es_lv = E_es_lv*(V_lv - V_0); % LV ESPVR

    P_ed_lv = A*(exp(B_lv*(V_lv - V_0)) - 1); % LV EDPVR

    P_es_rv = E_es_rv*(V_rv - V_0); % RV ESPVR

    P_ed_rv = A*(exp(B_rv*(V_rv - V_0)) - 1); %RV EDPVR

    switch k
        case 1
            e = 0.5*(1-cos(2*pi*t/T_es));
        case 2
            e = 0; 
        case 3
            e = 0;
    end 

    P_lv = (P_es_lv - P_ed_lv)*e + P_ed_lv; %left ventricular pressure

    P_rv = (P_es_rv - P_ed_rv)*e + P_ed_rv; %right ventricular pressure
        

The above code represents this section (except that I have implemented a simpler version of the piecewise e(t) function:

    %systemic and pulmonary circulation

    P_C_as = V_C_as/C_as; % systemic arteries
    P_C_vs = V_C_vs/C_vs; % systemic veins
    P_C_ap = V_C_ap/C_ap; % pulmonary arteries
    P_C_vp = V_C_vp/C_vp; % pulmonary veins

The above code represents this section:

    %blood volume

    V_T = V_C_as + V_C_vs + V_C_ap + V_C_vp + V_U + V_lv + V_rv;

    %differential equations

    ode1 = diff(V_lv) == ((P_C_vp - P_lv)/R_vp)*alpha_lv - ((P_lv - P_C_as)/R_cs)*beta_lv;

    ode2 = diff(V_rv) == ((P_C_vs - P_rv)/R_vp)*alpha_rv - ((P_rv - P_C_ap)/R_cp)*beta_rv;

    ode3 = diff(V_C_as) == ((P_lv - P_C_as)/R_cs)*beta_lv - ((P_C_as - P_C_vs)/R_as);

    ode4 = diff(V_C_vs) == ((P_C_as - P_C_vs)/R_as) - ((P_C_vs - P_rv)/R_vs)*alpha_rv;

    ode5 = diff(V_C_ap) == ((P_rv - P_C_ap)/R_cp)*beta_rv - ((P_C_ap - P_C_vp)/R_ap);

    ode6 = diff(V_C_vp) == ((P_C_ap - P_C_vp)/R_ap) - ((P_C_vp - P_lv)/R_vp)*alpha_lv;

    odes = [ode1; ode2; ode3; ode4; ode5; ode6];

    [VF, Sbs] = odeToVectorField(odes);
  
    %solve numerically

    M = matlabFunction(VF,'vars',{'t','Y'});

    [T, Y] = ode45(M, [tPool(k:k+1)], yInit);

    tAll  = cat(1, tAll, T);  % Append new solution to total solution
    yAll  = cat(1, yAll, Y);
    
    yInit = Y(end, :);  % Update the initial value
    

The above code represents this section:

    %set ejection/filling coefficients
end

plot(tAll, yAll)

As seen in equations A11 to A16, there are four coefficients that are set to either 1 or 0 depending on whether ventricular ejection or filling is occuring:

I know that I can't compare these pressure values directly as they are symbolic. I have tried using subs, but that only works for the pressures that don't contain time-varying terms. For example:

 test1 = subs(P_lv, V_lv, Y(end, 1))  
 test2 = subs(P_C_vp, V_C_vp, Y(end, 6))

test1 can be converted into a double, but test2 can't.

I implemented the piecewise e(t) function based on an answer from this post, but I can't figure out how to change the values of alpha and beta at each t value based on logical comparisons of symbolic variables.

Any help is appreciated.

Answer 1

[Star Strider]

It can, actually, using the vpa function:

test1 = vpa(subs(P_lv, V_lv, Y(end, 1)), 8)

test2 = vpa(subs(P_C_vp, V_C_vp, Y(end, 6)), 8)

however something is missing here, so I cannot run the code —

%create symbols

syms V_lv(t) V_rv(t) V_C_as(t) V_C_vs(t) V_C_ap(t) V_C_vp(t) e(t)

tAll  = [];
yAll  = [];

tPool = [0 2*T_es T]; % time periods in the cardiac cycle 

Unrecognized function or variable 'T_es'.

yInit = [12 16 32 23 40 50]; % initial conditions (volumes)

%initialise ejection/filling coefficients
alpha_rv = 0;
alpha_lv = 1;
beta_rv = 0;
beta_lv = 1;

for k = 1:numel(tPool) - 1

    %ventricles (time varying elastance)
    
    P_es_lv = E_es_lv*(V_lv - V_0); % LV ESPVR

    P_ed_lv = A*(exp(B_lv*(V_lv - V_0)) - 1); % LV EDPVR

    P_es_rv = E_es_rv*(V_rv - V_0); % RV ESPVR

    P_ed_rv = A*(exp(B_rv*(V_rv - V_0)) - 1); %RV EDPVR

    switch k
        case 1
            e = 0.5*(1-cos(2*pi*t/T_es));
        case 2
            e = 0; 
        case 3
            e = 0;
    end 

    P_lv = (P_es_lv - P_ed_lv)*e + P_ed_lv; %left ventricular pressure

    P_rv = (P_es_rv - P_ed_rv)*e + P_ed_rv; %right ventricular pressure
    
    
        %set ejection/filling coefficients
end

plot(tAll, yAll)

 test1 = vpa(subs(P_lv, V_lv, Y(end, 1)), 8)  
 test2 = vpa(subs(P_C_vp, V_C_vp, Y(end, 6)), 8)

.

[Hannah Goldblatt]

Thank you for the response. However, vpa still doesn't allow me to evaluate P_lv at time t. So the code:

 test1 = vpa(subs(P_lv, V_lv, Y(end, 1)), 8)
 test2 = vpa(subs(P_C_vp, V_C_vp, Y(end, 6)), 8)

    if test1 < test2
        alpha_lv = 1;
    else 
        alpha_lv = 0;
    end     

has the output:

[Star Strider]

If an expression contains symbolic variables, vpa will simply convert the numeric values to decimal fractions. The symbolic variables need to have numeric values assigned to them in order to produce a completely numeric result (and possibly be converted using the double function). I can’t run your code, so I can’t specifically test this with it.

Symbolic variables will not work with logical tests, unless the variables are assigned values and so produce a completely numeric result. Use the isAlways function to prove a symbolic result in a specific situation.

[Hannah Goldblatt]

I'll include the entirety of the code below. If I want to evaluate and compare symbolic variables at each point in time, should I set up a for loop and run ode45 for a duration of 1 within each run of the loop? The issue is that I need to know the pressures (and therefore the volumes) at each point in time in order to change the alpha and beta values.

It's likely that I'm misunderstanding the ode45 solver, but if it integrates the differential equations from t0 to tf in the given duration, is it not possible to change values at each point t within the [t0 tf] period? Should I rather create an anonymous function where the volumes and alpha and beta values are parameters?

The substitutions I did in test1 and test2 are not ideal, because they keep the V_lv(t) and V_C_vp(t) values constant over each interval in tPool, when these volumes should be updating at each point in time. In this case, I'm plotting one cardiac cycle, so every parameter should update every ms in order to properly track changes in volume and pressure throughout the cardiac system.

%Parameters from Burkhoff and Tyberg

%Heart

HR = 75; % heart rate (bpm)
HR_ms = HR/60000; % heart rate (bpms)
T = 1/HR_ms; % period of cardiac cycle (ms)
T_es = 175; % time to end systole (ms)
tau = 25; % time const of relaxation (ms)

V_0 = 0; % unstressed volume (ml) 
A = 0.35; % scaling factor for EDPVR (mmHg)
B_lv = 0.033; % exponent for LV EDPVR (1/ml)
B_rv = 0.023; % exponent for RV EDPVR (1/ml)
E_es_lv = 3; % LV end-systolic elastance (mmHg/ml)
E_es_rv = 0.7; % RV end-systolic elastance (mmHg/ml)

%Circulation

%Pulmonary
R_ap = 40; % arterial resistance (dyn.s/cm5)
R_cp = 27; % characteristic resistance (dyn.s/cm5)
R_vp = 20; % venous resistance (dyn.s/cm5)
R_Tp = 87; % total resistance (dyn.s/cm5)
C_ap = 13; % arterial capacitance (ml/mmHg)
C_vp = 8; % venous capacitance (ml/mmHg)

%Systemic
R_as = 1200; % arterial resistance (dyn.s/cm5)
R_cs = 40; % characteristic resistance (dyn.s/cm5)
R_vs = 20; % venous resistance (dyn.s/cm5)
R_Ts = 1260; % total resistance (dyn.s/cm5)
C_as = 1.32; % arterial capacitance (ml/mmHg)
C_vs = 70; % venous capacitance (ml/mmHg)

%Common 

V_T = 5500; % total blood volume (ml)
V_S = 750; % total stressed volume (ml)
V_U = 4750; % total unstressed volume (ml)

%create symbols

syms V_lv(t) V_rv(t) V_C_as(t) V_C_vs(t) V_C_ap(t) V_C_vp(t) e(t)

tAll  = [];
yAll  = [];

tPool = [0 2*T_es T]; % time periods in the cardiac cycle 

yInit = [12 16 32 23 40 50]; % initial conditions (volumes)

%initialise ejection/filling coefficients
alpha_rv = 0;
alpha_lv = 1;
beta_rv = 0;
beta_lv = 1;

for k = 1:numel(tPool) - 1

    %ventricles (time varying elastance)
    
    P_es_lv = E_es_lv*(V_lv - V_0); % LV ESPVR

    P_ed_lv = A*(exp(B_lv*(V_lv - V_0)) - 1); % LV EDPVR

    P_es_rv = E_es_rv*(V_rv - V_0); % RV ESPVR

    P_ed_rv = A*(exp(B_rv*(V_rv - V_0)) - 1); %RV EDPVR

    switch k
        case 1
            e = 0.5*(1-cos(2*pi*t/T_es));
        case 2
            e = 0; 
        case 3
            e = 0;
    end 

    P_lv = (P_es_lv - P_ed_lv)*e + P_ed_lv; %left ventricular pressure

    P_rv = (P_es_rv - P_ed_rv)*e + P_ed_rv; %right ventricular pressure

    %systemic and pulmonary circulation

    P_C_as = V_C_as/C_as; % systemic arteries
    P_C_vs = V_C_vs/C_vs; % systemic veins
    P_C_ap = V_C_ap/C_ap; % pulmonary arteries
    P_C_vp = V_C_vp/C_vp; % pulmonary veins

    %blood volume

    V_T = V_C_as + V_C_vs + V_C_ap + V_C_vp + V_U + V_lv + V_rv;

    %differential equations

    ode1 = diff(V_lv) == ((P_C_vp - P_lv)/R_vp)*alpha_lv - ((P_lv - P_C_as)/R_cs)*beta_lv;

    ode2 = diff(V_rv) == ((P_C_vs - P_rv)/R_vp)*alpha_rv - ((P_rv - P_C_ap)/R_cp)*beta_rv;

    ode3 = diff(V_C_as) == ((P_lv - P_C_as)/R_cs)*beta_lv - ((P_C_as - P_C_vs)/R_as);

    ode4 = diff(V_C_vs) == ((P_C_as - P_C_vs)/R_as) - ((P_C_vs - P_rv)/R_vs)*alpha_rv;

    ode5 = diff(V_C_ap) == ((P_rv - P_C_ap)/R_cp)*beta_rv - ((P_C_ap - P_C_vp)/R_ap);

    ode6 = diff(V_C_vp) == ((P_C_ap - P_C_vp)/R_ap) - ((P_C_vp - P_lv)/R_vp)*alpha_lv;

    odes = [ode1; ode2; ode3; ode4; ode5; ode6];

    [VF, Sbs] = odeToVectorField(odes);
  
    %solve numerically

    M = matlabFunction(VF,'vars',{'t','Y'});

    [T, Y] = ode45(M, [tPool(k:k+1)], yInit);

    tAll  = cat(1, tAll, T);  % Append new solution to total solution
    yAll  = cat(1, yAll, Y);
    
    yInit = Y(end, :);  % Update the initial value

    %set ejection/filling coefficients

    test1_1 = vpa(subs(P_lv), 8)
    
    test2_1 = vpa(subs(P_C_vp), 8)

    test1 = vpa(subs(P_lv, V_lv, Y(end, 1)), 8)

    test2 = vpa(subs(P_C_vp, V_C_vp, Y(end, 6)), 8)

    if test1 < test2
        alpha_lv = 1;
    else 
        alpha_lv = 0;
    end 

    if P_rv < P_C_vs
        alpha_rv = 1;
    else 
        alpha_rv = 0;
    end 

    if P_lv > P_C_as
        beta_lv = 1;
    else 
        beta_lv = 0;
    end 
    
    if P_rv > P_C_ap
        beta_rv = 1;
    else 
        beta_rv = 0;
    end 

end

test1_1(t) = 

test2_1(t) = 

test1(t) = 

test2(t) = 

5.9300428

Conversion to logical from symfun is not possible.

plot(tAll, yAll)

[Star Strider]

You need to evaluate ‘test1_1’, ‘test2_1’, ‘test1’ and ‘test2’ at one or more specific values of ‘t’ each time. The ‘t’ values should be the same values and sizes if you want to compare the functions using them. They can’t be compared symboilically.

When you do that, it may also be necessary for you to use the any or all functions to do the logical comparisons. (That you’re doing cardiovascular modeling got my attention!)

%Parameters from Burkhoff and Tyberg

%Heart

HR = 75; % heart rate (bpm)
HR_ms = HR/60000; % heart rate (bpms)
T = 1/HR_ms; % period of cardiac cycle (ms)
T_es = 175; % time to end systole (ms)
tau = 25; % time const of relaxation (ms)

V_0 = 0; % unstressed volume (ml) 
A = 0.35; % scaling factor for EDPVR (mmHg)
B_lv = 0.033; % exponent for LV EDPVR (1/ml)
B_rv = 0.023; % exponent for RV EDPVR (1/ml)
E_es_lv = 3; % LV end-systolic elastance (mmHg/ml)
E_es_rv = 0.7; % RV end-systolic elastance (mmHg/ml)

%Circulation

%Pulmonary
R_ap = 40; % arterial resistance (dyn.s/cm5)
R_cp = 27; % characteristic resistance (dyn.s/cm5)
R_vp = 20; % venous resistance (dyn.s/cm5)
R_Tp = 87; % total resistance (dyn.s/cm5)
C_ap = 13; % arterial capacitance (ml/mmHg)
C_vp = 8; % venous capacitance (ml/mmHg)

%Systemic
R_as = 1200; % arterial resistance (dyn.s/cm5)
R_cs = 40; % characteristic resistance (dyn.s/cm5)
R_vs = 20; % venous resistance (dyn.s/cm5)
R_Ts = 1260; % total resistance (dyn.s/cm5)
C_as = 1.32; % arterial capacitance (ml/mmHg)
C_vs = 70; % venous capacitance (ml/mmHg)

%Common 

V_T = 5500; % total blood volume (ml)
V_S = 750; % total stressed volume (ml)
V_U = 4750; % total unstressed volume (ml)

%create symbols

syms V_lv(t) V_rv(t) V_C_as(t) V_C_vs(t) V_C_ap(t) V_C_vp(t) e(t)

tAll  = [];
yAll  = [];

tPool = [0 2*T_es T]; % time periods in the cardiac cycle 

yInit = [12 16 32 23 40 50]; % initial conditions (volumes)

%initialise ejection/filling coefficients
alpha_rv = 0;
alpha_lv = 1;
beta_rv = 0;
beta_lv = 1;

figure

for k = 1:numel(tPool) - 1

    %ventricles (time varying elastance)
    
    P_es_lv = E_es_lv*(V_lv - V_0); % LV ESPVR

    P_ed_lv = A*(exp(B_lv*(V_lv - V_0)) - 1); % LV EDPVR

    P_es_rv = E_es_rv*(V_rv - V_0); % RV ESPVR

    P_ed_rv = A*(exp(B_rv*(V_rv - V_0)) - 1); %RV EDPVR

    switch k
        case 1
            e = 0.5*(1-cos(2*pi*t/T_es));
        case 2
            e = 0; 
        case 3
            e = 0;
    end 

    P_lv = (P_es_lv - P_ed_lv)*e + P_ed_lv; %left ventricular pressure

    P_rv = (P_es_rv - P_ed_rv)*e + P_ed_rv; %right ventricular pressure

    %systemic and pulmonary circulation

    P_C_as = V_C_as/C_as; % systemic arteries
    P_C_vs = V_C_vs/C_vs; % systemic veins
    P_C_ap = V_C_ap/C_ap; % pulmonary arteries
    P_C_vp = V_C_vp/C_vp; % pulmonary veins

    %blood volume

    V_T = V_C_as + V_C_vs + V_C_ap + V_C_vp + V_U + V_lv + V_rv;

    %differential equations

    ode1 = diff(V_lv) == ((P_C_vp - P_lv)/R_vp)*alpha_lv - ((P_lv - P_C_as)/R_cs)*beta_lv;

    ode2 = diff(V_rv) == ((P_C_vs - P_rv)/R_vp)*alpha_rv - ((P_rv - P_C_ap)/R_cp)*beta_rv;

    ode3 = diff(V_C_as) == ((P_lv - P_C_as)/R_cs)*beta_lv - ((P_C_as - P_C_vs)/R_as);

    ode4 = diff(V_C_vs) == ((P_C_as - P_C_vs)/R_as) - ((P_C_vs - P_rv)/R_vs)*alpha_rv;

    ode5 = diff(V_C_ap) == ((P_rv - P_C_ap)/R_cp)*beta_rv - ((P_C_ap - P_C_vp)/R_ap);

    ode6 = diff(V_C_vp) == ((P_C_ap - P_C_vp)/R_ap) - ((P_C_vp - P_lv)/R_vp)*alpha_lv;

    odes = [ode1; ode2; ode3; ode4; ode5; ode6];

    [VF, Sbs] = odeToVectorField(odes);
  
    %solve numerically

    M = matlabFunction(VF,'vars',{'t','Y'});

    [T, Y] = ode45(M, [tPool(k:k+1)], yInit);
    
    subplot(3,1,k)
    plot(T, Y)
    grid
    title("k = "+string(k))

    tAll  = cat(1, tAll, T);  % Append new solution to total solution
    yAll  = cat(1, yAll, Y);
    
    yInit = Y(end, :);  % Update the initial value

    %set ejection/filling coefficients

    test1_1 = vpa(subs(P_lv), 8)
    
    test2_1 = vpa(subs(P_C_vp), 8)

    test1 = vpa(subs(P_lv, V_lv, Y(end, 1)), 8)

    test2 = vpa(subs(P_C_vp, V_C_vp, Y(end, 6)), 8)

    if test1 < test2
        alpha_lv = 1;
    else 
        alpha_lv = 0;
    end 

    if P_rv < P_C_vs
        alpha_rv = 1;
    else 
        alpha_rv = 0;
    end 

    if P_lv > P_C_as
        beta_lv = 1;
    else 
        beta_lv = 0;
    end 
    
    if P_rv > P_C_ap
        beta_rv = 1;
    else 
        beta_rv = 0;
    end 

end

test1_1(t) = 

test2_1(t) = 

test1(t) = 

test2(t) = 

5.9300428

Conversion to logical from symfun is not possible.

plot(tAll, yAll)

.

[Hannah Goldblatt]

Thank you! I have changed the k array so that it represents increments of 8ms, and I substitute k for t in order to evaluate the volumes at particular points in time. It's not very elegant but it seems to be working. For interest's sake, the paper I'm using is "Why does pulmonary venous pressure rise after onset of LV dysfunction: a theoretical analysis" by Burkhoff and Tyberg.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

It would be nice to have a PDF of the paper.

.

[Hannah Goldblatt]

Sure thing, here is the pdf.

[Star Strider]

Thank you!

I’m interestied in physiological models. I wasn’t aware of this one.