% The function resp_act.m computes respiratory-related waveforms
% over the entire integration period.
%
% Function arguments:
% th - current parameter values
% at - column vector containing the arrival times of each
% respiratory cycle
% deltaT - half of the integration step size (s)
% stime - start time
% etime - end time
%
% Function outputs:
% thbreathe - 4xtime vector containing respiratory-related data
% - [Qlu; Pth; dPth; Ppac], where time is length of
% the vector stime:deltaT:etime
%
function thbreathe = resp_act(th,at,deltaT,stime,etime)
% Establishing times in which the respiratory-related waveforms
% are to be computed.
time = stime:deltaT:etime;
% Pre-allocating memory for the function outputs.
thbreathe = zeros(4,length(time));
% No breathing.
if (length(at) == 1)
% Setting respiratory-related waveforms to their functional
% reserve values.
thbreathe(1,:) = th(98)*ones(1,length(time));
thbreathe(2,:) = th(23)*ones(1,length(time));
thbreathe(3,:) = zeros(1,length(time));
thbreathe(4,:) = th(25)*ones(1,length(time));
% Step change in breathing.
elseif (length(at) == 2)
% Parametrization for step Qlu function -- hyperbolic tangent.
% Width of transition band of step with small values
% indicating a small transition band.
tau = 0.25;
% Height of step.
sigmailv = at(2);
time1 = stime:deltaT:th(103);
time2 = time1(end)+deltaT:deltaT:etime;
Q = (sigmailv./(1+exp(-(time-th(103))/tau)))+th(98);
dQ = ((sigmailv/tau)*exp(-(abs(time-th(103)))/tau))./((1+exp(-(abs(time-th(103)))/tau)).^2);
d2Q1 = -((sigmailv/tau^2)*exp(-(abs(time1-th(103)))/tau).*(exp(-(abs(time1-th(103)))/tau)-1))./((1+exp(-(abs(time1-th(103)))/tau)).^3);
d2Q2 = ((sigmailv/tau^2)*exp(-(abs(time2-th(103)))/tau).*(exp(-(abs(time2-th(103)))/tau)-1))./((1+exp(-(abs(time2-th(103)))/tau)).^3);
d2Q = [d2Q1 d2Q2];
thbreathe(1,:) = Q;
thbreathe(2,:) = th(91)-th(100)*dQ-(1/th(101))*(Q-th(102));
thbreathe(3,:) = -th(100)*d2Q-(1/th(101))*dQ;
thbreathe(4,:) = th(91)-th(100)*dQ;
% Impulse change in breathing.
elseif (length(at) == 3)
% Parametrization for impulse Qlu function -- derivative of
% hyperbolic tangent.
% Width of transition band of step with small values
% indicating a small transition band.
tau = 0.25;
% Area of impulse.
sigmailv = at(2);
time1 = stime:deltaT:th(103);
time2 = time1(end)+deltaT:deltaT:etime;
Q = ((sigmailv/tau)*exp(-(abs(time-th(103)))/tau))./((1+exp(-(abs(time-th(103)))/tau)).^2)+th(98);
dQ1 = -((sigmailv/tau^2)*exp(-(abs(time1-th(103)))/tau).*(exp(-(abs(time1-th(103)))/tau)-1))./((1+exp(-(abs(time1-th(103)))/tau)).^3);
dQ2 = ((sigmailv/tau^2)*exp(-(abs(time2-th(103)))/tau).*(exp(-(abs(time2-th(103)))/tau)-1))./((1+exp(-(abs(time2-th(103)))/tau)).^3);
dQ = [dQ1 dQ2];
d2Q = (sigmailv/tau^3)*(exp(-3*(abs(time-th(103)))/tau)-4*exp(-2*(abs(time-th(103)))/tau)+exp(-(abs(time-th(103)))/tau))./((1+exp(-(abs(time-th(103)))/tau)).^4);
thbreathe(1,:) = Q;
thbreathe(2,:) = th(91)-th(100)*dQ-(1/th(101))*(Q-th(102));
thbreathe(3,:) = -th(100)*d2Q-(1/th(101))*dQ;
thbreathe(4,:) = th(91)-th(100)*dQ;
else
flag = at(1);
at = at(2:length(at));
% Fixed-rate breathing -- sinusoidal Qlu variations plus
% functional reserve volume.
if (flag == 0)
% Calculating the respiratory-related waveforms all at once.
thbreathe(1,:) = (th(77)/2)*(1-cos(2*pi*time/th(42)))+th(98);
thbreathe(2,:) = th(91)-th(100)*(th(77)*pi/th(42))*sin(2*pi*time/th(42)) - (1/th(101))*(thbreathe(1,:)-th(102));
thbreathe(3,:) = -th(100)*(th(77)/2)*(2*pi/th(42))^2*cos(2*pi*time/th(42))-(1/th(101))*(th(77)*pi/th(42))*sin(2*pi*time/th(42));
thbreathe(4,:) = th(91)-th(100)*(th(77)*pi/th(42))*sin(2*pi*time/th(42));
% Random-interval breathing -- each Qlu respiratory cycle is one period
% of a sinusoid plus the functional reserve volume, but the duration of
% the periods are random.
else
% Prior to initiation of first respiratory cycle. Setting
% respiratory-related waveforms to functional reserve values.
startm = 1;
m = 1;
while (time(m) < at(1))
m=m+1;
end
endm=m-1;
thbreathe(1,startm:endm) = th(98)*ones(1,endm);
thbreathe(2,startm:endm) = th(23)*ones(1,endm);
thbreathe(3,startm:endm) = zeros(1,endm);
thbreathe(4,startm:endm) = th(25)*ones(1,endm);
% After initiation of first respiratory cycle, calculating
% respiratory-related waveforms one respiratory cycle at
% a time.
for counting = 2:length(at)
startm = endm+1;
m = startm;
while (m <= length(time) & time(m) < at(counting))
m=m+1;
end
endm=m-1;
period = at(counting)-at(counting-1);
% Establishing tidal volume for each respiratory cycle.
% Constant instantaneous alveolar ventilation rate
% set to mean rate of fixed-rate breathing (ml/s).
mbrealscalar(flag == 1);
if (flag == 1)
meanavr = (th(77)-th(99))/th(42);
Qtr = meanavr*period + th(99);
% Constant tidal volume resulting in a nearly constant
% instantaneous alveolar ventilation rate.
else
Qtr = 236.36*(period^(1/2));
end
thbreathe(1,startm:endm) = (Qtr/2)*(1-cos(2*pi*(time(startm:endm)-at(counting-1))/period))+th(98);
thbreathe(2,startm:endm) = th(91)-th(100)*(Qtr*pi/period)*sin(2*pi*(time(startm:endm)-at(counting-1))/period) - (1/th(101))*(thbreathe(1,startm:endm)-th(102));
thbreathe(3,startm:endm) = -th(100)*(Qtr/2)*(2*pi/period)^2*cos(2*pi*(time(startm:endm)-at(counting-1))/period)-(1/th(101))*(Qtr*pi/period)*sin(2*pi*(time(startm:endm)-at(counting-1))/period);
thbreathe(4,startm:endm) = th(91)-th(100)*(Qtr*pi/period)*sin(2*pi*(time(startm:endm)-at(counting-1))/period);
end
end
end