- Author:
- aram148 <a.rampadarath@auckland.ac.nz>
- Date:
- 2022-04-29 10:49:49+12:00
- Desc:
- Adding cross-bridge DI working model
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/7e5/rawfile/231928ae37dcd5eeda4e48113b576d08ef94cfd4/Peripheral_matlab_dell_updated/generalised_model_HM/seven_airway.m
% 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<N-1);
b(d) = 2*d;
bb = b(b<=N-1);
z = [aa bb];
end
% for i = 1: length(aa)
% % if (aa(i)) <=numBr && bb(i) <=numBr
% 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);
% % elseif (aa(i)) > 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];
