-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathin_circ.m
153 lines (153 loc) · 4.8 KB
/
in_circ.m
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
%filename: in_circ.m (initialization for circ)
T =0.0125 %Duration of heartbeat (minutes)
TS=0.0050 %Duration of systole (minutes)
tauS=0.0025 %CLV time constant during systole (minutes)
tauD=0.0075 %CLV time constant during diastole (minutes)
Rs=17.5 %Systemic resistance (mmHg/(liter/minute))
Rp= 1.79 %Pulmonary resistance (mmHg/(liter/minute))
%Unrealistic valve resistances,
%Chosen small enough to be negligible.
RMi=0.01 %mitral valve resistance (mmHg/(liter/minute))
RAo=0.01 %aortic valve resistance (mmHg/(liter/minute))
RTr=0.01 %tricuspid valve resistance (mmHg/(liter/minute))
RPu=0.01 %pulmonic valve resistance (mmHg/(liter/minute))
%The following values of Csa and Cpa are approximate.
%They will need adjustment to make the systemic
%blood pressure be roughly 120/80 mmHg
%and to make the pulmonary
%blood pressure be roughly 25/8 mmHg.
Csa=0.0006519 %Systemic arterial compliance (liters/mmHg)
Cpa=0.00174 %Pulmonary arterial compliance (liters/mmHg)
Csv=1.75 %Systemic venous compliance (liters/mmHg)
Cpv=0.08 %Pulmonary venous compliance (liters/mmHg)
CLVS=0.00003 %Min (systolic) value of CLV (liters/mmHg)
CLVD=0.0146 %Max (diastolic) value of CLV (liters/mmHg)
CRVS=0.0002 %Min (systolic) value of CRV (liters/mmHg)
CRVD=0.0365 %Max (diastolic) value of CRV (liters/mmHg)
Vsad=0.825 %Systemic arterial volume at P=0 (liters)
Vpad=0.0382 %Pulmonary arterial volume at P=0 (liters)
Vsvd=0 %Systemic venous volume at P=0 (liters)
Vpvd=0 %Pulmonary venous volume at P=0 (liters)
VLVd=0.027 %Left ventricular volume at P=0 (liters)
VRVd=0.027 %Right ventricular volume at P=0 (liters)
dt=0.01*T %Time step duration (minutes)
%This choice implies 100 timesteps per cardiac cycle.
klokmax=30*T/dt %Total number of timesteps
%This choice implies simulation of 15 cardiac cycles.
%Assign an index to each compliance vessel
%of the model circulation:
iLV=1
isa=2
isv=3
iRV=4
ipa=5
ipv=6
N=6
%Enter parameters and initial values
%into correct slots in arrays.
%Note that the code that follows is independent
%of the specific numbering scheme chosen above.
%Compliance vector:
C=zeros(N,1);
%This makes C a column vector of length N.
C(iLV)=CV_now(0,CLVS,CLVD); %initial value
C(isa)=Csa;
C(isv)=Csv;
C(iRV)=CV_now(0,CRVS,CRVD); %initial value
C(ipa)=Cpa;
C(ipv)=Cpv;
C %This writes the result on the screen.
%Pressure vector (initial values) at end of diastole:
P=zeros(N,1);
%This makes P a column vector of length N.
P(iLV)= 5;
P(isa)=80;
P(isv)= 2;
P(iRV)= 2;
P(ipa)= 8;
P(ipv)= 5;
P %This writes the result on the screen.
%Vector of dead volumes (volume at zero pressure);
%Note: Vd is only needed for output purposes.
%It drops out of the equations we solve for P,
%but we need it if we want to output (e.g., plot)
%the volume of any compliance vessel.
Vd=zeros(N,1);
%This makes Vd a column vector of length N.
Vd(iLV)=VLVd;
Vd(isa)=Vsad;
Vd(isv)=Vsvd;
Vd(iRV)=VRVd;
Vd(ipa)=Vpad;
Vd(ipv)=Vpvd;
Vd
%This writes the results on the screen.
%Conductance matrix:
G=zeros(N,N);
%This makes G an NxN matrix filled with zeros.
%Any element of G that is not explicitly
%made nonzero remains zero,
%thus modeling an infinite resistance connection,
%that is, no connection at all.
G(iLV,isa)=1/RAo; %But G(isa,iLV)=0 (no leak)
G(isa,isv)=1/Rs; %no valve
G(isv,isa)=1/Rs; %no valve
G(isv,iRV)=1/RTr; %But G(iRV,isv)=0; (no leak)
G(iRV,ipa)=1/RPu; %But G(ipa,iRV)=0; (no leak)
G(ipa,ipv)=1/Rp; %no valve
G(ipv,ipa)=1/Rp; %no valve
G(ipv,iLV)=1/RMi; %But G(iLV,ipv)=0; (no leak)
G %This writes the result on the screen.
%Matrix of initial valve states:
S=zeros(N,N)
%This makes S an NxN matrix filled with zeros
%(and writes it on the screen).
%Start with all valves closed.
%Valves will adjust to pressures
%during first time step.
%Initialize arrays to store data for plotting:
t_plot=zeros(1,klokmax);
C_plot=zeros(N,klokmax);
P_plot=zeros(N,klokmax);
%Other variables that we might want to plot
%can be found from these.
%For self-checking in P_new, set CHECK=1.
%To skip self-checking set CHECK=0.
%(should be much faster with CHECK=0)
CHECK=0
%Initialize flow computation (for output purposes only)
%assign an index to each flow of interest:
jAo=1;
js =2;
jTr=3;
jPu=4;
jp =5;
jMi=6;
Nflows=6;
%note index of upstream and downstream chamber
%for each flow:
iU=zeros(Nflows,1);
iD=zeros(Nflows,1);
iU(jAo)=iLV;
iD(jAo)=isa;
iU(js )=isa;
iD(js )=isv;
iU(jTr)=isv;
iD(jTr)=iRV;
iU(jPu)=iRV;
iD(jPu)=ipa;
iU(jp )=ipa;
iD(jp )=ipv;
iU(jMi)=ipv;
iD(jMi)=iLV;
%extract the conductances from the matrix G:
Gf=zeros(Nflows,1);
Gr=zeros(Nflows,1);
for j=1:Nflows
Gf(j)=G(iU(j),iD(j)); %forward conductance
Gr(j)=G(iD(j),iU(j)); %reverse conductance
end
%create arrays to store current pressure differences
%and history over time of the net flows:
Pdiff=zeros(Nflows,1);
Q_plot=zeros(Nflows,klokmax);