% The function simulate.m implements various computational models of the
% human cardiovascular system. The models include the following three
% major components:
%
% 1) Lumped parameter, pulsatile heart and circulation
% (various preparations)
% 2) Short-term regulatory system
% a) Arterial baroreflex system
% b) Cardiopulmonary baroreflex system
% c) Direct neural coupling mechanism between Qlu and F
% 3) Resting physiologic perturbations
% a) respiratory activity
% b) eExogenous disturbances to
% i) Ra
% ii) F
% iii) Qvo
% iv) Pasp
%
% Function arguments:
% flag - 8x1 vector indicating the status of the model
% (overrides parameter values assigned in th)
% - [preparation breathing dncm baro dra dpasp df dqvo]
%
% - preparation - status of pulsatile heart and circulation
% - 0 - intact circulation (nominal)
% - 1 - heart-lung unit designed for measuring cardiac function curves
% - 2 - systemic circulation designed for measuring venous return curves
% - 3 - left ventricle
% (automated step changes in Ppv and Pa NOT implemented)
% - 4 - intact circulation for measurement of single ventricular
% contraction Pa response
% - 5 - intact circulation with only linear elements
% - 6 - intact circulation with third-order systemic arteries
% - 7 - intact circulation with nonlinear systemic arterial compliance
% - 8 - intact circulation with third-order systemic arteries
% and nonlinear large, elastic arterial compliance
% - 9 - intact circulation with an arterial pressure reservoir
% - 10 - intact circulation with contracting atria
% - breathing - status of breathing (Qlu)
% - 0 - no breathing (Qlu constant)
% - 1 - fixed-rate breathing (Qlu periodic)
% - 2 - random-interval breathing (Qlu random; 5 s mean period)
% instantaneous alveolar ventilation rate constant over time
% - 3 - random-interval breathing (Qlu random; 5 s mean period)
% tidal volume fixed; instantaneous alveolar ventilation rate
% not constant
% - 4 - step increase in breathing (Qlu step)
% - 5 - impulse increase in breathing (Qlu impulse)
% - dncm - status of direct neural coupling mechanism between Qlu and F
% - 0 - dncm off
% - 1 - dncm on (not applicable for step increase in Qlu)
% - baro - status of baroreflexes
% - 0 - arterial and cardiopulmonary baroreflexes off
% - 1 - arterial baroreflex on only
% - 2 - cardiopulmonary baroreflex on only
% - 3 - arterial and cardiopulmonary baroreflexes on
% - dra - status of bandlimited disturbance to Ra
% - 0 - dra off
% - 1-3 - dra on
% - 1 - 1/f fluctuations
% - 2 - bandlimited, white fluctuations
% - 3 - sinusoidal fluctuations
% - dpasp - status of sinusoidal disturbance to Pasp
% - 0 - dpasp off
% - 1 - dpasp on
% - df - status of 1/f disturbance to F
% - 0 - df off
% - 1 - df on
% - dqvo - status of bandlimited disturbance to Qvo
% - 0 - dqvo off
% - 1 - dqvo on
% th - 113x1 vector containing the parameter values characterizing
% the simulation
% outputfile - prefix of the output files in MIT format generated by the model
% parameterfile - name of a file which contains the parameter values
% characterizing the model and its execution
% signals - the signals to be viewed by WAVE as the data are generated
%
% NOTE: If this program is executed in the MATLAB environment, it will
% not be possible to implement on-line viewing and parameter updating
% Thus, the last three arguments should be omitted.
%
% Function outputs:
% P - matrix representing simulated pressures [mmHg]
% = [Pl Pe Pv Pr Ppa Ppv Pth Palv Pra Pm] for preparation == 6,8
% = [Pl Pa Pv Pr Ppa Ppv Pth Palv Pra Pla] for preparation == 10
% = [Pl Pa Pv Pr Ppa Ppv Pth Palv Pra] otherwise
% Q - matrix representing simulated volumes [ml]
% = [Ql Qe Qv Qr Qpa Qpv Qlu Qm] for preparation == 6,8
% = [Ql Qa Qv Qr Qpa Qpv Qlu Qra Qla] for preparation == 10
% = [Ql Qa Qv Qr Qpa Qpv Qlu] otherwise
% q - matrix representing simulated flow rates [ml/s]
% = [qpv ql qm qv qr qpa qe] for preparation == 6,8
% = [qla ql qa qra qr qpa qv qpv] for preparation == 10
% = [qpv ql qa qv qr qpa] otherwise
% ap - matrix representing simulated adjustable parameters
% = [Cls Crs Qvo Ra F Tsys]
% ve - matrix representing the ventricular variable elastances
% [mmHg/ml]
% = [El Er Ela Era] for preparation == 10
% = [El Er] otherwise
% qrs - vector representing times of ventricular contractions [s]
% t - vector representing the associated time [s]
% num - matrix representing either cardiac function or venous
% return curve numerics
% = [mra mpa mco] for preparation == 1
% = [mra mco] for preparation == 2
% = [0, 0] otherwise
%
% Function usage (in MATLAB):
%
% [P,Q,q,ap,ve,qrs,t,num] = simulate(flag,th);
%
function [P,Q,q,ap,ve,qrs,t,num] = simulate(flag,th,outputfile,parameterfile,signals)
% Declaring C functions.
%#function wave_remote read_key write_param
% Checking number of arguments. If there are only two arguments,
% this indicates MATLAB execution and thus on-line viewing is
% turned off. (If it were left on, MATLAB would produce an error.)
if (nargin == 2)
signals = '-1';
end
% Checking dimensionality of the argument flag -- status parameters.
% If flag is of incorrect dimension, the program will exit in error.
if (length(flag) ~= 8)
error('flag is not of correct dimension');
end
% Assigning status parameters to new variable names.
preparation = flag(1);
breathing = flag(2);
dncm = flag(3);
baro = flag(4);
dra = flag(5);
dpasp = flag(6);
df = flag(7);
dqvo = flag(8);
% Assigning numerical integration variables.
% Number of integration steps.
tics = th(106);
% Integration step size (s).
deltaT = 1/th(113);
% Assigning integrate and fire variables.
% Fraction of current cardiac cycle that integration step is at.
ipfm = 0;
% Current number of ventricular contractions.
qrs_index = 1;
% Time in current cardiac cycle that integration is at.
sumdeltaT = 0;
% Short-term regulatory systems and resting physiologic
% perturbation variables.
% Note there are restrictions on altering Sgran and Sratio.
% 1) Sgran*Sratio <= 1.0.
% 2) (1.5/Sgran) must be an integer.
% 3) ((2/Sgran)-((Sratio*Sgran/2)/Sgran)) must be an integer.
% 4) (0.5-(Sratio*Sgran/2))/Sgran must be an integer.
% Granularity (s).
Sgran = 0.0625;
% Ratio of length of averaging window to Sgran
% (must be an integer).
Sratio = 4;
% Duration (samples).
Slength = 40/Sgran;
% Generation of a vector of arrival times for each respiratory cycle.
% at(1) = 0 & length(at) ~= 1: Qlu periodic
if (breathing == 1)
at = zeros(ceil((deltaT*tics)+1)+1,1);
i = 2;
x = (at(i) <= deltaT*tics+1);
mbintscalar(x);
while (x)
i = i+1;
at(i) = at(i-1)+th(42);
x = (at(i) <= deltaT*tics+1);
mbintscalar(x);
end
at = at(1:i);
% at(1) = 1: Qlu random with constant alveolar ventilation rate (2)
% at(1) = 2: Qlu random with fixed tidal volume (3)
elseif (breathing == 2 | breathing == 3)
at = zeros(ceil((deltaT*tics)+1)+1,1);
if (breathing == 2)
at(1) = 1;
else
at(1) = 2;
end
i = 2;
at(i) = rand_int_breath;
x = (at(i) <= deltaT*tics+1);
mbintscalar(x);
while (x)
i = i+1;
at(i) = at(i-1) + rand_int_breath;
x = (at(i) <= deltaT*tics+1);
mbintscalar(x);
end
at = at(1:i);
% at(1) = 0 & length(at) == 1: Qlu constant
elseif (breathing == 0)
at = 0;
% at(1) = 3 (length(at) == 2): step Qlu (cannot implement dncm)
elseif (breathing == 4)
at = [3; th(36)];
% at(1) = 4 (length(at) == 3): impulse Qlu
elseif (breathing == 5)
at = [4; th(36); (((deltaT*tics)+1)+1)];
% Exiting program with error if breathing is not valid.
else
error('breathing is not a valid option');
end
% Pre-calculating Qlu, Pth, dPth, and Ppac (Palv) over the entire
% integration period from the generated arrival time vector.
thbreathe = resp_act(th,at,deltaT/2,0,(tics-1)*deltaT);
% Adding a step change to Qlu, Pth, dPth, and Ppac, if mandated via Qfrs.
if ((breathing == 0 | breathing == 1 | breathing == 2) & th(36) ~= 0)
th(103) = 0;
thbreathes = resp_act(th,[3; th(36)],deltaT/2,0,(tics-1)*deltaT);
dthbreathe = [-th(98)*ones(1,length(thbreathes)); -th(23)*ones(1,length(thbreathes)); zeros(2,length(thbreathes))];
thbreathe = thbreathe+thbreathes+dthbreathe;
end
% Creating the filter characterizing direct neural coupling mechanism
% betweeen Qlu and F at a sampling period of Sgran. Decimating the
% above Qlu to a sampling period of Sgran in order to implement this filter.
% (not applicable to step input of Qlu)
if (dncm == 1)
hilvhr = dncm_filt(th,[Sgran; Slength; Sratio]);
ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran);
% Taking care of filtering edge effects.
ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];
rsastep = 1;
% Exiting program with error if dncm is not valid.
elseif (dncm ~= 0)
error('dncm is not a valid option');
end
% Autoregulation of local vascular beds -- exogenous disturbance to Ra.
% Generating sinusoidal Ra disturbance at a sampling period of Sgran
% for entire integration period.
if (dra == 3)
amp = th(61);
freq = th(62);
n = floor(((deltaT*tics)+1)/Sgran)+(Sratio-1);
i = 1:n;
tr = (i-1)*Sgran;
Ralvrp = (amp/(Sgran*2*pi*freq))*(cos(2*pi*freq*tr)-cos(2*pi*freq*(tr+Sgran)));
Ralvr = zeros(1,length(Ralvrp)-Sratio+1);
for kt = Sratio:length(Ralvrp)
Ralvr(kt-(Sratio-1)) = sum(Ralvrp(kt-(Sratio-1):kt))/Sratio;
end
Ralvr = Ralvr+((0.01*amp)/sqrt(Sgran))*randn(1,length(Ralvr));
lvrstep = 2;
% Generating bandlimited 1/f Ra disturbance at a sampling period of Sgran
% for entire integration period.
elseif (dra == 1)
hra = bl_filt(th(44),Sgran);
n = floor(((deltaT*tics)+1)/Sgran)+length(hra);
whitenoise = sqrt(1/Sgran)*th(47)*randn(1,n);
dmin = -4;
dmax = 0;
honeoverfra = oneoverf_filt(th,Sgran,dmin,dmax);
Ralvrw = conv(honeoverfra,whitenoise);
Ralvrw = Ralvrw(length(Ralvrw)-length(n)+1:length(Ralvrw));
% Filtering 1/f disturbance to Ra.
Ralvr = conv(hra,Ralvrw);
Ralvr = Ralvr(length(hra)+1:length(Ralvrw));
lvrstep = 2;
% Generating white noise at a sampling period of Sgran for bandlimited
% disturbance to Ra.
elseif (dra == 2)
wra = sqrt(1/Sgran)*th(47)*randn(length(-60:Sgran:60),1);
% Exiting program with error if dra is not valid.
elseif (dra ~= 0)
error('dra is not a valid option');
end
% Generating sinusoidal Pasp disturbance at a sampling period of Sgran
% for entire integration period.
if (dpasp == 1);
amp = th(63);
freq = th(64);
n = floor(((deltaT*tics)+1)/Sgran)+(Sratio-1)+Slength;
i = 1:n;
tr = (i-1)*Sgran;
tempPsp = (amp/(Sgran*2*pi*freq))*(cos(2*pi*freq*tr)-cos(2*pi*freq*(tr+Sgran)));
deltaPsp = zeros(length(tempPsp)-Sratio+1,1);
for kt = Sratio:length(tempPsp)
deltaPsp(kt-(Sratio-1)) = sum(tempPsp(kt-(Sratio-1):kt))/Sratio;
end
deltaPsp = deltaPsp+((0.1*amp)/sqrt(Sgran))*randn(length(deltaPsp),1);
ceostep = 2;
elseif (dpasp == 0)
deltaPsp = 0;
% Exiting program with error if dpasp is not valid.
else
error('dpasp is not valid');
end
% Creating 1/f filter and white noise at a sampling period of Sgran
% for the exogenous disturbance to F.
if (df == 1)
dmin = -4; % Desired minimum frequency decade with 1/f character
dmax = 0; % Desired maximum frequency decade with 1/f character
hans = ans_filt(th,[Sgran; Slength; Sratio]);
honeoverf = oneoverf_filt(th,Sgran,dmin,dmax);
hf = conv(honeoverf,hans);
wf = sqrt(1/Sgran)*th(48)*randn(length(hf),1);
% Exiting program with error if df is not valid.
elseif (df ~= 0)
error('df is not a valid option');
end
% Generating bandlimited, white disturbance to Qvo at a sampling period of Sgran
% for entire integration period.
if (dqvo == 1)
hvvo = bl_filt(th(65),Sgran);
n = floor(((deltaT*tics)+1)/Sgran)+length(hvvo);
Vvow = sqrt(1/Sgran)*th(66)*randn(1,n);
Vvo = conv(hvvo,Vvow);
Vvo = Vvo(length(hvvo)+1:length(Vvow));
nvvostep = 2;
% Exiting program with error if dqvo is not valid.
elseif (dqvo ~= 0)
error('dqvo is not a valid option');
end
% Pre-allocating memory for function output waveforms which is
% necessary for dramatic improvement in execution speed.
P = zeros(11,tics);
Q = zeros(9,tics);
q = zeros(8,tics);
t = zeros(1,tics);
ap = zeros(9,tics);
ve = zeros(4,tics);
% Pre-allocating memory and initializing variables for
% cardiac function/venous return curve generation.
if (preparation == 1)
% Pre-allocating memory.
Crd = ((th(27)-th(12))/th(32))*th(6);
final = round((th(27)-th(12)-20)/Crd)+3;
rarange = [final:-th(88):0, th(23)+1];
lengthra = length(rarange);
if (lengthra == 2);
lengthra = 1;
rarange = final;
end
lengthsa = length(30:th(90):th(31)-10);
lengtht = lengthra*lengthsa;
num = zeros(3,lengtht);
% Initializing variables for family of cardiac function curves.
if (lengthra > 1 & lengthsa > 1)
countn = 0;
countj = 0;
th(30) = final;
th(29) = 30;
tmp = 0;
% Initializing variables for a single cardiac function curve.
elseif (lengthra > 1 & lengthsa == 1)
countn = 0;
countj = 0;
th(30) = final;
tmp = 0;
% Initializing variables for curve of average ql versus Pa.
elseif (lengthra == 1 & lengthsa > 1)
final = th(30);
rarange = final;
th(29) = 30;
tmp = 0;
end
% Initializing variables for averaging Pra and ql over each cardiac cycle.
summrap = 0;
summco = 0;
% Initializing variable representing number of beats until voltage
% source(s) is varied.
fb_int = 10;
elseif (preparation == 2)
% Pre-allocating memory.
Crs = ((th(27)-th(12))/th(32))*th(5);
lengthra = length([Crs:th(89):60]);
num = zeros(2,lengthra);
if (length(num(1,:)) > 1)
countn = 0;
Crd = Crs;
th(6) = (th(32)/(th(27)-th(12)))*Crd;
end
% Initializing variables for integrating Pra and qv over each cardiac cycle.
summrap = 0;
summco = 0;
% Initializing variable representing number of beats until Crd
% is varied.
fb_int = 10;
else
num = zeros(2,1);
end
% Resetting ventricular compliance values for the intact circulatory
% preparation with only linear elements.
if (preparation == 5)
th(1) = th(1)*((th(26)-th(9))/th(31));
th(2) = th(2)*((th(26)-th(9))/th(31));
th(5) = th(5)*((th(27)-th(12))/th(32));
th(6) = th(6)*((th(27)-th(12))/th(32));
th(31) = th(26)-th(9);
th(32) = th(27)-th(12);
end
% Initializing the function output waveforms.
% Adjustable parameters.
ap(1:6,1) = [th(1)*((th(26)-th(9))/th(31)) th(5)*((th(27)-th(12))/th(32)) th(11) th(16) th(22) .3*sqrt(th(24))]';
% Pulsatile heart and circulation preparation.
if (preparation == 0 | preparation == 4 | preparation == 5 | preparation == 7)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = intact_init_cond(th);
elseif (preparation == 1)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = hlu_init_cond(th);
elseif (preparation == 2)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = sc_init_cond(th);
elseif (preparation == 3)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = lv_init_cond(th);
elseif (preparation == 5)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = init_cond(th);
elseif (preparation == 6)
[Ptemp,Q(1:7,1),q(1:7,1),ve(1:2,1)] = third_init_cond(th);
P(1:6,1) = Ptemp(1:6);
P(10:11,1) = Ptemp(7:8);
Q(8,1) = Q(7,1);
elseif (preparation == 7)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = nac_init_cond(th);
elseif (preparation == 8)
[Ptemp,Q(1:7,1),q(1:7,1),ve(1:2,1)] = third_nac_init_cond(th);
P(1:6,1) = Ptemp(1:6);
P(10:11,1) = Ptemp(7:8);
Q(8,1) = Q(7,1);
elseif (preparation == 9)
[P(1:6,1),Q(1:6,1),q(1:6,1),ve(1:2,1)] = apr_init_cond(th);
elseif (preparation == 10)
[tempp,tempv,q(1:8,1),ve(:,1)] = a_init_cond(th);
P(1:6,1) = tempp(1:6);
P(9,1) = tempp(7);
P(10:11,1) = tempp(7:8);
Q(1:6,1) = tempv(1:6);
Q(8:9,1) = tempv(7:8);
% Exiting program with error if preparation is not valid,
else
error('preparation is not a valid option');
end
if (preparation ~= 10)
P(9,1) = P(4,1);
end
% Respiratory-related waveforms.
P(7,1) = thbreathe(2,1);
P(8,1) = thbreathe(4,1);
Q(7,1) = thbreathe(1,1);
% Assigning setpoint variables
% Setpoint pressure for arterial baroreflex (mmHg).
Pasp = th(40);
% Setpoint pressure for cardiopulmonary baroreflex (mmHg).
Pratrsp = th(41);
% F setpoint value (Hz).
Fsp = th(22);
% Ra setpoint value (mmHg-s/ml).
Rasp = th(16);
% Cls, Crs setpoint value (ml/mmHg).
Clssp = th(1);
Crssp = th(5);
% Qvo setpoint value (ml).
Qvosp = th(11);
% 7x1 vector of setpoint values.
thsp = [Pasp; Fsp; Rasp; Clssp; Crssp; Qvosp; Pratrsp];
% Initializing baroreflex-related variables.
if (baro > 0 & baro < 4)
% Initializing input pressures. (The lowest
% index indicates the most recent pressure.)
bP = zeros(Slength,1);
bPratr = zeros(Slength,1);
for i = 1:Slength;
bP(i) = P(2,1);
bPratr(i) = P(9,1)-P(7,1);
end
bPtemp = zeros(Sratio,1);
bPratrtemp = zeros(Sratio,1);
for i = 1:Sratio
bPtemp(i) = P(2,1);
bPratrtemp(i) = P(9,1)-P(7,1);
end
% Initializing variables for integrating Pa, Pratr,
% and t over each Sgran interval.
sumP = 0;
sumPratr = 0;
sumt = 0;
% Exiting program with error if baro is not valid,
elseif (baro ~= 0)
error('baro is not a valid option');
end
% Initializing parameter vector and other variables in order to
% implement linear interpolation of the controllable parameters.
thn = th;
thp = th;
tn = 0;
tp = 0;
% Initializing parameters for writing data and annotations to file
% in MIT format. Calling the WAVE window.
if (strcmp(signals,'-1') == 0)
% Initializing waveform file variables.
kstart = 1;
qrs_index_start = 1;
fid = fopen([outputfile '.dat'],'wa');
% Initializing annotation file.
fid2 = fopen([outputfile '.qrs'],'wa');
fwrite(fid2,1,'bit10');
fwrite(fid2,1,'ubit6');
fwrite(fid2,-1,'bit10');
fwrite(fid2,62,'ubit6');
Nbytes = length(parameterfile)+2;
fwrite(fid2,Nbytes,'bit10');
fwrite(fid2,63,'ubit6');
fwrite(fid2,[parameterfile '.' num2str(0)],[num2str(Nbytes) '*uchar']);
% Assigning annotation flag.
if (th(105) == 0)
annotator = '-1';
else
annotator = 'qrs';
end
% Initializing parameter update counter.
count_pc = 0;
% Calling the WAVE window.
wave_remote( outputfile, annotator, 0, signals );
% Initializing variables which permit annotation flag to be
% updated on-line.
wrflag = 1;
outputfile2 = ['./' outputfile];
end
% Generating the function output waveforms for each integration step.
for k = 2:tics
% Generating the pressure waveforms and associated time for
% the current integration step -- fourth-order Runge-Kutta
% integration.
Ptempn = rk4([P(1:6,k-1); P(10:11,k-1)],[Q(1,k-1);Q(4,k-1)],thbreathe(:,2*(k-1)+1:2*(k-1)+3),th,ipfm,sumdeltaT,deltaT,preparation);
P(1:6,k) = Ptempn(1:6);
if (preparation == 10)
P(9,k) = Ptempn(7);
end
P(10,k) = Ptempn(7);
P(11,k) = Ptempn(8);
% Updating time.
t(k) = t(k-1)+deltaT;
% Computing fraction of cardiac cycle that was surpassed
% during the integration step.
frac = deltaT*th(22);
% Obtaining Pth, Ppac (Palv), and Qlu at the current integration step
P(7,k) = thbreathe(2,2*(k-1)+1);
P(8,k) = thbreathe(4,2*(k-1)+1);
Q(7,k) = thbreathe(1,2*(k-1)+1);
% Integrating Pa data using the Backward Euler method and
% denoting the time into the integration.
sumt = sumt+deltaT;
sumP = sumP+P(2,k)*deltaT;
% Computing Pra and integrating Pratr using the Backward
% Euler method.
if (preparation ~= 10)
x = (P(4,k) >= P(3,k));
mbintscalar(x);
if (x)
P(9,k) = P(3,k);
else
P(9,k) = P(4,k);
end
end
sumPratr = sumPratr+(P(9,k)-P(7,k))*deltaT;
% Implementing the short-term regulatory system and/or resting
% physiologic perturbations each Sgran s.
x = (sumt > Sgran);
mbintscalar(x);
if (x)
% Setting final th values of previous Sgran to initial th
% value of current Sgran.
thp = thn;
tp = tn;
% Initializing thc -- a vector which contains the mandated
% adjustments to controllable parameters of sampling period
% Sgran.
thc = zeros(5,1);
% Implementing the arterial baroreflex.
if (baro == 1 | baro == 3)
% Calculating the integral of Pa precisely
% over the Sgran interval.
sumP = sumP-P(2,k)*deltaT;
sumP = sumP+P(2,k)*(deltaT-(sumt-Sgran));
% Placing the averaged Pa over the Sgran*Sratio
% interval in the most recent slot of the arterial
% baroreflex input.
bPtemp = [(sumP/Sgran); bPtemp(1:Sratio-1)];
avebP = sum(bPtemp)/Sratio;
bP = [avebP; bP(1:Slength-1)];
% Calculating the new baroreflex adjusted parameters
% for the next Sgran interval.
% For constant Pasp.
if (length(deltaPsp) == 1)
thc = thc+abreflex([Sgran; Slength; Sratio],bP,th,thsp,deltaPsp);
% For sinusoidally varying Pasp.
else
thc = thc+abreflex([Sgran; Slength; Sratio],bP,th,thsp,deltaPsp(dpasp:dpasp+length(bP)-1));
dpasp=dpasp+1;
end
end
% Implementing the cardiopulmonary baroreflex.
if (baro == 2 | baro == 3)
% Calculating the integral of Pratr precisely
% over the Sgran interval.
sumPratr = sumPratr-(P(9,k)-P(7,k))*deltaT;
sumPratr = sumPratr+(P(9,k)-P(7,k))*(deltaT-(sumt-Sgran));
% Placing the averaged Pratr over the Sgran*Sratio interval
% in the most recent slot of the cardiopulmonary baroreflex
% input.
bPratrtemp = [(sumPratr/Sgran); bPratrtemp(1:Sratio-1)];
avebPratr = sum(bPratrtemp)/Sratio;
bPratr = [avebPratr; bPratr(1:Slength-1)];
% Calculating the new baroreflex adusted parameters
% for the next Sgran interval.
thc = thc+cpreflex([Sgran; Slength; Sratio],bPratr,th,thsp);
end
% Implementing the direct neural coupling mechanism between
% Qlu and F.
if (dncm == 1)
% Creating the dncm filter.
hilvhr = dncm_filt(th,[Sgran; Slength; Sratio]);
% Creating the input Qlu vector to be convolved with the filter.
tempilv = ilv(rsastep:rsastep+length(hilvhr)-1);
tempilv = tempilv-th(98);
% Compute the mandated change in F from the convolution.
thc(1) = thc(1)+sum(tempilv.*hilvhr(length(hilvhr):-1:1));
rsastep = rsastep+1;
end
% Implementing the exogenous disturbance to Ra.
if (dra ~= 0)
% Bandlimited, white disturbance
if (dra == 2)
% Creating the bandlimited filter.
hra = bl_filt(th(44),Sgran);
% Creating white noise input to be convolved with the filter.
wra = [sqrt(1/Sgran)*th(47)*randn(1); wra(1:length(wra)-1)];
% Computing the mandated change to Ra.
thc(2) = thc(2)+sum(wra.*hra');
% Other Ra disturbance.
else
% Assigning mandated change to Ra which was pre-computed.
thc(2) = thc(2)+Ralvr(lvrstep);
lvrstep = lvrstep+1;
end
end
% Implementing the 1/f disturbance to F.
if (df == 1)
% Creating a filter which is a cascade combination of
% an autonomic filter and 1/f filter.
hans = ans_filt(th,[Sgran; Slength; Sratio]);
hf = conv(honeoverf,hans);
% Creating white noise input to be convolved with the filter.
wf = [sqrt(1/Sgran)*th(48)*randn(1); wf(1:length(wf)-1)];
% Computing mandated change to F.
thc(1) = thc(1)+sum(wf.*hf);
end
% Implementing disturbance to Qvo.
if (dqvo == 1)
% Assigning mandated change to Qvo which was pre-computed.
thc(5) = thc(5)+Vvo(nvvostep);
nvvostep = nvvostep+1;
end
% Calculating the values of the adjustable parameters at the
% end of the next Sgran interval and assigning these values
% to the thn vector.
thn = th;
thn(1) = thsp(4)+thc(3);
thn(5) = thsp(5)+thc(4);
thn(11) = thsp(6)+thc(5);
thn(16) = thsp(3)+thc(2);
thn(22) = thsp(2)+thc(1);
% Denoting the time at the start of the next Sgran interval.
tn = t(k)-(sumt-Sgran);
% Resetting the initial values for the integration
% of Pa, Pratr, and t over the next Sgran interval.
sumP = P(2,k)*(sumt-Sgran);
sumPratr = (P(9,k)-P(7,k))*(sumt-Sgran);
sumt = sumt-Sgran;
end
% Calculating the new values of the adjustable parameters
% at the current time step via linear interpolation and assigning
% the values to the vector tha. (Note that the mandated changes
% to Cls and Crs are implemented at the start of each beat.)
tha = th;
tha(11) = ((thn(11)-thp(11))/Sgran)*(t(k)-(tp+Sgran))+thp(11);
tha(16) = ((thn(16)-thp(16))/Sgran)*(t(k)-(tp+Sgran))+thp(16);
tha(22) = ((thn(22)-thp(22))/Sgran)*(t(k)-(tp+Sgran))+thp(22);
% Adjusting Pv in response to any mandated change in Qvo -- for
% constant volume considerations.
P(3,k) = P(3,k)+(th(11)-tha(11))/th(4);
% Setting the th vector to its current values.
th = tha;
% Computing Pra and qpv or qv over each cardiac cycle in order
% to establish cardiac function/venous return curves.
% (Note that ql is not computed, because it is relatively
% inaccurate for large integration steps.)
if (preparation == 1)
summrap = summrap+0.5*deltaT*(P(9,k-1)+P(9,k));
x = (P(6,k) > P(1,k));
mbintscalar(x);
if (x)
qpv = (P(6,k)-P(1,k))/th(20);
else
qpv = 0;
end
summco = summco+0.5*deltaT*(q(1,k-1)+qpv);
elseif (preparation == 2)
summrap = summrap+0.5*deltaT*(P(9,k-1)+P(9,k));
x = (P(3,k) > P(4,k) & P(4,k) >= th(91));
y = (P(3,k) > th(91) & th(91) > P(4,k));
mbintscalar(x);
mbintscalar(y);
if (x)
qv = (P(3,k)-P(4,k))/th(17);
elseif (y)
qv = (P(3,k)-th(91))/th(17);
else
qv = 0;
end
summco = summco+0.5*deltaT*(q(4,k-1)+qv);
end
% Recalculating the fraction advance in the cardiac cycle to account
% for the mandated change to F.
frac = ((deltaT-(sumt-Sgran))*th(22)) + ((sumt-Sgran)*tha(22));
% Updating time in and fraction of current cardiac cycle.
% These values are associated with the newly calculated pressure
% and time data. The integration of the ipfm parameter is achieved
% with the Backward Euler method.
sumdeltaT = sumdeltaT + deltaT;
% Stopping the initiation of ventricular contractions after a
% desired number of beats.
if (preparation == 4 & qrs_index > th(43))
sumdelta = 0;
frac = 0;
end
ipfm = ipfm + frac;
% Completing a cardiac cycle when ipfm is greater than or equal to one.
ipfmcond = (ipfm >= 1);
mbintscalar(ipfmcond)
if (ipfmcond)
% Resetting the ipfm model once the current cardiac contraction
% is complete.
sumdeltaT = (ipfm - 1)/th(22);
ipfm = ipfm - 1;
qrs_index = qrs_index+1; % Number of beats.
ap(7,qrs_index) = t(k)-sumdeltaT; % Time of each beat.
ap(8,qrs_index) = k; % Sample number of each beat.
% Setting next systolic interval based on previous
% cardiac cycle length.
th(24) = t(k)-sumdeltaT-ap(7,qrs_index-1);
% Implementing mandated changes to Cls and Crs.
th(1) = thn(1);
th(5) = thn(5);
% Varying parameters and computing averages for generation
% of cardiac function/venous return curves
if (preparation == 1)
% Calculating the integral of Pra and qpv
% precisely over the previous cardiac cycle.
summrap = summrap-0.5*deltaT*(P(9,k-1)+P(9,k));
summrap = summrap+0.5*(deltaT-sumdeltaT)*(P(9,k-1)+P(9,k));
summco = summco-0.5*deltaT*(q(1,k-1)+qpv);
summco = summco+0.5*(deltaT-sumdeltaT)*(q(1,k-1)+qpv);
% Computing the average Pra and qpv over the
% previous cardiac cycle. .
mrap = summrap/(t(k)-sumdeltaT-ap(7,qrs_index-1));
mco = summco/(t(k)-sumdeltaT-ap(7,qrs_index-1));
% Recording mean Pra, Pa, and qpv to the num matrix
% and varying Pv and/or Pa every fb_int beats until
% all of the pre-determined variations have been
% implemented.
if (qrs_index == fb_int & countn < lengtht)
countn = countn+1;
countj = countj+1;
num(1,countn) = mrap;
num(2,countn) = P(2,k);
num(3,countn) = mco;
if (countj ~= lengthra)
th(30) = rarange(countj+1);
end
if (countj == lengthra & countn ~= lengtht)
th(29) = th(29)+th(90);
P(2,k) = th(29);
% Adjusting range of Pv for large Pa values.
if (th(29) > 180)
tmp = tmp+1;
final = final-0.5*tmp;
rarangei = final:-th(88):0;
while (length(rarangei) < lengthra-1)
rarangei = [rarangei, rarangei(length(rarangei))-th(88)];
end
rarange = [rarangei, th(23)+1];
if (lengthra == 2);
lengthra = 1;
rarange = final;
end
th(30) = final;
end
fb_int = fb_int+5;
countj = 0;
P(3,k) = rarange(countj+1);
elseif (countn ~= lengtht)
P(3,k) = th(30);
fb_int = fb_int+5;
end
end
% Resetting the initial values for the integration
% of Pra and qpv over the next cardiac cycle.
summrap = 0.5*sumdeltaT*(P(9,k-1)+P(9,k));
summco = 0.5*sumdeltaT*(q(1,k-1)+qpv);
end
if (preparation == 2)
% Calculating the integral of Pra and qv
% precisely over the previous cardiac cycle.
summrap = summrap-0.5*deltaT*(P(9,k-1)+P(9,k));
summrap = summrap+0.5*(deltaT-sumdeltaT)*(P(9,k-1)+P(9,k));
summco = summco-0.5*deltaT*(q(4,k-1)+qv);
summco = summco+0.5*(deltaT-sumdeltaT)*(q(4,k-1)+qv);
% Computing mean Pra and ql over the previous cardiac cycle.
mrap = summrap/(t(k)-sumdeltaT-ap(7,qrs_index-1));
mco = summco/(t(k)-sumdeltaT-ap(7,qrs_index-1));
% Recording mean Pra and ql to the num matrix
% and varying Crd after every fb_int beats
% until all the pre-determined variations have
% been implemented.
if (qrs_index == fb_int & Crd <= 60 & lengthra ~= 1)
countn = countn+1;
num(1,countn)=mrap;
num(2,countn)=mco;
Crdo = Crd;
Crd = Crd+th(89);
% Adjustment of Pr in response to variation in Crd -- for
% constant volume considerations.
th(6) = (th(32)/(th(27)-th(12)))*Crd;
ypres = (P(4,k)-P(7,k))/th(32);
xvol = vent_vol(Q(4,k-1),ypres,1/Crdo);
alpha = 1/th(6);
x0 = 1/(alpha+1);
kscale = 50;
c1 = 1+exp(-kscale*(1-x0));
d1 = 1+exp(kscale*x0);
denb1 = (1+(1/kscale)*log(c1/d1));
beta = (1-alpha)/denb1;
c = 1+exp(-kscale*(xvol-x0));
d = 1+exp(kscale*x0);
denb = (xvol+(1/kscale)*log(c/d));
P(4,k) = th(32)*(alpha*xvol+beta*denb)+P(7,k);
% Incrementing fb_int.
fb_int = fb_int+5;
end
% Resetting the initial values for the integration
% of Pra and qv over the next cardiac cycle.
summrap = 0.5*sumdeltaT*(P(9,k-1)+P(9,k));
summco = 0.5*sumdeltaT*(q(4,k-1)+qv);
end
end
% Recording the adjustable parameter values at the current integration step.
ap(1:6,k) = [th(1)*((th(26)-th(9))/th(31)) th(5)*((th(27)-th(12))/th(32)) th(11) th(16) th(22) .3*sqrt(th(24))]';
% Recording the variable elastance values at the current integration step.
if (preparation ~= 10)
[El,dEl] = var_cap(th(1),th(2),th(24),sumdeltaT);
[Er,dEr] = var_cap(th(5),th(6),th(24),sumdeltaT);
ve(1:2,k) = [El*(th(31)/(th(26)-th(9))); Er*(th(32)/(th(27)-th(12)))];
else
[El,dEl] = var_vcap(th(1),th(2),th(24),sumdeltaT);
[Er,dEr] = var_vcap(th(5),th(6),th(24),sumdeltaT);
[Ela,dEla] = var_acap(th(80),th(81),th(24),sumdeltaT);
[Era,dEra] = var_acap(th(84),th(85),th(24),sumdeltaT);
ve(:,k) = [El*(th(31)/(th(26)-th(9))); Er*(th(32)/(th(27)-th(12))); Ela; Era];
end
% Obtaining the volume values at the current integration step. Equations depend
% on the preparation being implemented.
if (preparation ~= 10)
Q(1:6,k) = [(1/El) th(3) th(4) (1/Er) th(7) th(8)]'.*(P(1:6,k)-P(7,k)*[1 0 0 1 1 1]') + [th(9:14)];
else
Q(1:6,k) = [(1/El) th(3) th(4) (1/Er) th(7) th(8)]'.*([P(1:6,k)]-P(7,k)*[1 0 0 1 1 1]')+[th(9:14)];
Q(8:9,k) = [(1/Era) (1/Ela)]'.*([P(10:11,k)]-P(7,k)*[1 1]')+[th(86); th(82)];
end
if (preparation == 6)
Q(2,k) = th(72)*P(2,k)+th(74);
Q(8,k) = th(73)*P(10,k)+th(75);
elseif (preparation == 1)
Q(2,k) = 0;
Q(3,k) = 0;
elseif (preparation == 3)
Q(2,k) = 0;
Q(6,k) = 0;
elseif (preparation == 7)
alphatau = -th(70)/th(3);
Qamax = th(71)-th(3)^2/th(70);
Qao = th(71)+(th(3)^2/th(70))*exp(-th(70)*th(40)/th(3))*(1-exp(th(70)*th(40)/th(3)));
Q(2,k) = (Qamax-Qao)*(1-exp(-alphatau*P(2,k)))+Qao;
elseif (preparation == 8)
alphatau = -th(70)/th(72);
Qamax = th(76)-th(72)^2/th(70);
Qao = th(76)+(th(72)^2/th(70))*exp(-th(70)*th(40)/th(72))*(1-exp(th(70)*th(40)/th(72)));
Q(2,k) = (Qamax-Qao)*(1-exp(-alphatau*P(2,k)))+Qao;
Q(8,k) = th(73)*P(10,k)+th(75);
end
if (preparation ~= 5)
x = (P(1,k) < P(7,k));
mbintscalar(x);
if (x)
Q(1,k) = th(9)*(2/pi)*atan(((P(1,k)-P(7,k))/((2/pi)*th(9)*El))) + th(9);
else
yl = (P(1,k)-P(7,k))/th(31);
xl = vent_vol(Q(1,k-1),yl,El);
Q(1,k) = (th(26)-th(9))*xl+th(9);
end
x = (P(4,k) < P(7,k));
mbintscalar(x);
if (x)
Q(4,k) = th(12)*(2/pi)*atan(((P(4,k)-P(7,k))/((2/pi)*th(12)*Er))) + th(12);
else
yr = (P(4,k)-P(7,k))/th(32);
xr = vent_vol(Q(4,k-1),yr,Er);
Q(4,k) = (th(27)-th(12))*xr+th(12);
end
end
% Correcting Qv and Pv in order to keep total blood volume in the
% intact circulatory preparations, which may vary due to integration
% error, constant.
if (preparation == 0 | preparation == 4 | preparation == 5 | preparation == 6 | preparation == 7 | preparation == 8 | preparation == 9 | preparation == 10)
if (preparation == 6 | preparation == 8)
Qtot = sum(Q([1:6 8],k));
elseif (preparation == 10)
Qtot = sum(Q([1:6 8 9],k));
else
Qtot = sum(Q(1:6,k));
end
Q(3,k) = Q(3,k) - (Qtot-th(21));
P(3,k) = (Q(3,k)-th(11))/th(4);
end
% Checking for a request to pause simulation after the completion
% of every beat provided that waveforms are being viewed on-line.
if (strcmp(signals,'-1') == 0 & ipfmcond)
flagnew = read_key;
if (flagnew == 1)
% Re-reading parameter file, in the event of a possible update.
thspold = thsp;
thold = th;
fidp = fopen(parameterfile,'r');
[th,signals] = read_param(fidp);
fclose(fidp);
% Precluding adjustment of unchangeable parameters.
th(98) = thold(98); % Qfr
th(33) = thold(33); % Pms
th(88) = thold(88); % Pvs
th(89) = thold(89); % Crds
th(90) = thold(90); % Pas
% Precluding adjustment of Pv, Pa, and Crd if they are being
% automatically varied.
if (preparation == 1 & lengthra > 1)
th(30) = thold(30);
end
if (preparation == 1 & lengthsa > 1)
th(29) = thold(29);
end
if (preparation == 2 & lengthra > 1)
th(6) = thold(6);
end
% Setting annotator flag.
if (th(105) == 0)
annotator = 'wqrs';
else
annotator = 'qrs';
end
% Determining if parameter update is relevant to
% status of current simulation.
parameter_change = param_change(th,thold,flag);
mbintscalar(parameter_change == 1);
if (parameter_change == 1)
% Annotating and documenting the relevant parameter update.
count_pc = count_pc+1;
ap(9,qrs_index) = count_pc;
write_param(parameterfile,num2str(count_pc));
% Adjusting setpoint pressures.
thsp(1) = th(40);
thsp(7) = th(41);
% Conserving Qfr by altering Pth at the functional reserve
% volume of the lungs.
th(23) = th(25)-((th(98)-th(102))/th(101));
% Updating respiratory-related pressures and volumes
% over the remainder of the integration period, if necessary.
% No breathing.
if (breathing == 0 & th(23) ~= thold(23))
thbreathe = resp_act(th,at,deltaT/2,0,(tics-1)*deltaT);
% Fixed-rate breathing.
elseif (breathing == 1 & (th(42) ~= thold(42) | th(77) ~= thold(77) | th(23) ~= thold(23) | th(100) ~= thold(100)))
atcount = 2;
x = (at(atcount) < t(k)+1.5);
mbintscalar(x);
while (x)
atcount=atcount+1;
x = (at(atcount) < t(k)+1.5);
mbintscalar(x);
end
start_index = k+ceil((at(atcount)-t(k))/deltaT);
start_time = t(k)+deltaT*ceil((at(atcount)-t(k))/deltaT)-at(atcount);
end_time = (length(thbreathe)-(2*(start_index-1)+1))*(deltaT/2)+start_time;
if (start_index < tics)
at = [at(1:atcount); zeros(length(start_index:tics)+1,1)];
x = (at(atcount) <= deltaT*tics+1);
mbintscalar(x);
while (x)
atcount=atcount+1;
at(atcount) = at(atcount-1)+th(42);
x = (at(atcount) <= deltaT*tics+1);
mbintscalar(x);
end
at = at(1:atcount);
thbreathe(:,2*(start_index-1)+1:end) = resp_act(th,zeros(4,1),deltaT/2,start_time,end_time);
ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran);
ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];
end
% Random-interval breathing.
elseif (breathing == 2 & (th(42) ~= thold(42) | th(77) ~= thold(77) | th(99) ~= thold(99) | th(23) ~= thold(23) | th(100) ~= thold(100)))
atcount = 2;
x = (at(atcount) < t(k)+1.5);
mbintscalar(x);
while (x)
atcount=atcount+1;
x = (at(atcount) < t(k)+1.5);
mbintscalar(x);
end
start_index = k+ceil((at(atcount)-t(k))/deltaT);
start_time = t(k)+deltaT*ceil((at(atcount)-t(k))/deltaT)-at(atcount);
end_time = (length(thbreathe)-(2*(start_index-1)+1))*(deltaT/2)+start_time;
if (start_index < tics)
thbreathe(:,2*(start_index-1)+1:end) = resp_act(th,[1; (at(atcount+1:end)-at(atcount))],deltaT/2,start_time,end_time);
ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran);
ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];
end
end
% Implementing a step change to Qlu, if mandated.
if (th(36) ~= 0)
if (breathing == 1 | breathing == 2)
atcount = 2;
x = (at(atcount) < t(k)+1.5);
mbintscalar(x);
while (x)
atcount=atcount+1;
x = (at(atcount) < t(k)+1.5);
mbintscalar(x);
end
start_index = k+ceil((at(atcount)-t(k))/deltaT);
start_time = (start_index-1)*deltaT;
th(103) = start_time+1.5;
elseif (breathing == 0)
th(103) = t(k)+1.5;
end
if (th(103) < (tics-1)*deltaT)
thbreathes = resp_act(th,[3; th(36)],deltaT/2,0,(tics-1)*deltaT);
dthbreathe = [-th(98)*ones(1,length(thbreathes)); -th(23)*ones(1,length(thbreathes)); zeros(2,length(thbreathes))];
thbreathe = thbreathe+thbreathes+dthbreathe;
ilv = ilv_dec(thbreathe(1,:),deltaT/2,Sgran);
ilv = [th(98)*ones(length(hilvhr)-(1.5/Sgran),1); ilv'; th(98)*ones(1.5/Sgran,1)];
end
end
% Assigning values to Pth, Ppac (Palv), and Qlu
P(7,k) = thbreathe(2,2*(k-1)+1);
P(8,k) = thbreathe(4,2*(k-1)+1);
Q(7,k) = thbreathe(1,2*(k-1)+1);
% Conserving volume in each compartment of the pulsatile
% heart and circulation through pressure adjustments.
P(1:6,k) = conserve_vol(th,sumdeltaT,Q(1:6,k),P(7,k));
if (preparation == 1)
P(2,k) = th(29);
P(3,k) = th(30);
end
% Updating adjustable parameters via instantaneous values
% and setpoint values.
if (thspold(3) ~= th(16))
thn(16) = th(16);
thp(16) = th(16);
thsp(3) = th(16);
ap(4,k) = th(16);
else
th(16) = thold(16);
end
if (thspold(2) ~= th(22))
thn(22) = th(22);
thp(22) = th(22);
thsp(2) = th(22);
ap(5,k) = th(22);
else
th(22) = thold(22);
end
if (thspold(6) ~= th(11))
thn(11) = th(11);
thp(11) = th(11);
thsp(6) = th(11);
ap(3,k) = th(11);
else
th(11) = thold(11);
if (preparation == 0)
P(3,k) = (Q(3,k)-th(11))/th(4);
end
end
thn([1:10 12:15 17:21 23:end]) = th([1:10 12:15 17:21 23:end]);
thp([1:10 12:15 17:21 23:end]) = th([1:10 12:15 17:21 23:end]);
if (thspold(4) ~= th(1))
thsp(4) = th(1);
else
th(1) = thold(1);
end
if (thspold(5) ~= th(5))
thsp(5) = th(5);
else
th(5) = thold(5);
end
end
end
end
% Obtaining the flow rate values at the current integration step. Equations
% are dependent on the preparation being implemented.
if (preparation == 0 | preparation == 1 | preparation == 2 | preparation == 4 | preparation == 5 | preparation == 6 | preparation == 7 | preparation == 8)
x = (P(6,k) > P(1,k));
mbintscalar(x);
if (x)
q(1,k) = (P(6,k)-P(1,k))/th(20);
else
q(1,k) = 0;
end
elseif (preparation == 3)
x = (th(34) > P(1,k));
mbintscalar(x);
if (x)
q(1,k) = (th(34)-P(1,k))/th(20);
else
q(1,k) = 0;
end
elseif (preparation == 10)
x = (P(11,k) > P(1,k));
mbintscalar(x);
if (x)
q(1,k) = (P(11,k)-P(1,k))/th(83);
else
q(1,k) = 0;
end
end
if (preparation == 0 | preparation == 1 | preparation == 4 | preparation == 5 | preparation == 6 | preparation == 7 | preparation == 8 | preparation == 9 | preparation == 10)
x = (P(1,k) > P(2,k));
mbintscalar(x);
if (x)
q(2,k) = (P(1,k)-P(2,k))/th(15);
else
q(2,k) = 0;
end
elseif (preparation == 2)
x = (P(1,k) > th(33));
mbintscalar(x);
if (x)
q(2,k) = (P(1,k)-th(33))/th(15);
else
q(2,k) = 0;
end
elseif (preparation == 3)
x = (P(1,k) > th(29));
mbintscalar(x);
if (x)
q(2,k) = (P(1,k)-th(29))/th(15);
else
q(2,k) = 0;
end
end
if (preparation == 0 | preparation == 2 | preparation == 4 | preparation == 5 | preparation == 7 | preparation == 9 | preparation == 10)
q(3,k) = (P(2,k)-P(3,k))/th(16);
elseif (preparation == 1 | preparation == 3)
q(3,k) = 0;
elseif (preparation == 6 | preparation == 8)
q(3,k) = (P(10,k)-P(3,k))/th(16);
end
if (preparation ~= 10 & preparation ~= 5)
x = (P(3,k) > P(4,k) & P(4,k) >= th(91));
y = (P(3,k) > th(91) & th(91) > P(4,k));
mbintscalar(x);
mbintscalar(y);
if (x)
q(4,k) = (P(3,k)-P(4,k))/th(17);
elseif (y)
q(4,k) = (P(3,k)-th(91))/th(17);
else
q(4,k) = 0;
end
elseif (preparation == 5)
x = (P(3,k) >= P(4,k));
mbintscalar(x);
if (x)
q(4,k) = (P(3,k)-P(4,k))/th(17);
else
q(4,k) = 0;
end
else
x = (P(9,k) >= P(4,k));
mbintscalar(x);
if (x)
q(4,k) = (P(9,k)-P(4,k))/th(87);
else
q(4,k) = 0;
end
end
if (preparation ~= 2)
x = (P(4,k) > P(5,k));
mbintscalar(x);
if (x)
q(5,k) = (P(4,k)-P(5,k))/th(18);
else
q(5,k) = 0;
end
elseif (preparation == 2)
x = (P(4,k) > th(28));
mbintscalar(x);
if (x)
q(5,k) = (P(4,k)-th(28))/th(18);
else
q(5,k) = 0;
end
end
if (preparation == 2 | preparation == 5)
q(6,k) = (P(5,k)-P(6,k))/th(19);
elseif (preparation == 0 | preparation == 1 | preparation == 4 | preparation == 6 |preparation == 7 | preparation == 8 | preparation == 9 | preparation == 10)
h1 = (1330/(1.06*980))*(P(6,k)-P(8,k));
x = h1 < -10;
mbintscalar(x);
y = h1 > 20;
mbintscalar(y);
if (x)
h1 = -10;
elseif (y)
h1 = 20;
end
h2 = (1330/(1.06*980))*(P(5,k)-P(8,k));
x = h2 < -10;
mbintscalar(x);
y = h2 > 20;
mbintscalar(y);
if (x)
h2 = -10;
elseif (y)
h2 = 20;
end
q(6,k) = (((P(5,k)-P(6,k))/(30*th(19)))*(h1+10)) + (((P(5,k)-P(8,k))/(30*th(19)))*(h2-h1)) - (((1.06*980)/(2660*(30*th(19))))*(h2^2-h1^2));
elseif (preparation == 3)
q(6,k) = 0;
end
if (preparation == 10)
x = (P(3,k) >= th(91) & P(9,k) >= th(91));
y = (P(3,k) > th(91) & th(91) > P(9,k));
z = (P(9,k) > th(91) & th(91) > P(3,k));
mbintscalar(x);
mbintscalar(y);
mbintscalar(z);
if (x)
q(7,k) = (P(3,k)-P(9,k))/th(17);
elseif (y)
q(7,k) = (P(3,k)-th(91))/th(17);
elseif (z)
q(7,k) = (th(91)-P(9,k))/th(17);
else
q(7,k) = 0;
end
elseif (preparation == 6 | preparation == 8)
q(7,k) = P(11,k);
end
if (preparation == 10)
q(8,k) = (P(6,k)-P(11,k))/th(20);
end
% Writing waveforms and annotations to files in MIT format every step parameter
% seconds once a window parameter duration of data have been calculated.
if (strcmp(signals,'-1') == 0 & t(k) > th(104) & (t(k)-t(kstart)) > th(78))
% Writing newly generated waveform data to MIT format files.
X = [P(1:9,kstart:k); Q(1:6,kstart:k); Q(7,kstart:k)*0.1; q(1:6,kstart:k); ap(1:2,kstart:k)*100; ap(3,kstart:k); ap(4,kstart:k)*100; ap(5,kstart:k)*600; ve(1:2,kstart:k)*10]*10;
fwrite(fid,X,'short');
% Updating variables necessary for writing waveforms.
kstart = k+1;
start_time = t(k)-th(104);
% Writing newly generated annotations to MIT format files.
qrso = [];
qrso = diff(ap(8,qrs_index_start:qrs_index))';
for i = 1:length(qrso)
fwrite(fid2,qrso(i),'bit10');
fwrite(fid2,1,'ubit6');
% If parameter update has been made, writing auxillary information.
x = (ap(9,qrs_index_start+i) ~= 0);
mbintscalar(x)
if (x)
suffix = num2str(ap(9,qrs_index_start+i));
fwrite(fid2,Nbytes+(length(suffix)-1),'bit10');
fwrite(fid2,63,'ubit6');
fwrite(fid2,[parameterfile '.' suffix],[num2str(Nbytes+length(suffix)-1) '*uchar']);
end
end
qrs_index_start = qrs_index;
% Alternating outputfile signaled to WAVE in order to permit
% on-line updating of the annotation flag.
if (wrflag == 0)
wave_remote( outputfile, annotator, start_time, signals );
wrflag = 1;
else
wave_remote( outputfile2, annotator, start_time, signals );
wrflag = 0;
end
end
% If cardiac function/venous return curves are complete, there
% is no need for further integration. That is, for-loop is exited.
if (preparation == 1 | preparation == 2)
if (countn >= length(num'))
tics = k;
break;
end
end
end
% Writing remaining portion of simulated data to the MIT format files.
if (strcmp(signals,'-1') == 0)
% Writing last portion of generated waveform data.
X = [P(1:9,kstart:tics); Q(1:6,kstart:tics); Q(7,kstart:tics)*0.1; q(1:6,kstart:tics); ap(1:2,kstart:tics)*100; ap(3,kstart:tics); ap(4,kstart:tics)*100; ap(5,kstart:tics)*600; ve(1:2,kstart:tics)*10]*10;
fwrite(fid,X,'short');
% Writing last annotations.
qrso = [];
qrso = diff(ap(8,qrs_index_start:qrs_index))';
for i = 1:length(qrso)
fwrite(fid2,qrso(i),'bit10');
fwrite(fid2,1,'ubit6');
% If parameter updates have been made, writing auxillary info.
x = (ap(9,qrs_index_start+i) ~= 0);
mbintscalar(x)
if (x)
suffix = num2str(ap(9,qrs_index_start+i));
fwrite(fid2,Nbytes+(length(suffix)-1),'bit10');
fwrite(fid2,63,'ubit6');
fwrite(fid2,[parameterfile '.' suffix],[num2str(Nbytes+length(suffix)-1) '*uchar']);
end
end
fwrite(fid2,0,'bit10');
fwrite(fid2,0,'ubit6');
fclose(fid2);
fclose(fid);
end
% If MATLAB execution, writing qrs vector as times of ventricular contractions.
% If Linux execution, writing qrs vector as samples of ventricular contractions.
if (nargin == 2)
qrs = ap(7,1:qrs_index);
else
qrs = ap(8,1:qrs_index);
end
% Only supplying as output the pressures that are generated by
% the model, for MATLAB execution only.
if (preparation == 6 | preparation == 8)
P = P(1:10,:);
q = q(1:7,:);
Q = Q(1:8,:);
ve = ve(1:2,:);
elseif (preparation == 10 | preparation == 13)
P = P([1:9 11],:);
else
P = P(1:9,1:tics);
q = q(1:6,1:tics);
Q = Q(1:7,1:tics);
ve = ve(1:2,1:tics);
ap = ap(:,1:tics);
t = t(1:tics);
end