-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgetSS.m
129 lines (104 loc) · 4.18 KB
/
getSS.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
function [SS, exitflag, residual] = getSS(IC, params, varargin)
% default settings, varargin is used to change settings
do_figs = false;
do_insulin = true;
do_FF = true;
MKX = 0; MKXslope = 0;
alt_sim = false; % use alternate equations
for i = 1:2:length(varargin)
temp = varargin{i+1};
if strcmp(varargin{i}, 'alt_sim')
alt_sim = temp;
elseif strcmp(varargin{i}, 'do_MKX')
MKX = temp(1);
MKXslope = temp(2);
elseif strcmp(varargin{i}, 'do_insulin')
do_insulin = temp(1);
elseif strcmp(varargin{i}, 'do_FF')
do_FF = temp(1);
elseif strcmp(varargin{i}, 'do_figs')
do_figs = temp(1);
else
disp('WRONG VARARGIN INPUT')
fprintf('What is this varargin input? %s \n', varargin{i})
error('wrong varargin input')
end % if
end %for
% ODE options
tspan = [0 4000];
options = odeset('RelTol',1.0e-6,'AbsTol',1e-9);
%% solve ODE system
%fprintf('solving ODEs \n')
[t,y] = ode15s(@(t,y) kreg_eqns(t,y,params,...
'SS', true, ...
'alt_sim', alt_sim,...
'do_MKX', [MKX, MKXslope],...
'do_insulin', do_insulin,...
'do_FF', do_FF), ...
tspan, IC, options);
if do_figs
%% make figures
fprintf('making figures \n')
% figure specs
lw = 3;
f.xlab = 16; f.ylab = 16; f.title = 18;
cmap = parula(5);
c1 = cmap(3,:);
cgraymap = gray(5);
cgray = cgraymap(3,:);
lwgray = 2; lsgray = '--';
% variables
figure(1)
clf
nrows = 2; ncols = 3;
subplot(nrows,ncols,1)
plot(t,y(:,1),'linewidth',lw,'color',c1)
ylabel('M_{Kgut}', 'fontsize', f.ylab)
xlabel('t', 'fontsize', f.xlab)
title('Gut K', 'fontsize', f.title)
grid on
subplot(nrows,ncols,2)
plot(t,y(:,2),'linewidth',lw,'color',c1)
ylabel('M_{Kplas}', 'fontsize', f.ylab)
xlabel('t', 'fontsize', f.xlab)
title('Plasma K', 'fontsize', f.title)
grid on
subplot(nrows,ncols,3)
plot(t,y(:,3),'linewidth',lw,'color',c1)
ylabel('M_{Kinter}', 'fontsize', f.ylab)
xlabel('t', 'fontsize', f.xlab)
title('Interstitial K', 'fontsize', f.title)
grid on
subplot(nrows,ncols,4)
plot(t,y(:,4),'linewidth',lw,'color',c1)
ylabel('M_{Kmuscle}', 'fontsize', f.ylab)
xlabel('t', 'fontsize', f.xlab)
title('Muscle K', 'fontsize', f.title)
grid on
end
%% take endpoint and make initial guess for fsolve
fprintf('fsolve for SS \n')
opts = optimoptions('fsolve', 'Display', 'off', 'MaxFunEvals', 10000, 'MaxIter', 10000);
IG = y(end,:);
[SS, residual, exitflag, output] = fsolve(@(y) kreg_eqns(0, y, params, ...
'SS', true, ...
'alt_sim', alt_sim,...
'do_MKX', [MKX, MKXslope],...
'do_insulin', do_insulin,...
'do_FF', do_FF),...
IG, opts);
% check between ODE and fsolve result
fracchange = (SS - IG)./IG;
if max(abs(fracchange)) > 0.1
fprintf('WARNING: maximum ODE to fsolve change: %0.4f \n', max(abs(fracchange)))
fprintf('steady states \n')
fprintf(' ODE fsolve \n')
fprintf('M_Kgut %0.4f %0.4f \n', IG(1), SS(1))
fprintf('M_Kplas %0.4f %0.4f \n', IG(2), SS(2))
fprintf('M_Kinter %0.4f %0.4f \n', IG(3), SS(3))
fprintf('M_Kmuscle %0.4f %0.4f \n', IG(4), SS(4))
end
if exitflag < 1
fprintf('***WARNING: exitflag indicates error!***\n')
end
end