- 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_Crossbridge/Reduced_DI_tree/seven_airway2.m
% A function that solves the 7 coupled airway BG model
function [U,UVDOT,VVDOT,force,ptop,Rr,Ptm] = seven_airway2(t,y,order,parm)
gamma = 25;
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-7))';
RR = y(end-6:end)';
% RR = y(end-3:end)';
rr(1)=RR(1);
rr(2)=RR(2);
rr(3)=RR(3);
rr(4)=RR(4);
rr(5)=RR(5);
rr(6)=RR(6);
rr(7)=RR(7);
Rref=0.4140;
%% 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('Constructing 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 = 40;
b = 42.5;
t1 = 45;
t2 = 47.5;
P01 = 2;
Pmin = 10;
A = 5e5; % Used to be 5e5;
% Pref = 10+ 10*sin(2*pi*f*t)*exp(-2*pi*f*t);
% Pref = 10+ 10*sin(2*pi*f*t);
Pref = Pmin+P01*sin(2*pi*f*t).*(t<=a)+((10*P01))*sin(2*pi*f*t).*(t>a).* (t<b)+P01*sin(2*pi*f*t).*(t>=b).*(t<=t1)+...
((10*P01))*sin(2*pi*f*t).*(t>t1).*(t<t2)+(P01)*sin(2*pi*f*t).*(t>=t2);
ptop = Pref;
kappa = parm(1,:);
% kappa = 0.1.*parm(1,:).*(t<30) + 0.5*parm(1,:).*(t>=30).*(t<=60) + 0.9*parm(1,:).*(t>60);
% 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;
E1 = 0.25;
%% Setting up BG parameters
E = 25;
h1 = 0.9732;
rho1 = 1.225;
Rr = 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);
for i=1:N
air_ord = C.order(i);
CC(i) = (2*pi*C.L(air_ord).*r(i).^3)/(E*h1);
I(i) = rho1*C.L(air_ord)/(pi.*r(i).^2);
Rr(i) = C.alpha(air_ord).*r(i).^(-4);
Rv(i) = 0.01/CC(i);
pa_bar(i) = ((P01*E1)./sqrt(E1.^2 + (2*pi*f*Rr(i)).^2)).*(t<=a)+((10*P01*E1)./sqrt(E1.^2 + (2*pi*f*Rr(i)).^2)).*(t>a).* (t<b)+...
((P01*E1)./sqrt(E1.^2 + (2*pi*f*Rr(i)).^2)).*(t>=b).*(t<=t1)+((10*P01*E1)./sqrt(E1.^2 + (2*pi*f*Rr(i)).^2)).*(t>t1).* (t<t2)+...
((P01*E1)./sqrt(E1.^2 + (2*pi*f*Rr(i)).^2)).*(t>=t2);%(10*E1)./sqrt(E1.^2 + (2*pi*f*R(i)).^2);
alph(i) = atan((2*pi*f*Rr(i))./E1);
pbot(i) = 10+pa_bar(i)*sin(2*pi*f*t-alph(i));%*exp(-2*pi*f*t);
% pbot(i) = 10+pa_bar(i)*exp((-(CC(i).^-1)*t)/(R(i))-alph(i));
end
% Pref = 10*exp(-((CC(7).^-1)*t)/(R(7)));%+ 10*sin(2*pi*f*t);
% ptop = Pref;
%%
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)-(Rr(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: (numOrd1/2)
vv(aa(i)) = (u(numOrd1+i) - u(aa(i)) - (Rr(aa(i))*vv(aa(i)))/2)./(I(aa(i))/2);
vv(bb(i)) = (u(numOrd1+i) - u(bb(i)) - (Rr(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))*Rr(aa(i))))/(I(aa(i)));
vv(bb(i)) = (u(numOrd1+i) - u(bb(i)) - (vv(bb(i))*Rr(bb(i))))/(I(bb(i)));
end
% vv(7) = (ptop - u(7) - (vv(7)*R(7)))/(I(7));
vv(N) = (Pref - u(N) - Rr(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)
% return
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)
%% Distribution Moment Conditions
k2=0.1;
k1=0.35.*(t<5)+0.06.*(t>=5);
% k1=0.35;
k7=0.005;
% k2=0; k1=k2; k7 = k2;
g1=2.*k7;
g2=20.*g1;
g3=3.*g1;
fp1=0.88;
gp1=0.22;
gp2=4.*(fp1+gp1);
gp3=3.*gp1;
% 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');
% V=-((gamma)/(2*pi*Rref))*RDOT;
end
RDOT=rho*((R')-r');
V=-((gamma)/(2*pi*Rref))*RDOT;
% 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];
%% Distribution Moment stuff
%Mean for cdf
p1=rr(2)/rr(1);
p2=rr(5)/rr(4);
%Standard deviation for cdf
q1=(sqrt((rr(3)/rr(1))-((rr(2)/rr(1)).^(2))));
q2=(sqrt((rr(6)/rr(4))-((rr(5)/rr(4)).^(2))));
%r,phi and I values for 1st PDE M1_lambda
r0=-p1/q1;
r1=(1-p1)/q1;
phi0=0.5*(1+erf((r0-p1)/(q1*sqrt(2))));
phi1=0.5*(1+erf((r1-p1)/(q1*sqrt(2))));
I0=-(exp(-((-p1./q1).^(2))/2))/(sqrt(2*pi));
I1=-(exp(-((((1-p1)./q1)).^(2))/2))/(sqrt(2*pi));
%r,phi and I values for 2nd PDE M2_lambda
r20=-p2/q2;
r21=(1-p2)/q2;
phi20=0.5*(1+erf((r20-p2)/(q2*sqrt(2))));
phi21=0.5*(1+erf((r21-p2)/(q2*sqrt(2))));
I20=-(exp(-((-p2./q2).^(2))/2))/(sqrt(2*pi));
I21=-(exp(-((((1-p2)./q2)).^(2))/2))/(sqrt(2*pi));
%Functions for the rhs of the first PDE M1_lambda
J0=phi0;
% J01=phi1;
% J0inf=phinf;
J10=((p1.*phi0)+(q1.*I0));
J11=(p1.*phi1)+(q1.*I1);
% J12=(p1.*phinf)+(q1.*I2);
J20=((p1.^(2)).*phi0)+((2*p1.*q1).*I0)+((q1.^(2)).*(phi0+(r0*I0)));
J21=((p1.^(2)).*phi1)+((2*p1.*q1).*I1)+((q1.^(2)).*(phi1+(r1*I1)));
% J22=((p1.^(2)))+((2*p1.*q1).*I2)+((q1.^(2)));
J30=(p1.^(3).*phi0)+((3.*p1.^(2).*q1).*I0)+((3.*p1.*q1.^(2)).*((phi0)+(r0*I0)))+((q1.^(3).*(2+r0.^(2)).*I0));
J31=(p1.^(3).*phi1)+((3.*p1.^(2).*q1).*I1)+((3.*p1.*q1.^(2)).*(phi1+(r1*I1)))+(q1.^(3).*(2+(r1.^(2))).*I1);
% J32=(p1.^(3))+((3.*p1.^(2).*q1).*I2)+((3.*p1.*q1.^(2)));
%Functions defined for the RHS of the second PDE M2_lambda
K0=phi20;
% K01=phi21;
% K0inf=phi2inf;
K10=((p2.*phi20)+(q2.*I20));
K11=(p2.*phi21)+(q2.*I21);
% K12=(p2.*phi2inf)+(q2.*I2in);
K20=((p2.^(2)).*phi20)+((2*p2.*q2).*I20)+((q2.^(2)).*(phi20+(r20*I20)));
K21=((p2.^(2)).*phi21)+((2*p2.*q2).*I21)+((q2.^(2)).*(phi21+(r21*I21)));
% K22=((p2.^(2)))+((2*p2.*q2).*I2in)+((q2.^(2)));
K30=(p2.^(3).*phi20)+((3.*p2.^(2).*q2).*I20)+((3.*p2.*q2.^(2)).*((phi20)+(r20*I20)))+((q2.^(3).*(2+r20.^(2)).*I20));
K31=(p2.^(3).*phi21)+((3.*p2.^(2).*q2).*I21)+((3.*p2.*q2.^(2)).*(phi21+(r21*I21)))+(q2.^(3).*(2+(r21.^(2))).*I21);
% K32=(p2.^(3))+((3.*p2.^(2).*q2).*I2in)+((3.*p2.*q2.^(2)));
%Components for the matrix F that will represent each moment,
%M1_lambda and M2_lambda36.883 s
A0= ((fp1*(1-rr(7)))/2)-(fp1*(J11-J10)*rr(1));%-2*(fp1*(K11-K10)*r(4));
A1=((fp1*(1-rr(7)))/3)-(fp1*(J21-J20)*rr(1));%-2*(fp1*(K21-K20)*r(4));
A2=((fp1*(1-rr(7)))/4)-(fp1*(J31-J30)*rr(1));%-2*(fp1*(K31-K30)*r(4));
B0=(gp2*J0)+gp1*(J11-J10)+(gp1+gp3)*(p1-J11);
B1=(gp2*J10)+gp1*(J21-J20)+(gp1+gp3)*((p1.^(2)+q1.^(2))-J21);
B2=(gp2*J20)+gp1*(J31-J30)+(gp1+gp3)*((p1.^(3)+3*p1*q1.^(2))-J31);
C0=k1;
C1=k1*p2;
C2=k1*(p2.^(2)+q2.^(2));
D0=k2;
D1=k2*p1;
D2=k2*(p1.^(2)+q1.^(2));
E0=(g2*K0)+g1*(K11-K10)+(g1+g3)*(p2-K11);
E1=(g2*K10)+g1*(K21-K20)+(g1+g3)*((p2.^(2)+q2.^(2))-K21);
E2=(g2*K20)+g1*(K31-K30)+(g1+g3)*(p2.^(3)+(3*p2*q2.^(2))-K31);
% RDOT= rho*(R-r);
% Putting together the matrix F for the RHS of the system of
% equations
F=[A0-B0*rr(1)+C0*rr(4)-k2*rr(1);A1-B1*rr(1)+C1*rr(4)-k2*rr(2)-V(i)*rr(1);A2-B2*rr(1)+C2*rr(4)-k2*rr(3)-2*V(i)*rr(2);D0*rr(1)-E0*rr(4)-k1*rr(4);D1*rr(1)-E1*rr(4)-k1*rr(5)-V(i)*rr(4);D2*rr(1)-E2*rr(4)-k1*rr(6)-2*V(i)*rr(5);-k1*rr(7)+(1-rr(7))*k2];
RRDOT=RDOT';
% size(RRDOT)
UVDOT = uv';
% size(UVDOT)
VVDOT = vv';
% size(VVDOT)
VDSDOT = vds';
% size(VDSDOT)
% NN = [RRDOT;UVDOT;VVDOT;VDSDOT];
U=u';
