Location: 12 L Platform 1 model codes @ a8a92308e217 / peripheral_airway / Peripheral_matlab / createConstSymm.m

Author:

aram148 <42922407+aram148@users.noreply.github.com>

Date:

2022-07-22 15:47:05+12:00

Desc:

Added documentation for VSM model

Permanent Source URI:

http://models.cellml.org/workspace/6b0/rawfile/a8a92308e217ac5626809237dd90a31240b22834/peripheral_airway/Peripheral_matlab/createConstSymm.m

% A function that will output all of the necessary constants for the
% generalised Type II model in a structure called C. This function can
% create these constants for any arbitrary symmetric tree, with the depth
% of the tree (i.e. its order) as input.

% This function calls on "createConnMatrices.m" to construct the
% connectivity matrices.

% DATE: 18/01/18
% AUTHOR: Austin Ibarra
% VERSION: 1.2

% == LOG CHANGE == %
% v1.1: - Added in the transformation matrix T, that maps mu_hat, a vector
% containing the values of mu for order 1 airways, to mu. This will make
% the calculation more efficient, and the reason we can do this is because
% the structure of our airway doesn't change as time goes on.

% v1.2: - Added in the Horsfield constants from orders 9 to 16.

function C = createConstSymm(order)
%% Create the Connectivity Matrices based on order

[C.Cjplus,C.Cjminus,C.Ciplus,C.Ciminus1,C.Ciminus2] = createConnMatrices(order);

%% Construct the vector that defines the order of each airway

numBr = 2^order - 1;
C.order = zeros(1,numBr);

count = 0;
for i = 1:order
    for j = count + 1:2^i - 1
        C.order(end-count) = order - i + 1;
        count = count + 1;
    end
end

%% Construct the vector vtop and vbot
% vtop is 1 for the top airway position and 0 elsewhere
% vbot is 1 for the terminal airways, and 0 elsewhere

C.vtop = zeros(numBr,1);
C.vtop(end) = 1;

C.vbot = zeros(numBr,1);
for i = 1:2^(order - 1)
    C.vbot(i) = 1;
end

%% v1d1: Construct the transformation T, that maps mu_hat onto mu

C.T = createSymmTMatrix(order);
          
%% Construct the constants based on Horsfield order

% Defined by Horsfield order. Constants taken from Politi et al (2010)
% From order 1 to 16

C.Ri = [0.058 0.065 0.073 0.083 0.096 0.113 0.132 0.156 0.185 0.222 ...
    0.269 0.326 0.395 0.475 0.569 0.686];
C.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];
C.P1 = [15.728 17.342 19.475 22.747 27.140 32.305 39.429 47.104 55.704 65.407 ...
    75.968 88.028 100.441 113.457 130.989 153.036];
C.n1 = [1 1 1 1 1 1 1 1 1 1 ...
    1 1 1 1 1 1];
C.n2 = [7 7 7.185 7.778 8 8 8 8.148 8.741 9.333 ...
    9.926 10 10 10 10 10];
C.L = [1.7 1.8778 2.0556 2.2333 2.4481 2.6852 3.0333 3.3889 3.7444 4.1667 ...
    4.6407 5.0630 5.5111 6.1037 6.7556 7.4667];
C.P2 = C.P1.*C.n2.*(C.Ri.^2 - C.rimax.^2)./(C.n1.*C.Ri.^2);


%% Convert the pressures from Pascals in cmH2O

convRate = 0.01019716213;
C.P1 = C.P1*convRate;
C.P2 = C.P2*convRate;

%% Calculate alpha for each airway

a = 10.19716213; % Conversion factor for kg*(mm)^(-1)*(s)^(-2) to cmH20
mu = 1.9008e-8; % dynamic viscosity kg(mm)^(-1)s^(-1)
C.alpha = zeros(1,numBr)';
for i = 1:numBr
    C.alpha(i) = (8*a*mu*C.L(C.order(i)))/pi;%pi/(8*a*mu*C.L(C.order(i)));    
end
        
end