OSAModel/OSA_O2_Model_Sim.m at main · Cardiovascular-Modeling-Laboratory/OSAModel · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
% OSA model code for simulations
% Total time for data (used for paper):
    % Simualtionx - 16 minutes
    % SevereOSA - 9.3 minutes
    % Normal - 10 minutes
clear; clc; close all;

% Load physiological parameters
load OSAParameters.mat
Alv_pars = [yO2, PB, PH2O, R, T, Beta_p, k_l];


% Lookup table approximation for oxygenation
Cd_T = (linspace(0,2e-6,8000))';
S_T = SatO2_2(Cd_T,Beta_p);
CT_T = 4*CHb.*S_T + Cd_T;

% Load breathing pattern file
prompt = "Enter file name for breathing pattern"; % Create a .mat file which contains time and alveolar volume data with variables named as below
file = input(prompt,'s');
Input = load(file,'t_vec','Alv_Vol');
C = struct2cell(Input);
t_Data = C{1};
V_Data = C{2};

% Solving parameters
dt = 0.01;
dz = vpci*dt;
k = floor(Lc/dz);
n_dt = floor(tcyclei/dt);
a_dt = floor(tsai/dt);
c_dt = floor(tsyscapi/dt);
prompt2 = "Enter total time for data in minutes";
total_t = input(prompt2); % Total time in seconds
s = floor(total_t*60/dt);

% Initialize variable vectors
t = zeros(1,s);
PAO2 = zeros(1,s);  CAO2 = zeros(1,s);
PpaO2 = zeros(1,s); CpaO2d = zeros(1,s); CpaO2T = zeros(1,s); SpaO2 = zeros(1,s);
PpvO2 = zeros(1,s); CpvO2d = zeros(1,s); CpvO2T = zeros(1,s); SpvO2 = zeros(1,s);
PsaO2 = zeros(1,s); CsaO2d = zeros(1,s); CsaO2T = zeros(1,s); SsaO2 = zeros(1,s);
PsvO2 = zeros(1,s); CsvO2d = zeros(1,s); CsvO2T = zeros(1,s); SsvO2 = zeros(1,s);
PcO2 = zeros(1,s); CcO2 = zeros(1,s);
CpcO2T = zeros(k+1,s); CpcO2d = zeros(k+1,s); SpcO2 = zeros(k+1,s); PpcO2 = zeros(k+1,s);
z = zeros(k+1,1);
CpcO2T1 = zeros(k+1,1); CpcO2d1 = zeros(k+1,1);
V = zeros(1,s);

% Initialize variables (t=0)
PAO2(1) = PAO2i;
CAO2(1) = PAO2(1)*Beta_p;

PpvO2(1) = PAO2i;
CpvO2d(1) = PpvO2(1)*Beta_p;
SpvO2(1) = SatO2_2(CpvO2d(1),Beta_p);
CpvO2T(1) = 4*CHb*SpvO2(1) + CpvO2d(1);

CpaO2T(1:n_dt) = CpvO2T(1)-MRO2/Qs;
CpaO2d(1:n_dt) = interp1(CT_T,Cd_T,CpaO2T(1));
PpaO2(1:n_dt) = CpaO2d(1)/Beta_p;
SpaO2(1:n_dt) = SatO2_2(CpaO2d(1),Beta_p);

PsaO2(1:a_dt+1) = PpvO2(1);
CsaO2d(1:a_dt+1) = PsaO2(1)*Beta_p;
SsaO2(1:a_dt+1) = SatO2_2(CsaO2d(1),Beta_p);
CsaO2T(1:a_dt+1) = 4*CHb*SsaO2(1) + CsaO2d(1);

PsvO2(1:c_dt+1) = PpaO2(1);
CsvO2d(1:c_dt+1) = CpaO2d(1);
SsvO2(1:c_dt+1) = SpaO2(1);
CsvO2T(1:c_dt+1) = CpaO2T(1);

for n=2:k+1
    CpcO2T1(1) = CpaO2T(1);
    CpcO2d1(1) = CpaO2d(1);
    z(n) = z(n-1)+dz;
    CpcO2T1(n) = CpcO2T1(n-1)+(k_l/Qp)*(dz/Lc)*(CAO2(1)-CpcO2d1(n-1));
    CpcO2d1(n) = interp1(CT_T,Cd_T,CpcO2T1(n));
end

CpcO2T(:,1) = CpcO2T1;
CpcO2d(:,1) = CpcO2d1;
SpcO2(:,1) = SatO2_2(CpcO2d(:,1),Beta_p);

CcO2(1) = trapz(z(:,1),CpcO2d(:,1))./Lc;
PcO2(1) = CcO2(1)/Beta_p;

V(1) = 2300;
t(1) = 0;

for i=2:s
            %Time
            t(i) = (i-1)*dt;

            %Alveolar Volume
            V(i) = interp1(t_Data,V_Data,t(i));

            %Calculate arterial oxygenation
            CpaO2T(i+n_dt-1) = CpvO2T(i-1)-MRO2/Qs;
            CpaO2d(i) = interp1(CT_T,Cd_T,CpaO2T(i));
            PpaO2(i) = CpaO2d(i)/Beta_p;
            SpaO2(i) = SatO2_2(CpaO2d(i),Beta_p);

            % Pulmonary capillary oxygenation
            CpcO2T(1,i) = CpaO2T(i);
            CpcO2T(2:end,i) = CpcO2T(1:(end-1),i-1) + (k_l/Qp)*(dz/Lc)*(CAO2(i-1)-CpcO2d(1:(end-1),i-1));
            CpcO2d(:,i) = interp1(CT_T,Cd_T,CpcO2T(:,i));
            PpcO2(:,i) = CpcO2d(:,i)/Beta_p;
            SpcO2(:,i) = SatO2_2(CpcO2d(:,i),Beta_p);

            % Venous oxygenation
            CpvO2T(i) = CpcO2T(end,i);
            CpvO2d(i) = interp1(CT_T,Cd_T,CpvO2T(i));
            PpvO2(i) = CpvO2d(i)/Beta_p;
            SpvO2(i) = SatO2_2(CpvO2d(i),Beta_p);

            % Systemic capillaries input
            CsaO2T(i+a_dt) = CpvO2T(i);
            CsaO2d(i) = interp1(CT_T,Cd_T,CsaO2T(i));
            PsaO2(i) = CsaO2d(i)/Beta_p;
            SsaO2(i) = SatO2_2(CsaO2d(i),Beta_p);

            % Systemic capillaries output
            CsvO2T(i+c_dt) = CsaO2T(i)-MRO2/Qs;
            CsvO2d(i) = interp1(CT_T,Cd_T,CsvO2T(i));
            PsvO2(i) = CsvO2d(i)/Beta_p;
            SsvO2(i) = SatO2_2(CsvO2d(i),Beta_p);

            % Alveolar oxygen partial pressure
            CcO2(i) = trapz(z(:,1),CpcO2d(:,i))./Lc;
            PcO2(i) = CcO2(i)/Beta_p;
            PAO2(i) = Alv_MassBalance(V(i),V(i-1),dt,PcO2(i-1),PAO2(i-1),Alv_pars);
            CAO2(i) = PAO2(i)*Beta_p;
end