% A function that solves the 7 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))'; % 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.33; A = 5e5; % Used to be 5e5; Pref = 10+ 10*sin(2*pi*f*t); % Rref = 0.2792; Rref=0.4140; ptop = Pref; 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 = 0.25; %% Setting up BG parameters E = 25; h1 = 0.9732; rho1 = 1.225; R = zeros(N,1); Rv = zeros(N,1); CC = zeros(N,1); I = zeros(N,1); pa_bar = zeros(N,1); alph = zeros(N,1); pbot = 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); I(i) = rho1*C.L(air_ord)/(pi.*r(i).^2); R(i) = C.alpha(air_ord).*r(i).^(-4); Rv(i) = 0.01/CC(i); pa_bar(i) = (10*E1)./sqrt(E1.^2 + (2*pi*f*R(i)).^2); alph(i) = atan((2*pi*f*R(i))./E1); pbot(i) = 10+pa_bar(i)*sin(2*pi*f*t-alph(i)); 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)-pbot(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)) = (uv(numOrd1+i) - u(aa(i)) - (R(aa(i))*vv(aa(i)))/2)./(I(aa(i))/2); vv(bb(i)) = (uv(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(7) - u(5) - (vv(5)*R(5)))/(I(5)); vv(6) = (u(7) - u(6) - (vv(6)*R(6)))/(I(6)); 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) 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*(-ptop*C.vtop - 1/lambda*qhat*C.vbot))); %optimising code gammajplus = C.Cjplus*p + ptop.*C.vtop; gammajminus = C.Cjminus*p + pbot.*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); % % %% Determine what Ptm is (transmural pressure) % % pmid = 0.5*(gammajplus + gammajminus); Ptm = zeros(N,1); for i=1:N air_ord = C.order(i); % Determine the order of the current airway Ptm(i) = pmid(i) - (kappa*Rref/r(i))+ tau(i); if Ptm(i) <= 0 R(i) = sqrt((C.Ri(air_ord).^2)*(1 - Ptm(i)./C.P1(air_ord)).^(-C.n1(air_ord))); elseif Ptm(i) >= 0 R(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)-r'); end RRDOT=RDOT; % size(RRDOT) UVDOT = uv'; % size(UVDOT) VVDOT = vv'; % size(VVDOT) VDSDOT = vds'; % size(VDSDOT) NN = [RRDOT;UVDOT;VVDOT;VDSDOT];