-
Notifications
You must be signed in to change notification settings - Fork 0
/
main_DOM_model.m
286 lines (227 loc) · 9.55 KB
/
main_DOM_model.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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
%% Main file for DOM-bacteria-network model
clear variables
close all
%% PARAMETERS AND MODEL SPECIFICATIONS
% Hold the random number generator seed for reproducible results
rng('default')
rng(1)
% Model parameters
numB = 35; % number of bacterial groups
numD = 100; % number of DOC compound groups
K = 10; % half-saturation constant
nexcr = 30; % number of excreted compounds per bacterial group
nsubs = 3; % number of compounds taken-up per bacterial group
eta = 0.2; % growth efficiency of bacteria
beta = 0.14/(1-eta); % fraction of uptake released back to DOC pool via excretion
r_max = 1; % maximum growth rate of bacteria
r_mort = [0.02 0]; % mortality rate of bacteria (second entry is for quadratic mortality term - not used here)
% Generate uptake- and release network
% consumption matrix (bacteria in rows, DOC in columns)
C = get_consumption_matrix(numB, numD, 'universe',...
nsubs, 'non-normalized');
% excretion matrix (bacteria in rows, DOC in columns)
E = get_excretion_matrix(C, nexcr);
% Supply
Stot = 0.08; % total default supply [mmolC/m^3/d]
S = repmat(Stot/numD, 1, numD); % even distribution of total supply across compounds
% Initial condition
sumB0 = 1; % total intital carbon concentration of microbes [mmolC/m^3]
sumD0 = 80; % total intital DOC concentration [mmolC/m^3]
B0 = repmat(sumB0/numB, numB, 1); % all bacterial groups start with the same share of carbon concentration
D0 = repmat(sumD0/numD, numD, 1); % all DOC compound groups start with the same share of carbon concentration
% Solver options
maxstep = 50; % maximum step length
reltol = 1e-06; % relative error tolerance
abstol = 1e-20; % absolute error tolerance
%% Minimal working example: Default simulation
% simulate for ten years
tspan = [0 10*365];
% Default simulation, returns:
% - t, the time vector (t x 1)
% - Bt, the carbon concentration of bacterial groups over time [mmolC/m^3] (t x numB)
% - Dt, the carbon concentration of DOC groups over time [mmolC/m^3] (t x numD)
% - It, the carbon concentration of the inorganic pool over time [mmolC/m^3] (t x 1)
% - At_D, the age of carbon in the DOC compound groups over time [days] (t x numD)
% - At_D, the age of carbon in the bacterial biomass over time [days] (t x numB)
[t, Bt, Dt, It, At_D, At_B] = wrap_ode_DOM_model(tspan, C, E, S, ...
beta, eta, r_mort, r_max, K, B0, D0, 'plot', 'on', 'reltol', reltol, 'maxstep', maxstep);
% Total DOC concentration at the end of the simulation in [mmolC/m^3]
sum(Dt(end,:))
% Total carbon concentration of bacteria at the end of the simulation in [mmolC/m^3]
sum(Bt(end,:))
% Total carbon concentration in the inorganic carbon pool (I) at the end of the simulation in [mmolC/m^3]
It(end)
% Concentration of first DOC compound at after 1 year in [mmolC/m^3]
idx = find(t>365, 1, 'first');
Dt(idx,1)
% Concentration of third bacterial group after 1 year in [mmolC/m^3]
Bt(idx,3)
% Age of second DOC compound group at the end of the simulation in years
% (The DOC is younger than the simulation time due to the supply of DOC)
At_D(end,2)/365
% Age of carbon in biomass of fifth bacterial group at the end of the simulation in years
At_B(end,5)/365
% Age distribution of individual compound groups at the end
figure(); histogram(At_D(end,:)/365)
xlabel('Age at the end of the simulation in years')
% Average concentration-weighted age of DOC in years
meanAge = sum(At_D.*Dt./repmat(sum(Dt,2), 1, numD),2)/365;
figure(); plot(t/365, meanAge)
xlabel('simulation time [years]'), ylabel('mean age of DOC [years]')
%% Same simulation without DOC supply (--> linear aging)
% no supply
S_none = zeros(1, numD);
% simulate for ten years
tspan = [0 10*365];
% Simulation without supply
[t, Bt, Dt, It, At_D, At_B] = wrap_ode_DOM_model(tspan, C, E, S_none, ...
beta, eta, r_mort, r_max, K, B0, D0, 'plot', 'on', 'reltol', reltol, 'maxstep', maxstep);
% Total DOC concentration at the end of the simulation in [mmolC/m^3]
sum(Dt(end,:))
% Total carbon concentration of bacteria at the end of the simulation in [mmolC/m^3]
sum(Bt(end,:))
% Total carbon concentration in the inorganic carbon pool (I) at the end of the simulation in [mmolC/m^3]
It(end)
% Concentration of first DOC compound at after 1 year in [mmolC/m^3]
idx = find(t>365, 1, 'first');
Dt(idx,1)
% Concentration of third bacterial group after 1 year in [mmolC/m^3]
Bt(idx,3)
% Age of second DOC compound group at the end of the simulation in years
% (The DOC is younger than the simulation time due to the supply of DOC)
At_D(end,2)/365
% Age distribution of individual compound groups at the end
figure(); histogram(At_D(end,:)/365)
xlabel('Age at the end of the simulation in years')
% Average concentration-weighted age of DOC in years
meanAge = sum(At_D.*Dt./repmat(sum(Dt,2), 1, numD),2)/365;
figure(); plot(t/365, meanAge)
xlabel('simulation time [years]'), ylabel('mean age of DOC [years]')
%% Working example II: Generate data for Figure 6 ("Toy model set-up")
% Simple toy set-up
numB_simp = 4;
numD_simp = 5;
% Prescribe consumption and excretion matrix
C_simp = [1 0 0 0 1; 0 1 0 0 1; 0 1 1 0 0; 0 0 0 1 1];
E_simp = [0 0.9 0.1 0 0; 0 0 0.9 0.1 0; 0 0 0 0.8 0.2; 0.6 0 0.4 0 0];
% Check matrices
assert(all(all(C_simp+E_simp<=1)), 'microbes should not take up compounds produced by themselves')
% No supply of DOC
S_simp = zeros(1, numD_simp);
% Increase fraction released for illustration purpose (to be able to see
% the re-working of DOC compounds through the network)
beta_simp = 2.2*0.14/(1-eta);
% leave the other parameters at default level
r_mort_simp = r_mort;
K_simp = K;
% solver time span
tspan_simp = [0, 600];
% start with equal amounts of all bacteria
B0_simp = repmat(0.02, 1, numB_simp);
% start with only the first DOC compound present at high concentrations
D0_simp = [20 zeros(1, numD_simp-1)];
% Simulate toy model
[t_simp, Bt_simp, Dt_simp, Ct_simp, At_simp] = wrap_ode_DOM_model(tspan_simp,...
C_simp, E_simp, S_simp, beta_simp, eta, r_mort_simp, r_max, K_simp, ...
B0_simp, D0_simp, 'plot', 'off', ...
'reltol', reltol, 'abstol', abstol, 'maxstep', maxstep);
%% Plot Figure 6 (Toy model)
figure('color', 'white', 'position', [328,294,365,480])
subplot(3,2,1)
imagesc(C_simp)
caxis([0 1])
ax = gca;
ax.XTick = 1:numD_simp;
ax.YTick = 1:numB_simp;
ax.YTickLabel = arrayfun(@(i) sprintf('B_%d', i), 1:numB_simp, 'Unif', false);
ax.XTickLabel = arrayfun(@(i) sprintf('D_%d', i), 1:numD_simp, 'Unif', false);
textStrings = num2str(C_simp(:),'%1.0f');
textStrings = strtrim(cellstr(textStrings));
[x,y] = meshgrid(1:numD_simp, 1:numB_simp);
hStrings = text(x(:),y(:),textStrings(:),...
'HorizontalAlignment','center');
set(hStrings(C_simp(:) > 0), 'FontWeight', 'bold');
set(hStrings(C_simp(:) == 0), 'Color', [.5 .5 .5]);
title('Uptake matrix {\bf{\itU}}')
axU = gca;
axis equal
a = gca;
b = copyobj(gca, gcf);
set(gca, 'Box', 'off', 'YColor', [1 1 1], 'XColor', [1 1 1], ...
'XTickLabel', [], 'YTickLabel', [])
set(b, 'Xcolor', 'k', 'YColor', 'k', 'Box', 'off')
b.Title.String = '';
uistack(b,'down')
subplot(3,2,2)
imagesc(E_simp)
caxis([0 1])
colormap([1 1 1]), axis square
ax = gca;
ax.XTick = 1:numD_simp;
ax.YTick = 1:numB_simp;
ax.YTickLabel = arrayfun(@(i) sprintf('B_%d', i), 1:numB_simp, 'Unif', false);
ax.XTickLabel = arrayfun(@(i) sprintf('D_%d', i), 1:numD_simp, 'Unif', false);
textStrings = num2str(E_simp(:),'%1.1f');
textStrings = strtrim(cellstr(textStrings));
for i = 1:length(textStrings)
if strcmp(textStrings(i), '0.0')
textStrings(i) = {'0'};
elseif strcmp(textStrings(i), '1.0')
textStrings(i) = {'1'};
end
end
[x,y] = meshgrid(1:numD_simp, 1:numB_simp);
hStrings = text(x(:),y(:),textStrings(:),...
'HorizontalAlignment','center');
midValue = mean(get(gca,'CLim'));
set(hStrings(E_simp(:) > 0), 'FontWeight', 'bold');
set(hStrings(E_simp(:) == 0), 'Color', [.5 .5 .5]);
title('Release matrix {\bf{\itR}}')
axR = gca;
axis equal
a = gca;
b = copyobj(gca, gcf);
set(gca, 'Box', 'off', 'YColor', [1 1 1], 'XColor', [1 1 1], ...
'XTickLabel', [], 'YTickLabel', [])
set(b, 'Xcolor', 'k', 'YColor', 'k', 'Box', 'off')
b.Title.String = '';
uistack(b,'down')
%%
subplot(3,2, [3 4])
hold(gca, 'all')
gcols = [1 0 0; 0 0 0; .3 .3 .3; .45 .45 .45; .7 .7 .7; .85 .85 .85];
ps = plot(t_simp, sum(Dt_simp,2), 'color', 'green', 'linewidth', 2);
set(gca, 'ColorOrder', gcols, 'YLim', [-0.5 sum(D0_simp)*1.22])
pi = plot(t_simp, Dt_simp, 'linewidth', 1.7);
uistack(pi(3), 'top')
uistack(pi(5), 'top')
[ld, ldt] = legend([pi; ps], [arrayfun(@(i) sprintf('D_%d', i), 1:numD_simp, 'Unif', false),...
'total']);
yd = ylabel('DOC concentration');
axD = gca;
axD.XTick = [];
axD.YTick = [];
subplot(3,2, [5 6])
hold(gca, 'all')
set(gca, 'ColorOrder', gray(numB_simp+2), 'YLim', [-0.05 3.1])
ps = plot(t_simp, sum(Bt_simp,2), 'color', 'red', 'linewidth', 2);
pi = plot(t_simp, Bt_simp, 'linewidth', 1.7);
[lb, lbt] = legend([pi;ps], [arrayfun(@(i) sprintf('B_%d', i), 1:numB_simp, 'Unif', false),...
'total']);
yb = ylabel('Biomass concentration');
xlabel('Time')
axB = gca;
axB.XTick = [];
axB.YTick = [];
set(findall(gcf,'-property','FontName'), 'FontName','Droid Sans', 'Fontsize', 8)
yb.Position(1) = yd.Position(1);
yd.Position(1) = yd.Position(1);
axB.XLabel.Position(2) = -0.4;
axB.YLabel.Position(1) = -25;
axD.YLabel.Position(1) = -25;
axD.Position(2) = 0.38;
set(ldt(1:6), 'Fontsize', 8, 'FontName', 'Droid Sans')
set(lbt(1:5), 'Fontsize', 8, 'FontName', 'Droid Sans')
set(ld, 'position', [0.7019 0.42 0.2027 0.2173]);
set(lb, 'position', [0.7019 0.1451 0.2027 0.1806]);
set(findall(gcf, '-property', 'fontsize'), 'fontsize', 8)