- 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_airway2.jl
# A function that solves the n coupled airway BG model
function seven_airway2(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 [U,UVDOT,VVDOT,force,phosphorylation,pmid]
end-4)]'
# RR = y[end-6:end]'
RR = y[end-3:end]'
M = RR[1]; Mp = RR[2]; AMp = RR[3]; AM = RR[4]
## 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("Constructing 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
ptop = 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
E1 = 2.5; #cmH20/ml
## Setting up BG parameters
E = 25
h1 = 0.9732
k11 = 0.0057
k22 = 0.2096
k33 = 0.0904
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)
pmid = 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]);#
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] = (P0*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<N-1];
b[d] = 2*d
bb = b[b<=N-1]
end
## Set up uv1 - uv[n]
## Set up v1 - v[n]
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)
# return
# Dalphar = diag(R)
W = Ci*Dalphar*(C.Cjplus - C.Cjminus)
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)
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)
# stiffness =rr[1]+rr[4]
phosphorylation = AMp+Mp
#
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
#
# drdt[i]=rho*(R[i]-(r[i]));
# V[i]=-((gamma)/(2*pi*Rref))*drdt[i];
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]'
U=u'
end
