# Why doesn't the code for cardiac cycle with varying resistance and capacitance run?

2 views (last 30 days)
Timothy Jen Roxas on 11 Nov 2021
Commented: William Rose on 11 Nov 2021
Whenever I attempt to run the code, this prompt appears:
"Function definitions in a script must appear at the end of the file.
Move all statements after the "sinusoid" function definition to before the first local
function definition"
I tried to add functions before the bracketed statements, but I can't figure out how.
Here is the code:
%% sinusoidal inflow
% clc;
% clear all
R = 1; % vascular resistance in mmHg-s/ml for healthy human
C = 1; % capacitance in ml/mmHg
tc = 0.8; % length of cardiac cycle (s)
[time_si,pressure_si] = sinusoid(R,C,tc);
%% Analyze
sys_p = max(pressure_si); % systolic
dia_p = min(pressure_si); % diastolic
map = 1/3*sys_p + 2/3*dia_p; % mean arterial pressure
pulse_pressure = sys_p - dia_p;
%% Change Resistance
clc
% increase resistance by 25%
R1 = 1.25;
C = 1;
[time_r1,pressure_r1] = sinusoid(R1,C,tc);
sys_p1 = max(pressure_r1)
dia_p1 = min(pressure_r1);
map1 = 1/3*sys_p1 + 2/3*dia_p1;
pulse_pressure_r1 = sys_p1 - dia_p1
% decrease resistance by 25%
R2 = 0.75;
C = 1;
[time_r2,pressure_r2] = sinusoid(R2,C,tc);
sys_p2 = max(pressure_r2)
dia_p2 = min(pressure_r2)
map2 = 1/3*sys_p2 + 2/3*dia_p2;
pulse_pressure_r2 = sys_p2 - dia_p2
lw = 3;
fs = 18;
figure
hold on
grid on
plot(time_si,pressure_si,'g','LineWidth',lw)
plot(time_r1,pressure_r1,'r','LineWidth',lw)
plot(time_r2,pressure_r2,'b','LineWidth',lw)
legend('Normal R','R increased 25%','R decreased 25%')
set(gca,'fontsize',fs)
ylim([0 150])
xlabel('Time [s]')
ylabel('Ventricular Pressure [mmHg]')
%% Change Capacitance
clc
% increase capacitance by 25%
R = 1;
C1 = 1.25;
[time_c1,pressure_c1] = sinusoid(R,C1,tc);
sys_p1 = max(pressure_c1)
dia_p1 = min(pressure_c1);
map1 = 1/3*sys_p1 + 2/3*dia_p1;
pulse_pressure_c1 = sys_p1 - dia_p1
% decrease capacitance by 25%
R = 1;
C2 = 0.75;
[time_c2,pressure_c2] = sinusoid(R,C2,tc);
sys_p2 = max(pressure_c2)
dia_p2 = min(pressure_c2)
map2 = 1/3*sys_p2 + 2/3*dia_p2;
pulse_pressure_c2 = sys_p2 - dia_p2
lw = 3;
fs = 18;
figure
hold on
grid on
plot(time_si,pressure_si,'g','LineWidth',lw)
plot(time_c1,pressure_c1,'r','LineWidth',lw)
plot(time_c2,pressure_c2,'b','LineWidth',lw)
legend('Normal C','C increased 25%','C decreased 25%')
set(gca,'fontsize',fs)
ylim([0 150])
xlabel('Time [s]')
ylabel('Ventricular Pressure [mmHg]')
%%
figure
subplot(121)
hold on
grid on
plot(time_si,pressure_si,'g','LineWidth',lw)
plot(time_r1,pressure_r1,'r','LineWidth',lw)
plot(time_r2,pressure_r2,'b','LineWidth',lw)
legend('Normal R','R increased 25%','R decreased 25%')
set(gca,'fontsize',fs)
ylim([0 150])
xlabel('Time [s]')
ylabel('Ventricular Pressure [mmHg]')
title('Effect of Resistance on Pressure')
subplot(122)
hold on
grid on
plot(time_si,pressure_si,'g','LineWidth',lw)
plot(time_c1,pressure_c1,'r','LineWidth',lw)
plot(time_c2,pressure_c2,'b','LineWidth',lw)
legend('Normal C','C increased 25%','C decreased 25%')
set(gca,'fontsize',fs)
ylim([0 150])
xlabel('Time [s]')
ylabel('Ventricular Pressure [mmHg]')
title('Effect of Capacitance on Pressure')

William Rose on 11 Nov 2021
I saved your program as script CVSimMatlabCentral.m. I ran it, with this result:
>> CVSimMatlabCentral
Unable to perform assignment because dot indexing is not supported for variables of this type.
Error in sinusoid (line 32)
mstruct.mapparallels = 0;
Error in CVSimMatlabCentral (line 11)
[time_si,pressure_si] = sinusoid(R,C,tc);
The problem here is that sinusoid() is a built-in Matlab function. It is a kind of map projection. It expects inputs different than those you have provided, hence the error.
I suspect you have written your own sinsoid() routine, but it is not in the code you posted. If you append your sinusoid() to your script, at the end, it will be used by your program, instead of the built-in sinusoid().
If you reply to this post, please attach your script as an .m file.
##### 2 CommentsShowHide 1 older comment
William Rose on 11 Nov 2021
@Timothy Jen Roxas, you estimate MAP by MAP=SysP/3+2*DiaP/3.
That estimate can be quite inaccurate, for a variety of reasons. In your simulation, that estimate is quite inaccurate. It significantly underestimates the true MAP. You can compute the mean pressure more accurately by averaging the pressure for one beat duration at the end of the simulation:
mapA = mean(pressure_si(time_si>=time_si(end)-tc)); %accurate MAP
Try it. In the results below (generated with attached script), notice the differece between MAP and the accurate MAP.
>> CVSimMatlabCentral
Normal__________: SysP=127.5, DiaP= 74.8, MAP= 92.4, MAPa=100.0, PP=52.7
Resistance +25%: SysP=152.4, DiaP= 99.5, MAP=117.1, MAPa=125.0, PP=52.9
Resistance -25%: SysP=102.7, DiaP= 50.4, MAP= 67.8, MAPa= 75.0, PP=52.3
Capacitance +25%: SysP=121.9, DiaP= 79.6, MAP= 93.7, MAPa=100.0, PP=42.4
Capacitance -25%: SysP=136.9, DiaP= 67.2, MAP= 90.4, MAPa=100.0, PP=69.7

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!