- 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/coup.jl
using PyPlot
## This model incorporates the Hai-Murphy crossbridge ODEs to the BG-AW-heterogeneity
## model being developed.
function coup(order,kappa)#,rinitial)
tic()
numBr = 2^order - 1; # number of airways
numOrd1 = (numBr+1)/2
# initial values for airways [mm]
rimax = [0.296 0.318 0.337 0.358 0.384 0.414 0.445 0.484 0.539 0.608 ...
0.692 0.793 0.913 1.052 1.203 1.374]
rinitial=[]
for k=1:order
rinitial=[rinitial rimax[k]*ones(2^(order-k),1)]'
end
rinitial=rinitial+0.01*rand(1,length(rinitial))'; # peturbation of rinitial
M = 1; Mp = 0; AMp = 0; AM = 0
RR0= [M Mp AMp AM]'
r[:,1]=rinitial; #airways, order 1 : order n
uv[:,1] = zeros(numBr,1); # pressure in airway
vv[:,1] = zeros(numBr,1); # airflow within the airway
vds[:,1] = zeros((numBr+1)/2,1); # airflow entering alveoli
RR[:,1] = RR0
y0 = [r[:,1] uv[:,1] vv[:,1] vds[:,1] RR[:,1]]'; #ICs vector
#Stepwise kappa value
tspan = range(0,100,length=1000)
parm = kappa
# Set numerical accuracy options for ODE solver
options = odeset("RelTol', 1e-08, 'AbsTol', 1e-08, 'MaxStep", 0.1)
[t,NN]=ode15s(@(t,y) seven_airway[t,y,order,parm], tspan, y0, options)
U = zeros(length(NN),length(r)); # pressure including viscous damping
UVDOT = zeros(length(NN),length(r))
VVDOT = zeros(length(NN),length(r)); #airflow
force = zeros(length(NN),length(r))
phosphorylation = zeros(length(NN),length(r))
pmid = zeros(length(NN),length(r))
for j = 1:length(NN)
[U[j,:],UVDOT[j,:],VVDOT[j,:],force[j,:],phosphorylation[j,:],pmid[j,:]] = seven_airway2[t[j],NN[j,:],order,parm]
end
toc()
# CM=jet[length(rinitial)]
r_new = NN[:,1:numBr]
air = zeros(length(NN),length(r))
for i=1:length(rinitial)
for j = 1:length(t)
if r_new[j,i] <= 0.01
air[j,i]=0
else()
air[j,i]=1
end
end
end
for i = 1:(numOrd1)
# figure(1)
# subplot(order,order,i)
# plot(t,r_new[:,i])
# title(sprintf("Radius of airway [#d]",i))
# #
# # figure(2)
# # subplot(order,10,i)
# # # M[i] = findpeaks[r_new[:,i]]
# # DataInv[:,i] = 1.01*max(r_new[:,i]) - r_new[:,i]
# # [Minima_P(i},MinIIdx{i)] = findpeaks[DataInv[:,i]]
# # abc = r_new[:,i];
# # Minima_P[i] = abc[MinIIdx[i]];
# # plot(t[MinIIdx(i}],Minima_P{i),"-o")
#
figure(2)
subplot(order,order,i)
plot(t,VVDOT[:,i])
title(sprintf("Flow in airway [#d]",i))
# # #
# figure(3)
# subplot(order,order,i)
# plot(t,U[:,i])
# title(sprintf("Pressure in airway [#d]",i))
#
# figure(5)
# subplot(order,order,i)
# plot(t,pmid[:,i])
end
# #
# #
# legendStr = cell(1,length(rinitial))
# for i = 1:length(legendStr)
# legendStr[i] = sprintf("r_[#d]",i)
#
# end
# legend(legendStr,"Location','BestOutside")
# ## Measure for state of heterogeniety
x = cell(1,numOrd1)
for i = 1:length(x)
x[i] = air[:,i]
end
summ = zeros(1,length(x[1]))
for i = 1:length(x)
switch i
case 1
neighbour_diff = sqrt((x(:,i) - x(:,end))^2 + (x(:,i) - x(:,i+1))^2)
case length(x)
neighbour_diff = sqrt((x(:,i) - x(:,i-1))^2 + (x(:,i) - x(:,1))^2)
otherwise()
neighbour_diff = sqrt((x(:,i) - x(:,i-1))^2 + (x(:,i) - x(:,i+1))^2)
end
summ = summ + neighbour_diff'
end
summ=summ/(2^(order-3)*4*sqrt(2))
y=cell(1,numOrd1)
for k=1:length(y)
y[k]=VVDOT[:,k]
end
flow_sum = zeros(1,length(y[1]))
for ii = 1:length(y)
switch ii
case 1
neighbour_diff1 = sqrt((y(:,ii) - y(:,end))^2 + (y(:,ii) - y(:,ii+1))^2)
case length(y)
neighbour_diff1 = sqrt((y(:,ii) - y(:,ii-1))^2 + (y(:,ii) - y(:,1))^2)
otherwise()
neighbour_diff1 = sqrt((y(:,ii) - y(:,ii-1))^2 + (y(:,ii) - y(:,ii+1))^2)
end
flow_sum = flow_sum + neighbour_diff1'
return [t,NN,air,U,UVDOT,VVDOT,force,phosphorylation,summ,flow_sum]
end
flow_sum=flow_sum./numOrd1
# ## graphing airways
node1 = zeros(1,numBr)
node2 = zeros(1,numBr)
# e.g. for an order 4: node1 = [1 2 2 3 3 4 4 5 5 6 6 7 7 8 8]
node1[1] = 1
for i = 1:(numBr-1)/2
node1[2*i:2*i+1] = [i+1 i+1]
end
# e.g. for an order 4: node2 = [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]
for i = 1:numBr
node2[i] = i+1
end
ind = zeros(1,numBr)
for i = 1:numBr
ind[i] = i
end
# figure()
A=reverse(air[end,:], dims = 2)
figure(4)
G = graph[node1,node2,A]
h=plot(G,"LineWidth",2)
h.EdgeCData=A
labelnode[h,[node1 node2],'']
labeledge[h,ind,G.Edges.Weight[ind]]
colormap jet
# for i = 1:8
# subplot(4,2,i)
# plot(time,VVDOT[:,i])
# title(sprintf("Flow of airway [#d]",i))
# end
end
