- Author:
- aram148 <a.rampadarath@auckland.ac.nz>
- Date:
- 2022-02-10 14:31:03+13:00
- Desc:
- updates to model
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/7e5/rawfile/9eee1d052a54eb1dbf926da85ada39f5e5544793/Peripheral_matlab/generalised_model_HM/mass.m
function [dX] = mass(t,M)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
T = t;
m = M;
% alpha = 1.3e-6; beta = 5.95e-5;
% DO = 3.36; POi = 98;
%%% gas exchange parameters
alpha = 1.3e-6; % O2 solubility by Henry's law(M/mmHg)
beta = 3.08e-5; % CO2 solubility by Hency's law (M/mmHg)
CHb = 0.021; % Hb binding site conc (mol/L of RBC's)
Hct = 0.40; % hematocrit (unitless)
C0 = CHb*Hct; % blood oxygen binding capacity (mol/L)
n = 2.7; % Hill exponent
P50 = 30; % half-max saturation of Hb
DO = 1; % (apparent) O2 diffusion coefficient (ml/s)
Pair = 150;
gamma = 16800; %mmHg * L / mol (convert M to mmHg, for gas)
R_R = 1.17; Q_v = 2.97;
RFt = 0.0154; V0 = 3;
Q_air = Q_v*sin(2*pi*R_R*t);
qi = RFt*Q_air;
A = qi/(2*pi*R_R);
V_alv = A.*(1-cos(2*pi*R_R*t))+ V0;
Pairway1 = M(1); %airway 1 - interacts with atmosphere
Pairway2 = M(2); %airway 2 - big ol' mixing chamber
Malv = M(3:end); %mass of O2 in the alveoli (mol)
PalvO = gamma*Malv./V_alv;
Vairway = 0.6;
f = 0.7;
Vairway1 = Vairway*(1-f);
Vairway2 = Vairway*f;
P_alv = (beta*m)./V_alv;
% m = ((qi/beta).*P_alv - alpha*DO.*(P_alv-POi)).*(Q_air<0) ;
PAO2 = 100; %resting systemic artial O2 partial pressure (mmHg)
PVO2 = 45; %resting systemic venous O2 partial pressure (mmHg)
CAO2 = alpha*PAO2 + C0*(PAO2.^n)./(PAO2.^n + P50.^n); %converted to mol/L
CVO2 = alpha*PVO2 + C0*(PVO2.^n)./(PVO2.^n + P50.^n); %coverted to mol/L
CO = CVO2;%.*ones(10,1); %initial oxygen concentration for capillaries (per capillary segment)
%%
Plookup = 0:1:10000; %look up table
Clookup = alpha*Plookup + C0*(Plookup.^n)./(Plookup.^n + P50.^n);
PO = interp1(Clookup,Plookup,CO); %capillary oxygen tension
%%
dV = zeros(4,1);
if Q_air > 0
dV(1) = (Q_air/Vairway1)*(Pair - Pairway1); %airway1 RHS
dV(2) = (Q_air/Vairway2)*(Pairway1 - Pairway2); %airway2 RHS
dV(3:end) = qi.*Pairway2/gamma - alpha*DO*(PalvO - PO); %Malv RHS
else
dV(1) = -(Q_air/Vairway1)*(Pairway2 - Pairway1); %airway1 RHS
dV(2) = -sum(qi.*(PalvO - Pairway2))./Vairway2; %airway2 RHS
dV(3:end) = qi.*PalvO/gamma - alpha*DO*(PalvO - PO); %Malv RHS
end
dX = dV;
end
