% The function sc_eval_deriv.m evaluates the right-hand side of a % set of first-order differential equations which characterize the % systemic circulation at a desired time ste[. % % Function arguments: % P - a 6x1 vector containing the current pressure values % = [Pl; Pa; Pv; Pr; Ppa; Ppv] % Qc - a 2x1 vector containing the current ventricular volume values % - [Ql; Qr] % Pbreathe - 4x1 vector containing current respiratory values % [Qlu; Pth; dPth; Ppac (Palv)] % th - current parameter values % sumdeltaT - time surpassed in current cardiac cycle % % Function outputs: % dP - a 6x1 vector containing the derivative of current % pressure values % = [dPl; dPa; dPv; dPr; dPpa; dPpv] % function dP = sc_eval_deriv(P,Qc,Pbreathe,th,sumdeltaT); % Computing flow rates from the current pressure values. if (P(6) >= P(1)) qli = (P(6)-P(1))/th(20); else qli = 0; end if (P(1) >= th(33)) qlo = (P(1)-th(33))/th(15); else qlo = 0; end qa = (P(2)-P(3))/th(16); if (P(3) > P(4) & P(4) >= th(91)) qri = (P(3)-P(4))/th(17); elseif (P(3) > th(91) & th(91) > P(4)) qri = (P(3)-th(91))/th(17); else qri = 0; end qlovr = qri; if (P(4) >= th(28)) qro = (P(4)-th(28))/th(18); else qro = 0; end h1 = (1330/(1.06*980))*(P(6)-Pbreathe(4)); if (h1 < -10) h1 = -10; elseif (h1 > 20) h1 = 20; end h2 = (1330/(1.06*980))*(P(5)-Pbreathe(4)); if (h2 < -10) h2 = -10; elseif (h2 > 20) h2 = 20; end qpa = (((P(5)-P(6))/(30*th(19)))*(h1+10)) + (((P(5)-Pbreathe(4))/(30*th(19)))*(h2-h1)) - (((1.06*980)/(2660*(30*th(19))))*(h2^2-h1^2)); % Computing ventricular elastances. [El,dEl] = var_cap(th(1),th(2),th(24),sumdeltaT); [Er,dEr] = var_cap(th(5),th(6),th(24),sumdeltaT); % Computing dP from the above calculations and th parameter vector. if (P(1) < Pbreathe(2)) g = ((P(1)-Pbreathe(2)))/(El*th(9)*(2/pi)); dPl = El*(qli-qlo)*(1+g^2) + (P(1)-Pbreathe(2))*(dEl/El) + Pbreathe(3); else y = (P(1)-Pbreathe(2))/th(31); x = vent_vol(Qc(1),y,El); x0 = 1/(El+1); k = 50; mbrealscalar(-k*(x-x0)); mbrealscalar(-k*(1-x0)); mbrealscalar(k*x0); dx = (qli-qlo)/(th(26)-th(9)); dx0 = -dEl/((El+1)^2); c = 1+exp(-k*(1-x0)); d = 1+exp(k*x0); g = 1+exp(-k*(x-x0)); dc = k*dx0*exp(-k*(1-x0)); dd = k*dx0*exp(k*x0); dg = -k*(dx-dx0)*exp(-k*(x-x0)); denb = (1+(1/k)*log(c/d)); b = (1-El)/denb; db = (-dEl*denb - (1-El)*(1/k)*((d*dc-c*dd)/(c*d)))/(denb^2); dy = dEl*x + El*dx + db*(x+(1/k)*(log(g/d))) + b*(dx+(1/k)*((d*dg-g*dd)/(g*d))); dPl = th(31)*dy + Pbreathe(3); end if (P(4) < Pbreathe(2)) g = ((P(4)-Pbreathe(2)))/(Er*th(12)*(2/pi)); dPr = Er*(qri-qro)*(1+g^2) + (P(4)-Pbreathe(2))*(dEr/Er) + Pbreathe(3); else y = (P(4)-Pbreathe(2))/th(32); x = vent_vol(Qc(2),y,Er); x0 = 1/(Er+1); k = 50; mbrealscalar(-k*(x-x0)); mbrealscalar(-k*(1-x0)); mbrealscalar(k*x0); dx = (qri-qro)/(th(27)-th(12)); dx0 = -dEr/((Er+1)^2); c = 1+exp(-k*(1-x0)); d = 1+exp(k*x0); g = 1+exp(-k*(x-x0)); dc = k*dx0*exp(-k*(1-x0)); dd = k*dx0*exp(k*x0); dg = -k*(dx-dx0)*exp(-k*(x-x0)); denb = (1+(1/k)*log(c/d)); b = (1-Er)/denb; db = (-dEr*denb - (1-Er)*(1/k)*((d*dc-c*dd)/(c*d)))/(denb^2); dy = dEr*x + Er*dx + db*(x+(1/k)*(log(g/d))) + b*(dx+(1/k)*((d*dg-g*dd)/(g*d))); dPr = th(32)*dy + Pbreathe(3); end dP = [0; ((qlovr-qa)/th(3)); (qa-qri)/th(4); dPr; 0; 0; 0; 0];