- 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/coup.m
%% This model incorporates the Hai-Murphy crossbridge ODEs to the BG-AW-heterogeneity
%% model being developed.
function [t,NN,air,U,UVDOT,VVDOT,force,phosphorylation,summ,pa_bar,R,Ptm]=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 = linspace(0,100,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,:),pa_bar(j,:),R(j,:),Ptm(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';
end
flow_sum=flow_sum./(2^(order-3)*4*sqrt(2));%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=fliplr(air(end,:));
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
