% The function third_nac_init_cond.m computes the initial values of the % intact circulation with third-order systemic arteries and nonlinear % large, elastic arterial compliance. The initial values are computed % based on conservation laws of a linearized model. % % Function arguments: % th - vector containing the initial parameter values % % Function outputs: % P - a 8x1 vector containing the six initial pressures % = [Pl(0); Pe(0); Pv(0); Pr(0); Ppa(0); Ppv(0); Pm(0); qe(0)] % Q - a 7x1 vector containing the six initial volumes % = [Ql(0); Qe(0); Qv(0); Qr(0); Qpa(0); Qpv(0); Qm(0)] % q - a 7x1 vector containing the six initial flow rates % = [qpv(0); ql(0); qm(0); qv(0); qr(0); qpa(0); qe(0)] % ve - a 2x1 vector containing the initial variable elastance values % = [El(0); Er(0)]; % function [P,Q,q,ve] = third_nac_init_cond(th) % Storing nonlinear ventricular compliance values for % subsequent initial ventricular volume calculation. Cld = th(2); Crd = th(6); % Converting ventricular compliances to linear values which % is necessary for estimating initial pressures. 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)); % Estimating initial pressures as well as qe which is also a state variable. Ts = .3*sqrt(1/th(22)); Td = 1/th(22) - Ts; temp = -th(2)*th(23) + th(1)*th(23); b = [temp+th(6)*th(23)-th(5)*th(23); temp; temp; temp; temp; temp; temp; th(21)-(th(14)+th(13)+th(12)+th(11)+th(10)+th(9))+th(2)*th(23)+th(6)*th(23)+th(7)*th(23)+th(8)*th(23)]; A = [th(1), -th(2), 0, 0, -th(5), th(6), 0, 0; th(1)+(Ts/th(15)), -th(2), -Ts/th(15), 0, 0, 0, 0, 0; th(1), -th(2), 1/(th(22)*th(16)), -1/(th(22)*th(16)), 0, 0, 0, 0; th(1), -th(2), 0, Td/th(17), 0, -Td/th(17), 0, 0; th(1), -th(2), 0, 0, Ts/th(18), 0, -Ts/th(18), 0; th(1), -th(2), 0, 0, 0, 0, 1/(th(22)*th(19)), -1/(th(22)*th(19)); th(1), -(th(2)+Td/th(20)), 0, 0, 0, 0, 0, Td/th(20); 0, th(2), (th(72)+th(73)), th(4), 0, th(6), th(7), th(8)]; x = A\b; P = [x(2:4)' x(6:8)' x(3) 0]'; % Establishing initial flow rates. q = zeros(7,1); if (P(6) > P(1)) q(1) = (P(6)-P(1))/th(20); else q(1) = 0; end if (P(1) > P(2)) q(2) = (P(1)-P(2))/th(15); else q(2) = 0; end q(3) = (P(7)-P(3))/th(16); if (P(3) > P(4)) q(4) = (P(3)-P(4))/th(17); else q(4) = 0; end if (P(4) > P(5)) q(5) = (P(4)-P(5))/th(18); else q(5) = 0; end q(6) = (P(5)-P(6))/th(19); % Establishing initial variable ventricular elastance values. ve = [1/th(2) 1/th(6)]'; % Establishing initial volumes. Q = [th(2) th(72) th(4) th(6) th(7) th(8) th(73)]'.*(P(1:7)-th(23)*[1 0 0 1 1 1 0]') + [th(9); th(74); th(11:14); th(75)]; if (P(1) < th(23)) Q(1) = th(9)*(2/pi)*atan(((P(1)-th(23))/((2/pi)*th(9)*ve(1)))) + th(9); else yl = (P(1)-th(23))/th(31); xl = vent_vol(0.5,yl,1/Cld); Q(1) = (th(26)-th(9))*xl+th(9); end if (P(4) < th(23)) Q(4) = th(12)*(2/pi)*atan(((P(4)-th(23))/((2/pi)*th(12)*ve(2)))) + th(12); else yr = (P(4)-th(23))/th(32); xr = vent_vol(0.5,yr,1/Crd); Q(4) = (th(27)-th(12))*xr+th(12); end 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))); v(2) = (Qamax-Qao)*(1-exp(-alphatau*P(2)))+Qao;