% A function that solves the n coupled airway BG model function [NN] = seven_airway(t,y,order,parm) numBr = 2^order - 1; r = y(1:numBr)'; uv = y(numBr+1: 2*numBr)'; vv = y(2*numBr+1:3*numBr)'; vds = y(3*numBr+1:(end-4))'; % RR = y(end-6:end)'; RR = y(end-3:end)'; M = RR(1); Mp = RR(2); AMp = RR(3); AM = RR(4); % r1 = [y(1) y(2) y(3) y(4) y(5) y(6) y(7)]; % uv1 = [y(8) y(9) y(10) y(11) y(12) y(13) y(14)]; % vv1 = [y(15) y(16) y(17) y(18) y(19) y(20) y(21)]; % vds1 =[y(22) y(23) y(24) y(25)]; %% Construct the constant array as a persistent variable % So that repeated calls to this function doesn't result in repeated % calculations of C persistent C if isempty(C) fprintf('Contrusting the structure array of constants...\n'); C = createConstSymm(order); fprintf('Finished.\n'); end %% Parameters for MoC model rho=2; %% Define the constants f = 0.2; A = 5e5; % Used to be 5e5; % Pref = 10+ 10*sin(2*pi*f*t).*(t<=60)+ 120*sin(2*pi*f*t).*(t>60 && t<62)+ 10*sin(2*pi*f*t).*(t>=62); % Rref = 0.2792; P0 = 10; Pref = 10 + P0*sin(2*pi*f*t); Rref=0.4140; u_top = Pref; % kappa = 0.25*parm(1,:).*(t<=60) + 0.5*parm(1,:).*(t>60 && t<=120) + 0.75*parm(1,:).*(t>120 && t<=180) + (parm(1,:)).*(t>=180); kappa = parm(1,:); % epsilon = 10^(parm(2,:)); % A parameter that controls sliding between Type II and Type I + II epsilon = 0; qhat = -50; N = 2^order - 1; % number of branches for symmetric trees numOrd1 = (N+1)/2; NT = (N-1)/2; E1 = 2.5; %cmH20/ml %% Setting up BG parameters E = 25;%kPa rho1 = 1.225; k11 = 0.0057; k22 = 0.2096; k33 = 0.0904; R = zeros(N,1); R_0 = zeros(N,1); Rv = zeros(N,1); CC = zeros(N,1); I = zeros(N,1); pa_bar = zeros(N,1); alph = zeros(N,1); u_bot = zeros(N,1); u=zeros(N,1); for i=1:N air_ord = C.order(i); %CC(i) = (2*pi*C.L(air_ord).*r(i).^3)/(E*h1); WT(i) = (k11.*(2*r(i)).^2)/4 + (k22.*(2.*r(i)))/2 +k33;%wall thickness CC(i) = (C.L(air_ord)*(2*r(i)).^3)/(4*E*WT(i));% taken from Yoon et al 2020 I(i) = rho1*C.L(air_ord)/(pi.*r(i).^2);%airway Inertance R(i) = C.alpha(air_ord).*r(i).^(-4); %airway resistance Rv(i) = 0.01/CC(i);%air viscoelastic effect pa_bar(i) = (P0*E1)./sqrt(E1.^2 + (2*pi*f*R(i)).^2); alph(i) = atan((2*pi*f*R(i))./E1); u_bot(i) = 10+pa_bar(i)*sin(2*pi*f*t-alph(i));%acinar pressure end %% for a = numOrd1+1:N uv(a) = (vv(a)-vv(2*(a-numOrd1))-vv(2*(a-numOrd1)-1))/CC(a); u(a) = uv(a) +Rv(a)*((2*(a-numOrd1))-vv(2*(a-numOrd1)-1)); end for b = 1:numOrd1 uv(b) = (vv(b)-vds(b))/CC(b); u(b) = uv(b) + Rv(b)*(vv(b)-vds(b)); vds(b) = (uv(b)-u_bot(b)-(R(b)/2)*vds(b))/(I(b)/2); end for d = 1:N-1 a(d) =((2*d) - 1); aa = a(a numBr && bb(i) > numBr % % vv(aa(i)) = (uv(numOrd1+i) - u(aa(i)) - (R(aa(i))*vv(aa(i))))./(I(aa(i))); % % vv(bb(i)) = (uv(numOrd1+i) - u(bb(i)) - (R(bb(i))*vv(bb(i))))./(I(bb(i))); % % end % end % vv(N) = (Pref - u(N) - R(N)*vv(N))/I(N); %% Set up uv1 - uv7 % uv(1) = (vv(1)-vds(1))/CC(1); % u(1) = uv(1)+Rv(1)*(vv(1)-vds(1)); % uv(2) = (vv(2)-vds(2))/CC(2); % u(2) = uv(2)+Rv(2)*(vv(2)-vds(2)); % uv(3) = (vv(3)-vds(3))/CC(3); % u(3) = uv(3)+Rv(3)*(vv(3)-vds(3)); % uv(4) = (vv(4)-vds(4))/CC(4); % u(4) = uv(4)+Rv(4)*(vv(4)-vds(4)); % uv(5) = (vv(5)-vv(2)-vv(1))/CC(5); % u(5) = uv(5)+Rv(5)*(vv(5)-vv(2)-vv(1)); % uv(6) = (vv(6)-vv(4)-vv(3))/CC(6); % u(6) = uv(6)+Rv(6)*(vv(6)-vv(4)-vv(3)); % uv(7) = (vv(7)-vv(6)-vv(5))/CC(5); % u(7) = uv(7)+Rv(7)*(vv(7)-vv(6)-vv(5)); %% Set up v1 - v7 % vv(1) = (u(5) - u(1) - (vv(1)*R(1))/2)/(I(1)/2); % vv(2) = (u(5) - u(2) - (vv(2)*R(2))/2)/(I(2)/2); % vv(3) = (u(6) - u(3) - (vv(3)*R(3))/2)/(I(3)/2); % vv(4) = (u(6) - u(4) - (vv(4)*R(4))/2)/(I(4)/2); for i = 1: (numOrd1/2) vv(aa(i)) = (u(numOrd1+i) - u(aa(i)) - (R(aa(i))*vv(aa(i)))/2)./(I(aa(i))/2); vv(bb(i)) = (u(numOrd1+i) - u(bb(i)) - (R(bb(i))*vv(bb(i)))/2)./(I(bb(i))/2); end for i = (numOrd1/2)+1 : length(aa) vv(aa(i)) = (u(numOrd1+i) - u(aa(i)) - (vv(aa(i))*R(aa(i))))/(I(aa(i))); vv(bb(i)) = (u(numOrd1+i) - u(bb(i)) - (vv(bb(i))*R(bb(i))))/(I(bb(i))); end % vv(7) = (ptop - u(7) - (vv(7)*R(7)))/(I(7)); vv(N) = (Pref - u(N) - R(N)*vv(N))/I(N); %% Set up vds1 - vds4 % vds(1) = (u(1)-pbot(1)-(R(1)/2)*vds(1))/(I(1)/2); % vds(2) = (u(2)-pbot(2)-(R(2)/2)*vds(2))/(I(2)/2); % vds(3) = (u(3)-pbot(3)-(R(3)/2)*vds(3))/(I(3)/2); % vds(4) = (u(4)-pbot(4)-(R(4)/2)*vds(4))/(I(4)/2); %% Determine what R is (airway radius) RDOT=zeros(N,1); %% Find expressions for the pressure and flows Ci = C.Ciplus - C.Ciminus1 - C.Ciminus2; % size(C.alpha.^(-1)) % size(r') Dalphar = diag(C.alpha.^(-1).*r'.^4); % size(Dalphar) % Dalphar1 = diag(R) % return W = Ci*Dalphar*(C.Cjplus - C.Cjminus); % size(W) % return lambda = dot(C.alpha.^(-1).*r'.^4, C.vbot); temp = (W\(Ci*Dalphar*(C.vbot*C.vbot')*Dalphar*C.Cjplus))/lambda; %optimising code Lambda = eye(size(temp)) - temp; p = Lambda\(W\(Ci*Dalphar*(-u_top*C.vtop - 1/lambda*qhat*C.vbot))); %optimising code gammajplus = C.Cjplus*p + u_top.*C.vtop; gammajminus = C.Cjminus*p + u_bot.*C.vbot; deltap = gammajplus - gammajminus; % T = C.T; % The transformation matrix that will map mu_hat to mu % % %% Determine what mu is (parenchymal shear modulus) % mu = zeros(N,1); % for i = 1:numOrd1 % Wanting to loop through only the order 1 airways % Add the Type 1 coupling by putting in the dependence on the % neighbours. if i == 1 mu(i) = abs(A/3*(deltap(i)*r(i)^4 + epsilon*(deltap(numOrd1)*r(numOrd1)^4 + deltap(i+1)*r(i+1)^4)))/2; elseif i == numOrd1 mu(i) = abs(A/3*(deltap(i)*r(i)^4 + epsilon*(deltap(i-1)*r(i-1)^4 + deltap(1)*r(1)^4)))/2; else mu(i) = abs(A/3*(deltap(i)*r(i)^4 + epsilon*(deltap(i-1)*r(i-1)^4 + deltap(i+1)*r(i+1)^4)))/2; end end % % % Use the transformation matrix, T, to create the full mu vector mu = T*mu; % % %% Determine what tau is (parenchymal tethering pressure) tau = 2*mu.*(((Rref - r')/Rref) + 1.5*((Rref - r')/Rref).^2); k2 = 0.5; k5 = 0.5; k3 = 0.4; k4 = 0.1; k7 = 0.01; k1 = 0.55.*(t<=5)+0.06.*(t>5); k6 =k1; % dt=0.5*t; % force=kappa*(rr(2)+rr(5)); force = kappa*(AMp+AM); % pmid = 0.5*(gammajplus + gammajminus); % Ptm = zeros(N,1); Ptm = pmid - (force*Rref./r')+ tau; for i=1:N air_ord = C.order(i); % Determine the order of the current airway % Ptm(i) = pmid(i) - (force*Rref/r(i))+ tau(i); if Ptm(i) <= 0 R_0(i) = sqrt((C.Ri(air_ord).^2)*(1 - Ptm(i)./C.P1(air_ord)).^(-C.n1(air_ord))); elseif Ptm(i) >= 0 R_0(i) = sqrt(C.rimax(air_ord).^2 - (C.rimax(air_ord).^2 - C.Ri(air_ord).^2)*(1 - Ptm(i)/C.P2(air_ord)).^(-C.n2(air_ord))); end RDOT=rho*((R_0)-r'); end F = [-k1*M+k2*Mp+k7*AM; k4*AMp+k1*M-(k2+k3)*Mp; k3*Mp+k6*AM-(k5+k4)*AMp; k5*AMp-(k6+k7)*AM]; RRDOT=RDOT; % size(RRDOT) UVDOT = uv'; % size(UVDOT) VVDOT = vv'; % size(VVDOT) VDSDOT = vds'; % size(VDSDOT) NN = [RRDOT;UVDOT;VVDOT;VDSDOT;F];