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

Question

Timothy Jen Roxas 2021-11-11

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1583894-why-doesn-t-the-code-for-cardiac-cycle-with-varying-resistance-and-capacitance-run

评论： William Rose 2021-11-11

采纳的回答： William Rose

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2021-11-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1583894-why-doesn-t-the-code-for-cardiac-cycle-with-varying-resistance-and-capacitance-run#answer_828969

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

2 个评论
显示无隐藏无

William Rose 2021-11-11

在 MATLAB Online 中打开

CVSimMatlabCentral.m

@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 2021-11-11

在 MATLAB Online 中打开

CVSimMatlabCentral.m

@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