- Author:
- aram148 <42922407+aram148@users.noreply.github.com>
- Date:
- 2022-07-26 11:49:42+12:00
- Desc:
- updated documentation and renamed model
- Permanent Source URI:
- https://models.fieldml.org/workspace/7e5/rawfile/30e6db7fb82a9bb715b0121cdfc131c19e5cb2ff/Peripheral_matlab_dell_updated/seven_airway.jl
# A function that solves the n coupled airway BG model
function 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:(
return [NN]
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)
@sprintf("Contrusting the structure array of constants...\n")
C = createConstSymm[order]
@sprintf("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 & 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)
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 & 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[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
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]'
end
