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

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.

Please help in fixing the code.

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')

Answer 1

[William Rose]

@Timothy Jen Roxas,

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.

[William Rose]

@Timothy Jen Roxas,

I wrote a "sinusoid" function and placed it at the end of your script, just to show that it works. It runs without error. See attached. Change the sinusoid function to make it do what you want.

I changed the line that estimates diastolic pressure. You had

dia_p = min(pressure_si); % diastolic

I changed it to

dia_p = min(pressure_si(time_si>0.85*time_si(end))); % diastolic

Your line of code searches the entire pressure signal to find the minimum, so it finds the initial value of pressure, p=0, and calls it diastolic pressure. But that is not the final diastolic pressure, at least not when my sinusoid() is used. My line for dia_p searches the final 15% of the pressure signal (about one and a half beats). Therefore it finds the minimum pressure within the last beat and a half of the simulation. That is better.

[William Rose]

@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