-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathDEMENT.0.7.2.R
1085 lines (968 loc) · 50.9 KB
/
DEMENT.0.7.2.R
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
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
###############################################################################
sum.grid <- function(x,rn,grid.size) {
y <- NULL
for (i in 1:dim(x)[2]) y <- cbind(y,rowSums(matrix(x[,i],ncol=grid.size)))
dimnames(y) <- list(rn,colnames(x))
y
}
###############################################################################
expand <- function(x,grid.size) {
matrix(rep(t(x),grid.size),ncol=dim(x)[2],byrow=T,dimnames=list(rep(rownames(x),grid.size),colnames(x)))
}
###################################### Plot grid concentrations; grid is a vector with x changing fastest
plot.grid <- function(grid,x,zlim=NULL,title="") {
grid.mat <- matrix(grid,byrow=T,ncol=x)
if(!is.null(zlim)) zlim<-c(0,zlim)
myImagePlot(grid.mat,zlim=zlim,title=title)
}
# ----- Define a function for plotting a matrix ----- #
myImagePlot <- function(x, ...){
min <- min(x)
max <- max(x)
yLabels <- rownames(x)
xLabels <- colnames(x)
title <-c()
# check for additional function arguments
if( length(list(...)) ){
Lst <- list(...)
if( !is.null(Lst$zlim) ){
min <- Lst$zlim[1]
max <- Lst$zlim[2]
}
if( !is.null(Lst$yLabels) ){
yLabels <- c(Lst$yLabels)
}
if( !is.null(Lst$xLabels) ){
xLabels <- c(Lst$xLabels)
}
if( !is.null(Lst$title) ){
title <- Lst$title
}
}
# check for null values
if( is.null(xLabels) ){
xLabels <- c(1:ncol(x))
}
if( is.null(yLabels) ){
yLabels <- c(1:nrow(x))
}
# Red and green range from 0 to 1 while Blue ranges from 1 to 0
ColorRamp <- rgb( seq(0,1,length=256), # Red
seq(0,1,length=256), # Green
seq(1,0,length=256)) # Blue
ColorLevels <- seq(min, max, length=length(ColorRamp))
# Reverse Y axis
reverse <- nrow(x) : 1
yLabels <- yLabels[reverse]
x <- x[reverse,]
# Data Map
image(1:length(xLabels), 1:length(yLabels), t(x), col=ColorRamp, xlab="",
ylab="", axes=FALSE, zlim=c(min,max))
if( !is.null(title) ){
title(main=title)
}
#rect(0.5,0.5,length(xLabels)+0.5,length(yLabels)+0.5)
}
############################################################################
MakePlots <- function(out) {
par(cex.axis=2,cex.lab=2.5,cex.main=3,font.lab=2,font.axis=2,lwd=2,mar=c(5,7,4,2)+0.1,tcl=0.4,las=1,mgp=c(3,0.75,0),mfrow=c(2,5))
RGB.palette <- colorRampPalette(c("red","yellow","blue"),space="rgb")
n_taxa <- dim(out$MicrobesSeries)[2]
colors.RGB <- RGB.palette(out$params["Enz_per_taxon_max",]+1)
plot(0:out$end_time,out$Mic_Sum[,"C"],ylim=c(0,max(out$MicrobesSeries/out$grid.size,na.rm=T)),ylab="",main="Microbial C",xlab="Day",type="n",axes=F,frame.plot=T)
for (i in 1:n_taxa){
lines(0:out$end_time,out$MicrobesSeries[,i]/out$grid.size,col=colors.RGB[(rowSums(out$EnzGenes)+1)[i]],lwd=2)
}
axis(1,lwd=2); axis(2,lwd=2)
n_subs <- length(dimnames(out$SubstratesSeries)[2][[1]])
color.sub <- c("brown","darkblue","goldenrod","salmon","orange","pink","gray","purple","blue","turquoise","forestgreen","green")
plot(0:out$end_time,out$Substrate_Sum,ylim=c(0,max(out$SubstratesSeries/out$grid.size,na.rm=T)),main=expression(bold(paste("Substrate C (mg ",cm^-3,")"))),ylab=NA,xlab="Day",type="n",axes=F,frame.plot=T)
axis(1,lwd=0,lwd.ticks=2); axis(2,lwd=0,lwd.ticks=2)
for (i in 1:n_subs){
lines(0:out$end_time,out$SubstratesSeries[,i]/out$grid.size,col=color.sub[i],lwd=2)
}
legend("topright",legend=c("Dead microbe","Inactive enzyme","Cellulose","Hemicellulose","Starch","Chitin","Lignin","Protein 1","Protein 2","Protein 3","Phospholipid","Nucleic acid"),lty=1,col=color.sub,box.lty=0,cex=1.1)
plot(0:out$end_time,out$Mic_Sum[,"C"]/out$Mic_Sum[,"P"],ylim=c(0,max(out$MicrobesPSeries,na.rm=T)),ylab=NA,main="Microbial C:P",xlab="Day",type="l",axes=F,frame.plot=T)
for (i in 1:n_taxa){
lines(0:out$end_time,out$MicrobesPSeries[,i],col=colors.RGB[(rowSums(out$EnzGenes)+1)[i]])
}
axis(1,lwd=2); axis(2,lwd=2)
plot(0:out$end_time,out$Mic_Sum[,"C"]/out$Mic_Sum[,"N"],ylim=c(0,max(out$MicrobesNSeries,na.rm=T)),ylab=NA,main="Microbial C:N",xlab="Day",type="l",axes=F,frame.plot=T)
for (i in 1:n_taxa){
lines(0:out$end_time,out$MicrobesNSeries[,i],col=colors.RGB[(rowSums(out$EnzGenes)+1)[i]])
}
axis(1,lwd=2); axis(2,lwd=2)
plot(0:out$end_time,out$Mic_Sum[,"N"]/out$Mic_Sum[,"P"],ylim=c(0,max(out$MicrobesNPSeries,na.rm=T)),ylab=NA,main="Microbial N:P",xlab="Day",type="l",axes=F,frame.plot=T)
for (i in 1:n_taxa){
lines(0:out$end_time,out$MicrobesNPSeries[,i],col=colors.RGB[(rowSums(out$EnzGenes)+1)[i]])
}
axis(1,lwd=2); axis(2,lwd=2)
# Color code enzymes according to the main substrate they degrade
enz.potential <- out$ReqEnz[[1]]*out$Vmax
colors <- color.sub[apply(enz.potential,2,which.max)]
plot(0:out$end_time,out$Enzyme_Sum,ylim=c(0,max(out$EnzymesSeries/out$grid.size,na.rm=T)),ylab=NA,main="Enzyme",xlab="Day",type="n",axes=F,frame.plot=T)
for (i in 1:length(colors)){
lines(0:out$end_time,out$EnzymesSeries[,i]/out$grid.size,col=colors[i],lwd=2)
}
axis(1,lwd=2); axis(2,lwd=2)
plot(0:out$end_time,out$NH4Series,ylim=c(0,max(cbind(out$NH4Series/out$grid.size,out$PO4Series/out$grid.size))),
ylab=NA,main="NH4 | PO4",xlab="Day",type="n",axes=F,frame.plot=T)
lines(0:out$end_time,out$NH4Series/out$grid.size,col="orange")
lines(0:out$end_time,out$PO4Series/out$grid.size,col="purple")
axis(1,lwd=2); axis(2,lwd=2)
n_monomers <- length(dimnames(out$MonomersSeries)[2][[1]])
colors <- c("black","black",color.sub)
plot(0:out$end_time,out$Monomer_Sum,ylim=c(0,max(out$MonomersSeries/out$grid.size,na.rm=T)),ylab=NA,main="Monomer C",xlab="Day",type="n",axes=F,frame.plot=T)
for (i in 1:n_monomers){
lines(0:out$end_time,out$MonomersSeries[,i]/out$grid.size,col=colors[i])
}
axis(1,lwd=2); axis(2,lwd=2)
plot(0:out$end_time,out$RespSeries/out$grid.size,ylim=c(0,max(out$RespSeries/out$grid.size)),ylab=NA,main="Respiration",xlab="Day",type="l",axes=F,frame.plot=T)
axis(1,lwd=2); axis(2,lwd=2)
index <- is.finite(log10(apply(out$MicrobesSeries,2,max)))
plot(rowSums(out$EnzGenes)[index],log10(apply(out$MicrobesSeries,2,max)/out$grid.size)[index],main=expression(bold(paste(Log[10]," max density"))),
ylab=NA,xlab="Number of enzyme genes",pch=19,cex=2,axes=F,frame.plot=T,#ylim=c(-2,1),
col=colors.RGB[1+as.vector(rowSums(out$EnzGenes)[index])])
axis(1,lwd=2); axis(2,lwd=2)
abline(lm(log10(apply(out$MicrobesSeries,2,max)/out$grid.size)[index]~rowSums(out$EnzGenes)[index]),lwd=1.75)
} # End of MakePlots() ########################################################
###############################################################################
RunPulse <- function(
params,
timestamp,
rng.seed,
grid.size,
Microbes,
Substrates,
SubInput,
Enzymes,
Monomers,
MonInput,
MonomersProduced,
ReprodNew,
Colonization.reset,
Ea,
Vmax0,
Km0,
ReqEnz,
EnzGenes,
EnzProdInduce,
EnzProdConstit,
UptakeGenes,
UptakeGenesForEnz,
Uptake_ReqEnz,
EnzAttrib,
Uptake_Ea,
Uptake_Vmax0,
Uptake_Km0,
CUE.ref,
OptimalRatios,
MinRatios,
RangeRatios,
fb,
Temp,
Psi,
Tolerance
) {
# Initialize names and indicies
n_taxa <- params["n_taxa",]
n_substrates <- params["n_substrates",]
n_enzymes <- params["n_enzymes",]
Mon.names <- c("NH4","PO4",rownames(Substrates[1:n_substrates,]))
n_monomers <- length(Mon.names)
Mon.names.G <- rownames(Monomers)
Enz.names <- colnames(Vmax0)
Sub.names <- rownames(Substrates[1:n_substrates,])
Mic.names <- rownames(Microbes[1:n_taxa,])
is.NH4 <- which(rownames(Monomers)=="NH4")
is.PO4 <- which(rownames(Monomers)=="PO4")
org <- which(rownames(Monomers) != c("NH4","PO4"))
mineral <- which(rownames(Monomers) == c("NH4","PO4"))
is.deadMic <- which(rownames(Substrates)=="DeadMic")
is.deadEnz <- which(rownames(Substrates)=="DeadEnz")
is.cellulose <- which(rownames(Substrates)=="Cellulose")
is.lignin <- which(rownames(Substrates)=="Lignin")
MonomerRatios <- Monomers
MonomerRatios[,] <- 0
MonomerRatios[is.NH4,"N"] <- 1
MonomerRatios[is.PO4,"P"] <- 1
Tref <- 293
# Create matrices to be filled in
tev <- matrix(rep(0,n_substrates*grid.size*n_enzymes),ncol=n_enzymes)
Max_Uptake1 <- matrix(rep(0,grid.size*n_taxa*n_monomers),nrow=n_taxa)
Max_Uptake2 <- matrix(Max_Uptake1,ncol=n_taxa,dimnames=list(Mon.names.G,Mic.names))
TUC.mat <- TUN.mat <- TUP.mat <- matrix(rep(0,grid.size*n_taxa*n_monomers),nrow=n_monomers)
EP.mat <- matrix(rep(0,grid.size*n_taxa*n_enzymes),nrow=n_taxa,dimnames=list(Mic.names,rownames(Enzymes)))
Death.mat <- matrix(rep(0,grid.size*n_taxa*3),nrow=n_taxa)
DeadEnz.mat <- matrix(rep(0,grid.size*n_enzymes*3),nrow=n_enzymes)
# Create index for matrix manipulation of taxon uptake
# index1 is the shape of the starting matrix
TU.index1 <- matrix(1:(grid.size*n_taxa*n_monomers),ncol=n_taxa)
TU.index <- matrix(TU.index1,nrow=n_monomers)
for(i in 1:grid.size) {
rows <- ((i-1)*n_monomers+1):(i*n_monomers)
cols <- ((i-1)*n_taxa+1):(i*n_taxa)
TU.index[,cols] <- TU.index1[rows,]
}
# Create index for matrix manipulation of taxon enzyme production
EP.index1 <- matrix(1:(grid.size*n_taxa*n_enzymes),ncol=n_enzymes)
EP.index <- matrix(EP.index1,nrow=n_taxa)
for(i in 1:grid.size) {
rows <- ((i-1)*n_taxa+1):(i*n_taxa)
cols <- ((i-1)*n_enzymes+1):(i*n_enzymes)
EP.index[,cols] <- EP.index1[rows,]
}
# Index to manipulate Vmax values
vm.index1 <- matrix(1:(grid.size*n_substrates*n_enzymes),ncol=n_substrates)
vm.index <- matrix(vm.index1,ncol=n_enzymes)
for(i in 1:grid.size) {
vm.index[((i-1)*n_substrates+1):(i*n_substrates),] <- t(vm.index1[((i-1)*n_enzymes+1):(i*n_enzymes),])
}
# Set grid sums to initial values
Monomers.grid <- sum.grid(Monomers,Mon.names,grid.size)
Enzymes.grid <- sum.grid(Enzymes,Enz.names,grid.size)
Substrates.grid <- sum.grid(Substrates,Sub.names,grid.size)
Microbes.grid <- sum.grid(Microbes,Mic.names,grid.size)
n_fungi <- sum(fb)/grid.size
Fungi <- which(fb[1:n_taxa]==1)
is.fungi <- fb==1
Fungi.count.zero <- Fungi.count <- rep(0,n_taxa)
Fungi.vec <- Fungi.vec.zero <- matrix(rep(0,n_taxa*grid.size),ncol=1,dimnames=list(rownames(Microbes),"Count"))
Fungi.count[Fungi] <- Microbes.grid[Fungi,"C"]/(0.5*params["max_size_f",])
z <- rep(1:n_taxa,grid.size)
kill.reset <- rep(0,grid.size*n_taxa)
# Initialize the time series of data to hold
RespSeries <- 0
Cum_Leaching_N <- 0
Cum_Leaching_P <- 0
EnzymesSeries <- t(Enzymes.grid[,"C"])
SubstratesSeries <- Substrates.grid[,"C"]
MonomersSeries <- Monomers.grid[,"C"]
NH4Series <- Monomers.grid["NH4","N"]
PO4Series <- Monomers.grid["PO4","P"]
MicrobesSeries <- Microbes.grid[,"C"]
MicrobesNSeries <- Microbes.grid[,"C"]/Microbes.grid[,"N"]
MicrobesPSeries <- Microbes.grid[,"C"]/Microbes.grid[,"P"]
MicrobesNPSeries <- Microbes.grid[,"N"]/Microbes.grid[,"P"]
# Sum microbe, substrate, monomer, enzyme pools
Mic_Sum <- Cum_Microbe <- colSums(Microbes.grid)
Substrate_Sum <- sum(Substrates.grid[,"C"])
Cum_Substrate <- colSums(Substrates.grid)
Monomer_Sum <- sum(Monomers.grid[,"C"])
Cum_Monomer <- colSums(Monomers.grid)
Enzyme_Sum <- Cum_Enzyme_C <- sum(Enzymes.grid)
Cum_Enzyme_N <- sum(Enzymes.grid*EnzAttrib[,"N_cost"])
Cum_Enzyme_P <- sum(Enzymes.grid*EnzAttrib[,"P_cost"])
# Relative enzyme carbon cost for each enzyme gene
Induce_Enzyme_C <- t(EnzProdInduce)*EnzAttrib[,"C_cost"]
Constit_Enzyme_C <- t(EnzProdConstit)*EnzAttrib[,"C_cost"]
# Relative uptake enzyme carbon cost for each enzyme gene
Uptake_Cost <- t(UptakeGenes)
rownames(Uptake_Cost) <- sprintf("%s%03d","Upt",1:params["n_uptake",])
# Expand variables to the size of the grid
Induce_Enzyme_C <- expand(t(Induce_Enzyme_C),grid.size)
Constit_Enzyme_C <- expand(t(Constit_Enzyme_C),grid.size)
Uptake_Cost <- expand(t(Uptake_Cost),grid.size)
tUptakeGenes <- expand(t(UptakeGenes),grid.size)
UptakeGenes <- expand(UptakeGenes,grid.size)
# Create variable to hold grid at maximum enzyme biomass
grid.max <- NULL
for(i_t in 1:params["end_time",]) { # Loop through time
# Mortality rate increases linearly with -Psi times slope beta
# which can differ for bacteria vs fungi; tolerance of 1 eliminates
# dependence on Psi. Fungi and bacteria can also differ in baseline
# mortality according to Death_Ratio
r.death <- (1-fb)*params["Death_Rate",]*(1 - (params["beta.bac",]*Psi[i_t])*(1 - Tolerance)) +
fb*params["Death_Rate",]*((1 - (params["beta.fungi",]*Psi[i_t])*(1 - Tolerance))*params["Death_Ratio",])
# Calculate Vmax values with Arrhenius equation
# Gas constant = 0.008314 kJ/(mol K)
# Arrhenius equation for Vmax multiplied by exponential decay for Psi sensitivity
Vmax <- Vmax0*exp((-Ea/0.008314)*(1/(Temp[i_t]+273)-1/Tref))*exp(params["Psi.slope.Vmax",]*Psi[i_t])
Uptake_Vmax <- Uptake_Vmax0*exp((-Uptake_Ea/0.008314)*(1/(Temp[i_t]+273)-1/Tref))*exp(params["Psi.slope.uptake",]*Psi[i_t])
# Expand to the size of the grid after calculating
etVmax <- expand(t(Vmax),grid.size)
Uptake_Vmax <- expand(Uptake_Vmax,grid.size)
# Recalculate Km values with Arrhenius equation
Km <- Km0*exp((-params["Km_Ea",]/0.008314)*(1/(Temp[i_t]+273)-1/Tref))
Uptake_Km <- Uptake_Km0*exp((-params["Km_Ea",]/0.008314)*(1/(Temp[i_t]+273)-1/Tref))
# Carbon use efficiency dependent on temperature and number of enzyme genes
CUE <- CUE.ref + (Temp[i_t] - (Tref-273))*params["CUE_temp",]
# Reset the reproduction matrix
Reprod <- ReprodNew
# Reset the colonization matrix
Colonization <- Colonization.reset
# Fungal translocation: calculate average biomass within fungal taxa
Mean.fungi <- Microbes.grid/as.vector(Fungi.count)
Mean.fungi[!is.finite(Mean.fungi)] <- 0
# Expand the fungal average across the grid
eMF <- expand(Mean.fungi,grid.size)
# Reset the fungal count to zero
Fungi.count <- Fungi.count.zero
# Select the daughter cells that are fungi versus bacteria
daughters.b <- which(Reprod[,"C"]>0 & fb==0)
daughters.f <- which(Reprod[,"C"]>0 & fb==1)
num.b <- length(daughters.b)
num.f <- length(daughters.f)
shift_x <- shift_y <- rep(0,grid.size*n_taxa)
shift_x[daughters.b] <- sample(c(-params["dist",]:params["dist",]),num.b,replace=T) # vector of dispersal movements in x direction
shift_y[daughters.b] <- sample(c(-params["dist",]:params["dist",]),num.b,replace=T) # vector of dispersal movements in y direction
shift_x[daughters.f] <- 1 # Fungi always move positively in x direction
# vector of dispersal movements in y direction; constrained to one box away determined by probability "direct"
shift_y[daughters.f] <- sample(c(-1:1),num.f,replace=T,prob=c(0.5*(1-params["direct",]),params["direct",],0.5*(1-params["direct",])))
new_x <- (rep(1:params["x",],each=n_taxa,times=params["y",]) + shift_x + params["x",])%%params["x",] # calculate x coordinates of dispersal destinations
new_x[new_x==0] <- params["x",] # Substitute coordinates when there is no shift
new_y <- (rep(1:params["y",],each=params["x",]*n_taxa) + shift_y + params["y",])%%params["y",] # calculate y coordinates of dispersal destinations
new_y[new_y==0] <- params["y",]
Reprod[daughters.f,] <- eMF[daughters.f,] # set all fungi equal to their grid averages for translocation before colonization
index_col <- n_taxa*((new_y-1)*params["x",]+(new_x-1))+z # convert x,y coordinates to a vector of destination locations
# Transfer cells to new locations and sum when two or more of the same taxa go to same location
for(i in c(daughters.b,daughters.f)) {
Colonization[index_col[i],] <- Colonization[index_col[i],] + Reprod[i,]
}
# Set Monomers to the grid average
Monomers <- expand(Monomers.grid/grid.size,grid.size)
# Translocate nutrients within fungal taxa
i <- Microbes[fb==1,"C"]>0
Microbes[fb==1,][i,] <- eMF[fb==1,][i,]
# Colonization of dispersing microbes
Microbes <- Microbes + Colonization
# Multiply Vmax values for each substrate by the quantity of each enzyme
tev[,] <- (as.vector(Enzymes)*etVmax)[vm.index]
# Equation for Michaelis-Menten enzyme catalysis
rss <- rowSums(Substrates)
Decay <- tev*rss/(Km+rss)
# Loop that pulls out each batch of required enzymes and sums across redundant enzymes
DecaySums <- NULL
for(i in 1:dim(ReqEnz)[3]) {
DecaySums <- cbind(rowSums(as.matrix(ReqEnz[,,i])*Decay),DecaySums)
}
# Assess the rate-limiting enzyme and set decay to that rate
# Compare to substrate available and take the min, allowing for a tolerance of 1e-9
# Link cellulose degradation to lignin concentration (LCI)
DecayRates <- pmin(pmin(DecaySums[,1],DecaySums[,2],na.rm=T),rss-1e-9*rss,na.rm=T)
ss7 <- rowSums(Substrates[is.lignin,])
DecayRates[is.cellulose] <- DecayRates[is.cellulose]*(1+params["LCI_slope",]*ss7/(Substrates[is.cellulose,1]+ss7))
# Monomer stoichiometry: organic monomers follow substrate
SubstrateRatios <- Substrates/rss
SubstrateRatios[rss==0,] <- 0
MonomerRatios[org,] <- SubstrateRatios
# Update monomer pools
# Decayed substrate stays in the grid box for the whole iteration
# Preferential access by resident microbes
# Any leftovers diffuse (get averaged) across grid at next iteration
Monomers[org,] <- Monomers[org,] + (DecayRates+MonInput[org])*MonomerRatios[org,]
Monomers[mineral,] <- Monomers[mineral,] + MonInput[mineral]*MonomerRatios[mineral,]
# Keep track of mass balance for inputs
Cum_Monomer <- Cum_Monomer + colSums(as.vector(MonInput)*MonomerRatios)
# Recalculate Monomer stoichiometry after changes due to dead microbial biomass
rsm <- rowSums(Monomers)
MonomerRatios[org,] <- Monomers[org,]/rsm[org]
MonomerRatios[org,][rsm[org]==0,] <- 0
# Multiply microbial biomass by each taxon's uptake allocation to get biomass of each uptake enzyme by taxon
rsi <- rowSums(Microbes)
MicCXGenes <- rsi*UptakeGenes
# Equation for hypothetical potential uptake (per unit of compatible uptake protein)
Potential_Uptake <- Uptake_ReqEnz*Uptake_Vmax*rsm/(Uptake_Km+rsm)
# Matrix multiplication to get max possible uptake by taxon and monomer
# Must extract each grid point separately for operation
for(i in 1:grid.size) {
i.a <- ((i-1)*n_monomers+1):(i*n_monomers)
i.b <- ((i-1)*n_taxa+1):(i*n_taxa)
Max_Uptake <- MicCXGenes[i.b,]%*%t(Potential_Uptake[i.a,])
Max_Uptake1[,i.a] <- Max_Uptake # wide format; mon*grid cols
Max_Uptake2[i.a,] <- t(Max_Uptake) # long format; mon*grid rows
}
# Sum the total potential uptake of each monomer
csmu <- colSums(Max_Uptake1)
# Take the min of the monomer available and the max potential uptake
Min_Uptake <- pmin(csmu,rsm,na.rm=T)
# Scale the uptake to what's available
Uptake <- Max_Uptake2*Min_Uptake/csmu
Uptake[csmu==0,] <- 0
# Prevent total uptake from getting too close to zero
Uptake <- Uptake - 1e-9*Uptake
# Calculate total uptake by monomer
Monomer_Uptake <- rowSums(Uptake)*MonomerRatios
# Calculate total uptake by taxon
TUC.mat[,] <- (Uptake*MonomerRatios[,"C"])[TU.index]
TUN.mat[,] <- (Uptake*MonomerRatios[,"N"])[TU.index]
TUP.mat[,] <- (Uptake*MonomerRatios[,"P"])[TU.index]
# Sum across monomers
Taxon_Uptake_C <- colSums(TUC.mat)
Taxon_Uptake_N <- colSums(TUN.mat)
Taxon_Uptake_P <- colSums(TUP.mat)
# Enzyme production: rows are taxa, cols are enzyme genes
# Two alternatives: proportional to biomass or proportional to uptake C
Taxon_Enzyme_Production <- Taxon_Uptake_C*Induce_Enzyme_C + rsi*Constit_Enzyme_C
# Constrain enzyme production if it will cost more N than currently available in microbial biomass
Enzyme_Cost_N <- colSums(t(Taxon_Enzyme_Production)*EnzAttrib[,"N_cost"])
# Only select entries with microbes present
i <- which(Enzyme_Cost_N>0)
Taxon_Enzyme_Production[i,] <- Taxon_Enzyme_Production[i,]*pmin(Microbes[i,"N"],Enzyme_Cost_N[i])/Enzyme_Cost_N[i]
# Total enzyme production across all taxa
EP.mat[,] <- Taxon_Enzyme_Production[EP.index]
Enzyme_Production <- colSums(EP.mat)
# Total enzyme carbon cost for each taxon
ttep <- t(Taxon_Enzyme_Production)
Enzyme_Maint <- colSums(ttep*EnzAttrib[,"Maint_cost"])
Enzyme_Cost <- colSums(ttep) + Enzyme_Maint
Enzyme_Cost_N <- colSums(ttep*EnzAttrib[,"N_cost"])
Enzyme_Cost_P <- colSums(ttep*EnzAttrib[,"P_cost"])
# Uptake enzyme production: cols are taxa, rows are enzymes
# Total uptake enzyme carbon cost for each taxon
Uptake_Maint <- colSums(t(rsi*Uptake_Cost))*params["Uptake_Maint_cost",]
# Calculate loss rates and random death
Enzyme_Loss <- params["Enzyme_Loss_Rate",]*Enzymes
# Update enzyme pools
Enzymes <- Enzymes + Enzyme_Production - Enzyme_Loss
# Update microbe pools
Microbes[,"C"] <- Microbes[,"C"] + Taxon_Uptake_C*CUE - Enzyme_Cost - Uptake_Maint
Microbes[,"N"] <- Microbes[,"N"] + Taxon_Uptake_N - Enzyme_Cost_N
Microbes[,"P"] <- Microbes[,"P"] + Taxon_Uptake_P - Enzyme_Cost_P
# Kill microbes that are starving and transfer biomass to substrate
Death <- Colonization.reset
if (is.na(sum(Microbes))) {break}
starve_index <- (Microbes[,"C"]<params["C_min",] | Microbes[,"N"]<params["N_min",] | Microbes[,"P"]<params["P_min",]) & Microbes[,"C"]>0
Death[starve_index,] <- Microbes[starve_index,]
Microbes[starve_index,] <- 0
# Determine which locations have microbes
is.m <- which(Microbes[,"C"]>0)
# Zero the death index
kill <- kill.reset
# Sample for death events (=1)
kill[is.m] <- runif(length(is.m)) < r.death[is.m]
# Kill microbes
killed <- kill*Microbes
Death <- Death + killed
Microbes <- Microbes - killed
# Adjust stoichiometry of microbial biomass
# Index locations of microbial cells
mic.index <- as.numeric(which(Microbes[,"C"]>0))
# drop=F retains matrix structure if only one row is selected
Mic.subset <- Microbes[mic.index,,drop=F]
# Calculate microbial ratios
MicrobeRatios <- (Mic.subset/rowSums(Mic.subset))
# Select the corresponding minimum quotas
MinRat <- MinRatios[mic.index,,drop=F]
# Index only microbes that have below-minimum quotas
rat.index <- mic.index[MicrobeRatios[,"C"] < MinRat[,"C"] | MicrobeRatios[,"N"] < MinRat[,"N"] | MicrobeRatios[,"P"] < MinRat[,"P"] ]
Mic.subset <- Microbes[rat.index,,drop=F]
StartMicrobes <- Mic.subset
MicrobeRatios <- Excess <- (Mic.subset/rowSums(Mic.subset))
MinRat <- MinRatios[rat.index,,drop=F]
# Calculate difference between min and actual ratios
Deficit <- Deficit.0 <- MicrobeRatios - MinRat
# Determine the limiting nutrient that will be conserved
Limiting <- max.col(-Deficit/MinRat,ties.method="first")
# Set all deficient ratios to their minima
MicrobeRatios[Deficit<0] <- MinRat[Deficit<0]
# Reduce the mass fractions for non-deficient elements in proportion to the distance from the minimum
# Calculate how far above the minimum each non-deficient element is
Excess[Deficit>0] <- Deficit[Deficit>0]
# Set deficient element fractions to zero
Excess[Deficit<0] <- 0
# Partition the total deficit to the excess element(s) in proportion to their distances from their minima
Deficit.0[Deficit>0] <- 0
MicrobeRatios[Deficit>0] <- MicrobeRatios[Deficit>0] + (rowSums(Deficit.0)*Excess/rowSums(Excess))[Deficit>0]
# Construct hypothetical nutrient quotas for each possible minimum nutrient
MC <- Mic.subset[,"C"]; MN <- Mic.subset[,"N"]; MP <- Mic.subset[,"P"]
MRC <- MicrobeRatios[,"C"]; MRN <- MicrobeRatios[,"N"]; MRP <- MicrobeRatios[,"P"]
new.C <- cbind(MC,MN*MRC/MRN,MP*MRC/MRP)
new.N <- cbind(MC*MRN/MRC,MN,MP*MRN/MRP)
new.P <- cbind(MC*MRP/MRC,MN*MRP/MRN,MP)
# Insert the appropriate set of nutrient quotas scaled to the minimum nutrient
select <- cbind(1:length(rat.index),Limiting)
Microbes[rat.index,] <- cbind(new.C[select],new.N[select],new.P[select])
# Sum up the element losses from biomass across whole grid and calculate average loss
MicLoss <- StartMicrobes - Microbes[rat.index,]
# Update substrate pools
# Substrate stoichiometry
rss <- rowSums(Substrates)
SubstrateRatios <- Substrates/rss
SubstrateRatios[rss==0,] <- 0
# Add inputs and remove decay
Substrates <- Substrates + SubInput - DecayRates*SubstrateRatios
# Account for inputs in mass balance
Cum_Substrate <- Cum_Substrate + colSums(SubInput)
# Add dead microbial biomass
Death.mat[,] <- Death
Substrates[is.deadMic,] <- Substrates[is.deadMic,] + colSums(Death.mat)
# Add dead enzymes
DeadEnz.mat[,] <- c(Enzyme_Loss,Enzyme_Loss*EnzAttrib[,"N_cost"],Enzyme_Loss*EnzAttrib[,"P_cost"])
Substrates[is.deadEnz,] <- Substrates[is.deadEnz,] + colSums(DeadEnz.mat)
# Update monomer pools
Monomers[is.NH4,"N"] <- Monomers[is.NH4,"N"] + sum(MicLoss[,"N"])/grid.size
Monomers[is.PO4,"P"] <- Monomers[is.PO4,"P"] + sum(MicLoss[,"P"])/grid.size
Monomers <- Monomers - Monomer_Uptake
# Sum microbes prior to reproduction
Microbes.grid <- sum.grid(Microbes,Mic.names,grid.size)
# Reset the vector of fungal locations
Fungi.vec <- Fungi.vec.zero
# Add one or two fungi to the count vector based on size
Fungi.vec[fb==1][Microbes[fb==1,"C"]>0] <- 1
Fungi.vec[fb==1][Microbes[fb==1,"C"]>params["max_size_f",]] <- 2
# Sum up the count vector
Fungi.count <- sum.grid(Fungi.vec,Mic.names,grid.size)
# Cell division
MicrobesBeforeDivision <- Microbes
Microbes[fb==0,][Microbes[fb==0,"C"]>params["max_size_b",],] <- Microbes[fb==0,][Microbes[fb==0,"C"]>params["max_size_b",],]/2
Microbes[fb==1,][Microbes[fb==1,"C"]>params["max_size_f",],] <- Microbes[fb==1,][Microbes[fb==1,"C"]>params["max_size_f",],]/2
# Add daughter cells to matrix of new reproduction
ReprodNew <- MicrobesBeforeDivision-Microbes
# Sum pools across grid
Respiration <- sum(Taxon_Uptake_C*(1-CUE)) + sum(Enzyme_Maint) + sum(Uptake_Maint) + sum(MicLoss[,"C"])
Monomers.grid <- sum.grid(Monomers,Mon.names,grid.size)
Enzymes.grid <- sum.grid(Enzymes,Enz.names,grid.size)
Substrates.grid <- sum.grid(Substrates,Sub.names,grid.size)
# Update monomer leaching
Leaching_N <- Monomers.grid["NH4","N"]*params["Leaching",]*exp(params["Psi.slope.leach",]*Psi[i_t])
Leaching_P <- Monomers.grid["PO4","P"]*params["Leaching",]*exp(params["Psi.slope.leach",]*Psi[i_t])
# Keep track of mass balance
Cum_Leaching_N <- Cum_Leaching_N + Leaching_N
Cum_Leaching_P <- Cum_Leaching_P + Leaching_P
Monomers.grid["NH4","N"] <- Monomers.grid["NH4","N"] - Leaching_N
Monomers.grid["PO4","P"] <- Monomers.grid["PO4","P"] - Leaching_P
# Record total pools for each time step
RespSeries <- c(RespSeries,Respiration)
MonomersSeries <- rbind(MonomersSeries,Monomers.grid[,"C"])
NH4Series <- rbind(NH4Series,Monomers.grid["NH4","N"])
PO4Series <- rbind(PO4Series,Monomers.grid["PO4","P"])
EnzymesSeries <- rbind(EnzymesSeries,t(Enzymes.grid))
SubstratesSeries <- rbind(SubstratesSeries,Substrates.grid[,"C"])
MicrobesSeries <- rbind(MicrobesSeries,Microbes.grid[,"C"])
MicrobesNSeries <- rbind(MicrobesNSeries,Microbes.grid[,"C"]/Microbes.grid[,"N"])
MicrobesPSeries <- rbind(MicrobesPSeries,Microbes.grid[,"C"]/Microbes.grid[,"P"])
MicrobesNPSeries <- rbind(MicrobesNPSeries,Microbes.grid[,"N"]/Microbes.grid[,"P"])
# Sum microbial pools
Mic_Sum <- rbind(Mic_Sum,colSums(Microbes.grid))
# Output grid if enzyme mass has peaked
if (i_t>3 & params["output.mid",]==1) {
if (Enzyme_Sum[i_t]<Enzyme_Sum[i_t-1] & Enzyme_Sum[i_t-1]>Enzyme_Sum[i_t-2]) {
grid.max <- list("Substrates"=Substrates,"Monomers"=Monomers,"Enzymes"=Enzymes,"Microbes"=Microbes)
grid.time <- i_t
}
}
# Output grid images
if (i_t == 1 || i_t%%params["print.grid",]==0){
apply(matrix(1:n_taxa),1,function(i)plot.grid(Microbes[n_taxa*(0:(grid.size-1))+i,"C"],params["x",],zlim=72*fb[1:n_taxa][i]+3,paste(i,":",fb[1:n_taxa][i],":",i_t,sep="")))
}
# Sum substrate, monomer, enzyme pools
Substrate_Sum <- c(Substrate_Sum,sum(Substrates.grid[,"C"]))
Monomer_Sum <- c(Monomer_Sum,sum(Monomers.grid[,"C"]))
Enzyme_Sum <- c(Enzyme_Sum,sum(Enzymes.grid))
print(i_t)
}
# Calculate recoveries
Recovered_C <- (sum(RespSeries) + sum(Substrates.grid[,"C"]) + sum(Monomers.grid[,"C"]) + sum(Enzymes.grid) + sum(Microbes.grid[,"C"]))-(Cum_Microbe["C"] + Cum_Substrate["C"] + Cum_Monomer["C"] + Cum_Enzyme_C)
Recovered_N <- (Cum_Leaching_N + sum(Substrates.grid[,"N"]) + sum(Monomers.grid[,"N"]) + sum(Enzymes.grid*EnzAttrib[,"N_cost"]) + sum(Microbes.grid[,"N"]))-(Cum_Microbe["N"] + Cum_Substrate["N"] + Cum_Monomer["N"] + Cum_Enzyme_N)
Recovered_P <- (Cum_Leaching_P + sum(Substrates.grid[,"P"]) + sum(Monomers.grid[,"P"]) + sum(Enzymes.grid*EnzAttrib[,"P_cost"]) + sum(Microbes.grid[,"P"]))-(Cum_Microbe["P"] + Cum_Substrate["P"] + Cum_Monomer["P"] + Cum_Enzyme_P)
index <- is.finite(log10(apply(MicrobesSeries,2,max)))
correl <- cor.test(rowSums(EnzGenes)[index],log(apply(MicrobesSeries,2,max)/grid.size)[index])
# Output grid at end if necessary
if (is.null(grid.max)) {
grid.max <- list("Substrates"=Substrates,"Monomers"=Monomers,"Enzymes"=Enzymes,"Microbes"=Microbes)
grid.time <- params["end_time",]
}
# Correlations of taxon abundance from the grid
corr.vec <- matrix(NA,ncol=2,nrow=n_taxa*n_taxa)
count <- 1
taxon.mat <- matrix(grid.max$Microbes[,"C"],ncol=n_taxa,byrow=T)
for (i in 1:n_taxa) {
for (j in 1:n_taxa) {
taxon.a <- taxon.mat[,i]
taxon.b <- taxon.mat[,j]
overlap <- which(taxon.a>0 & taxon.b>0)
taxon.1 <- taxon.a[overlap]
taxon.2 <- taxon.b[overlap]
if (length(taxon.1)>2) {
corr.list <- cor.test(taxon.1,taxon.2)
} else {
corr.list <- list("estimate"=NA,"p.value"=NA)
}
corr.vec[count,1] <- corr.list$estimate
corr.vec[count,2] <- corr.list$p.value*((params["n_taxa",]*params["n_taxa",]-params["n_taxa",])/2) # Bonferroni corrected p-values
count <- count+1
}
}
# Output pulse results
out.pulse <- list(
"params"=params,
"seed"=rng.seed,
"Ea"=Ea,
"Uptake_Ea"=Uptake_Ea,
"Vmax"=Vmax0,
"Uptake_Vmax"=Uptake_Vmax0,
"Km"=Km0[1:n_substrates,],
"Uptake_Km"=Uptake_Km0[1:n_monomers,],
"ReqEnz"=list(ReqEnz[,,1][1:n_substrates,],ReqEnz[,,2][1:n_substrates,]),
"Uptake_ReqEnz"=Uptake_ReqEnz[1:n_monomers,],
"OptimalRatios"=OptimalRatios,
"RangeRatios"=RangeRatios,
"EnzGenes"=EnzGenes,
"EnzProdInduce"=EnzProdInduce,
"EnzProdConstit"=EnzProdConstit,
"UptakeGenes"=UptakeGenes[1:n_taxa,],
"UptakeGenesForEnz"=UptakeGenesForEnz,
"EnzAttrib"=EnzAttrib,
"MonInput"=MonInput[1:n_monomers],
"SubInput"=SubInput[1:n_substrates,],
"CUE"=CUE.ref,
"RespSeries"=RespSeries,
"EnzymesSeries"=EnzymesSeries,
"SubstratesSeries"=SubstratesSeries,
"MonomersSeries"=MonomersSeries,
"NH4Series"=NH4Series,
"PO4Series"=PO4Series,
"MicrobesSeries"=MicrobesSeries,
"MicrobesNSeries"=MicrobesNSeries,
"MicrobesPSeries"=MicrobesPSeries,
"MicrobesNPSeries"=MicrobesNPSeries,
"grid"=grid.max,
"grid.time"=grid.time,
"grid.size"=grid.size,
"corr.vec"=corr.vec,
"cor.test"=correl,
"end_time"=params["end_time",],
"Mic_Sum"=Mic_Sum,
"Substrate_Sum"=Substrate_Sum,
"Monomer_Sum"=Monomer_Sum,
"Enzyme_Sum"=Enzyme_Sum,
"timestamp"=timestamp,
"fb"=fb[1:n_taxa],
"Temp"=Temp,
"Psi"=Psi,
"Tolerance"=Tolerance,
"RecoveredC_Resp_Sub_Mon_Enz_Mic"=
c(Recovered_C,sum(RespSeries),sum(Substrates.grid[,"C"]),sum(Monomers.grid[,"C"]),sum(Enzymes.grid),sum(Microbes.grid[,"C"])),
"RecoveredN_Leach_Sub_Mon_Enz_Mic"=
c(Recovered_N,Cum_Leaching_N,sum(Substrates.grid[,"N"]),sum(Monomers.grid[,"N"]),sum(Enzymes.grid*EnzAttrib[,"N_cost"]),sum(Microbes.grid[,"N"])),
"RecoveredP_Leach_Sub_Mon_Enz_Mic"=
c(Recovered_P,Cum_Leaching_P,sum(Substrates.grid[,"P"]),sum(Monomers.grid[,"P"]),sum(Enzymes.grid*EnzAttrib[,"P_cost"]),sum(Microbes.grid[,"P"]))
)
out.pulse
} # End of RunPulse() #########################################################
###############################################################################
TraitModel <- function(job.time,task.ID){
# Read in parameters
params.file <- paste("params/params",task.ID,".txt",sep="")
params1 <- read.table(params.file, header=F, row.names=1, sep="\t",stringsAsFactors=F)
params <- params1
grid.size <- params["x",]*params["y",]
seed.val <- as.POSIXct(job.time,format="%y%m%d%H%M%S") + task.ID # convert to date integer, add task id, and use as seed
set.seed(seed.val)
timestamp <- format(seed.val, "%y%m%d%H%M%S") # convert date seed to 12 digit number
rng.seed <- as.numeric(timestamp)
if(params["set.seed",]==1) {
set.seed(as.POSIXct(as.character(params["seed",]),format="%y%m%d%H%M%S"))
rng.seed <- params["seed",]
}
n_monomers <- params["n_substrates",]+2
n_genes <- params["n_enzymes",]
n_upgenes <- params["n_uptake",]
# Read in activation energies for substrates
Ea.frame <- read.table("Ea.txt", header=TRUE, row.names=1, sep="\t",stringsAsFactors=F)
# Choose taxa with "fungal" strategy, specifically ability to translocate C and nutrients
fb <- sample(c(1,0),params["n_taxa",],replace=T,prob=c(params["fb",],(1-params["fb",])))
# Set stress tolerance allocation for each taxon
Tolerance <- c(runif(params["n_taxa",],params["Tol_min",],params["Tol_max",]))
if(params["Tol_min",]==params["Tol_max",]) {dummy <- runif(params["n_taxa",])}
# Pre-exponential constants for enzymes
# Rows are substrates; cols are enzymes
Vmax0 <- c(runif(params["n_substrates",]*params["n_enzymes",],params["Vmax0_min",],params["Vmax0_max",]))
Vmax0 <- matrix(Vmax0,nrow=params["n_substrates",],ncol=params["n_enzymes",],byrow=T,dimnames=list(sprintf("%s%03d","Sub",1:params["n_substrates",]),sprintf("%s%03d","Enz",1:params["n_enzymes",])))
# Pre-exponential constants for uptake
# Rows are monomers; cols are uptake enzymes
Uptake_Vmax0 <- c(runif(params["n_uptake",]*n_monomers,params["Uptake_Vmax0_min",],params["Uptake_Vmax0_max",]))
Uptake_Vmax0 <- matrix(Uptake_Vmax0,nrow=n_monomers,ncol=params["n_uptake",],byrow=T,dimnames=list(sprintf("%s%03d","Mon",1:n_monomers),sprintf("%s%03d","Upt",1:params["n_uptake",])))
# Enzyme specificity matrix of activation energies
# Rows are substrates; cols are enzymes
Ea <- sapply(seq(len=dim(Ea.frame)[1]),function(i)runif(params["n_enzymes",],Ea.frame$Ea_min[i],Ea.frame$Ea_max[i]))
Ea <- matrix(Ea,nrow=params["n_substrates",],ncol=params["n_enzymes",],byrow=T,dimnames=list(sprintf("%s%03d","Sub",1:params["n_substrates",]),sprintf("%s%03d","Enz",1:params["n_enzymes",])))
# Uptake specificity matrix of activation energies
# Rows are monomers; cols are uptake enzymes
Uptake_Ea <- c(runif(params["n_uptake",]*n_monomers,params["Uptake_Ea_min",],params["Uptake_Ea_max",]))
Uptake_Ea <- matrix(Uptake_Ea,nrow=n_monomers,ncol=params["n_uptake",],byrow=T,dimnames=list(sprintf("%s%03d","Mon",1:n_monomers),sprintf("%s%03d","Upt",1:params["n_uptake",])))
# Enzymes required for substrate degradation
# Rows are substrates; cols are enzymes
# Same number within row implies redundancy
# Ensures each substrate is degraded by at least 1 enzyme and every enzyme degrades at least 1 substrate
probability_vector <- rep(0,params["n_enzymes",])
probability_vector[1:params["Enzymes_per_sub",]] <- 1
ReqEnz1 <- NULL
for(i in 1:params["n_substrates",]) {
ReqEnz1 <- c(ReqEnz1,sample(probability_vector,params["n_enzymes",],replace=F))
}
ReqEnz1 <- matrix(ReqEnz1,nrow=params["n_substrates",],ncol=params["n_enzymes",],byrow=T)
probability_vector <- rep(0,params["n_substrates",])
probability_vector[1] <- 1
for (i in 1:params["n_enzymes",]) {
if(sum(ReqEnz1[,i])==0) ReqEnz1[,i] <- sample(probability_vector,params["n_substrates",],replace=F)
}
# Choose some substrates that require multiple enzymes
probability_vector <- rep(0,params["n_enzymes",])
if(params["Avg_extra_req_enz",]>0) probability_vector[1:params["Avg_extra_req_enz",]] <- 1
ReqEnz2 <- sample(probability_vector,params["n_substrates",]*params["n_enzymes",],replace=T)
ReqEnz2[ReqEnz1==1] <- 0
ReqEnz2 <- matrix(ReqEnz2,nrow=params["n_substrates",],ncol=params["n_enzymes",],byrow=T)
for(j in 1:dim(ReqEnz2)[1]) {
# Put in NAs if the substrate does not require a second enzyme
if(rowSums(ReqEnz2)[j]==0) ReqEnz2[j,] <- NA
}
# Generate the correct structure for the required enzyme matrices
ReqEnz <- c(expand(ReqEnz1,grid.size),expand(ReqEnz2,grid.size))
# rows, cols, stacks; default assemble by cols
ReqEnz <- array(ReqEnz,dim=c(grid.size*params["n_substrates",],params["n_enzymes",],2),dimnames=list(rep(sprintf("%s%03d","Sub",1:params["n_substrates",]),grid.size),sprintf("%s%03d","Enz",1:params["n_enzymes",]),c("set1","set2")))
# Enzymes used for uptake
# Rows are monomers; cols are uptake enzymes
# Same number within row implies redundancy
# Make sure each monomer is taken up by at least one transporter and every transporter takes up at least one monomer
probability_vector <- rep(0,params["n_uptake",])
probability_vector[1:params["Uptake_per_monomer",]] <- 1
Uptake_ReqEnz1 <- NULL
for(i in 1:n_monomers) {
Uptake_ReqEnz1 <- c(Uptake_ReqEnz1,sample(probability_vector,params["n_uptake",],replace=F))
}
Uptake_ReqEnz <- matrix(Uptake_ReqEnz1,nrow=n_monomers,ncol=params["n_uptake",],byrow=T,dimnames=list(sprintf("%s%03d","Mon",1:n_monomers),sprintf("%s%03d","Upt",1:params["n_uptake",])))
probability_vector <- rep(0,n_monomers)
probability_vector[1] <- 1
for (i in 1:params["n_uptake",]) {
if(sum(Uptake_ReqEnz[,i])==0) Uptake_ReqEnz[,i] <- sample(probability_vector,n_monomers,replace=F)
}
# Substrate concentrations
substrates.frame <- read.table("substrates.txt", header=TRUE, sep="\t",stringsAsFactors=F,row.names=1)
Substrates <- data.matrix(substrates.frame)
# Substrate input rates
SubInputC <- read.table("inputs.txt", header=TRUE, sep="\t",stringsAsFactors=F,row.names=1)[,1]
SubInputN <- SubInputC*Substrates[,"N"]/Substrates[,"C"]
SubInputP <- SubInputC*Substrates[,"P"]/Substrates[,"C"]
SubInput <- c(SubInputC,SubInputN,SubInputP)
SubInput <- matrix(SubInput,nrow=params["n_substrates",],ncol=3,byrow=F,dimnames=list(sprintf("%s%03d","SubIn",1:params["n_substrates",]),c("C","N","P")))
SubInput["SubIn001",] <- SubInput["SubIn002",] <- 0
Substrates["DeadMic",] <- Substrates["DeadEnz",] <- 0
# Monomers produced by each substrate (binary)
# Rows are substrates; cols are monomers
# All rows should sum to 1
MonomersProduced <- cbind(rep(0,params["n_substrates",]),rep(0,params["n_substrates",]),diag(params["n_substrates",]))
MonomersProduced <- as.matrix(MonomersProduced)
dimnames(MonomersProduced) <- list(sprintf("%s%03d","Sub",1:params["n_substrates",]),sprintf("%s%03d","Mon",-1:(n_monomers-2)))
dimnames(MonomersProduced)[[2]][1:4]<- c("NH4","PO4","DeadMic","DeadEnz")
# Optimal stoichiometry of bacterial taxa
OptimalRatios <- rep(c(params["Cfrac_b",],params["Nfrac_b",],params["Pfrac_b",]),params["n_taxa",])
OptimalRatios <- matrix(OptimalRatios,nrow=params["n_taxa",],ncol=3,byrow=T,dimnames=list(sprintf("%s%03d","Tax",1:params["n_taxa",]),c("C","N","P")))
# Substitute fungal stoichiometry
OptimalRatios[fb==1,] <- matrix(rep(c(params["Cfrac_f",],params["Nfrac_f",],params["Pfrac_f",]),sum(fb)),ncol=3,byrow=T)
RangeRatios <- rep(c(params["Crange",],params["Nrange",],params["Prange",]),params["n_taxa",])
RangeRatios <- matrix(RangeRatios,nrow=params["n_taxa",],ncol=3,byrow=T,dimnames=list(sprintf("%s%03d","Tax",1:params["n_taxa",]),c("C","N","P")))
RangeWeight <- RangeRatios/OptimalRatios
# Calculate minimum cell quotas
MinRatios <- OptimalRatios-RangeRatios
# Enzyme genes possessed by taxa (binary)
# Rows are taxa; cols are genes
# Not all taxa need have an enzyme gene
# This method draws the number of genes per taxon from a uniform distribution
# and allows production rates between Enz_Prod_min and max
genes_per_taxon <- sample(params["Enz_per_taxon_min",]:params["Enz_per_taxon_max",],params["n_taxa",],replace=T)
EnzGenes <- sapply(seq(params["n_taxa",]),function(i){
probability_vector <- rep(0,n_genes)
probability_vector[0:genes_per_taxon[i]] <- 1
sample(probability_vector,n_genes,replace=F)})
EnzGenes <- matrix(EnzGenes,nrow=params["n_taxa",],ncol=n_genes,byrow=T,dimnames=list(sprintf("%s%03d","Tax",1:params["n_taxa",]),sprintf("%s%03d","Enz",1:n_genes)))
Normalize <- 1
if (params["NormalizeProd",]==1) {
Normalize <- rowSums(EnzGenes)/params["Enz_per_taxon_max",]
}
EnzProdInduce <- apply(EnzGenes,1,function(x)runif(1,params["Enz_Prod_min",],params["Enz_Prod_max",]))
EnzProdInduce <- EnzProdInduce*EnzGenes/Normalize
EnzProdInduce[is.na(EnzProdInduce)] <- 0
EnzProdConstit <- apply(EnzGenes,1,function(x)runif(1,params["Constit_Prod_min",],params["Constit_Prod_max",]))
EnzProdConstit <- EnzProdConstit*EnzGenes/Normalize
EnzProdConstit[is.na(EnzProdConstit)] <- 0
# Uptake genes possessed by taxa (binary)
# Rows are taxa; cols are genes
# All taxa must have at least one uptake gene
# For each taxon, check if it has the enz gene associated with the monomer(s) to be targeted by the uptake gene
# If so, assign that uptake gene to the taxon with probability p
# Enz-Taxon%*%Enz-Sub%*%Sub-Monomer to get Taxon-Monomer matrix
# Matrix multiply this by the required uptake enzyme matrix, then assign uptake genes
RE2 <- ReqEnz[,,2][1:params["n_substrates",],]
RE2[is.na(RE2)] <- 0
enz_sub <- ReqEnz[,,1][1:params["n_substrates",],] + RE2
# Matrix multiplication to relate taxa to the monomers they can generate with their enzymes
UptakeGenes <- EnzGenes%*%t(enz_sub)%*%MonomersProduced
UptakeGenes[,1:2] <- 1
UptakeGenes[UptakeGenes>0] <- 1
# Make sure every taxon is likely to have an uptake enzyme for at least 1 organic monomer
# Not guaranteed unless uptake_prob = 1
probability_vector <- rep(0,n_monomers-2)
probability_vector[1] <- 1
for (i in 1:params["n_taxa",]) {
if(sum(UptakeGenes[i,3:n_monomers])==0) UptakeGenes[i,3:n_monomers] <- sample(probability_vector,n_monomers-2,replace=F)
}
# Give taxa random number of additional uptake genes between the number they have and n_upgenes
for (i in 1:params["n_taxa",]) {
n.zero <- length(UptakeGenes[i,][UptakeGenes[i,]==0])
probability_vector <- rep(0,n.zero)
probability_vector[1:sample(0:n.zero,1)] <- 1
UptakeGenes[i,][UptakeGenes[i,]==0] <- sample(probability_vector,n.zero,replace=F)
}
UptakeGenesForEnz <- UptakeGenes
# If true then the uptake potential is normalized to the number of uptake genes
Normalize <- 1
if (params["NormalizeUptake",]==1) {
Normalize <- rowSums(UptakeGenes)/params["n_uptake",]
}
UptakeGenes <- UptakeGenes/Normalize
# Choose total amount of uptake allocation at random
UptakeProd <- c(runif(params["n_taxa",],params["Uptake_C_cost_min",],params["Uptake_C_cost_max",]))
UptakeGenes <- UptakeGenes*UptakeProd
# Calculate CUE as a function of the number of enzyme genes, uptake genes, and reference CUE
CUE.ref <- params["CUE_enz",]*rowSums(EnzGenes)/(params["Enz_per_taxon_max",]) +
params["CUE_uptake",]*rowSums(UptakeGenesForEnz)/(params["n_uptake",]) + params["CUE_ref",] + params["Tol_CUE",]*Tolerance
# Enzyme attributes
# Rows are enzymes; cols are attributes
EnzAttrib_C <- rep(params["Enz_C_cost",],params["n_enzymes",])
EnzAttrib <- rbind(EnzAttrib_C,params["Enz_N_cost",],params["Enz_P_cost",],params["Enz_Maint_cost",])
EnzAttrib <- matrix(EnzAttrib,nrow=params["n_enzymes",],ncol=4,byrow=T,dimnames=list(sprintf("%s%03d","Enz",1:params["n_enzymes",]),c("C_cost","N_cost","P_cost","Maint_cost")))
# Monomer pool sizes for all elements
Monomers <- rbind(c(0,params["Init_NH4",],0),c(0,0,params["Init_PO4",]),Substrates*(params["Monomer_Substrate_Ratio",]))
Monomers[is.na(Monomers)] <- 0
Monomers <- as.matrix(Monomers)
dimnames(Monomers) <- list(c("NH4","PO4","DeadMic","DeadEnz",sprintf("%s%03d","Mon",3:(n_monomers-2))),c("C","N","P"))
# Monomer input rates
MonInput <- c(params["Input_NH4",],params["Input_PO4",],read.table("inputs.txt", header=TRUE, sep="\t",stringsAsFactors=F,row.names=1)[,2])
MonInput <- as.matrix(MonInput)
rownames(MonInput) <- rownames(Monomers)
# Enzyme pool sizes for all elements
Enzymes <- matrix(c(runif(params["n_enzymes",],params["Enz_min",],params["Enz_max",])),nrow=params["n_enzymes",],ncol=1,dimnames=list(sprintf("%s%03d","Enz",1:params["n_enzymes",]),"C"))
# Microbial pool sizes for all elements
BacC <- 0.5*params["max_size_b",]
BacN <- BacC*params["Nfrac_b",]/params["Cfrac_b",]
BacP <- BacC*params["Pfrac_b",]/params["Cfrac_b",]
FunC <- 0.5*params["max_size_f",]
FunN <- FunC*params["Nfrac_f",]/params["Cfrac_f",]
FunP <- FunC*params["Pfrac_f",]/params["Cfrac_f",]
Microbes <- rep(c(BacC,BacN,BacP),params["n_taxa",]*grid.size)
Microbes <- matrix(Microbes,nrow=params["n_taxa",]*grid.size,ncol=3,byrow=T,dimnames=list(rep(sprintf("%s%03d","Tax",1:params["n_taxa",]),grid.size),c("C","N","P")))
fb <- rep(fb,grid.size)
Microbes[fb==1,"C"] <- FunC; Microbes[fb==1,"N"] <- FunN; Microbes[fb==1,"P"] <- FunP
# Randomly assign microbes to each grid box
p.b <- params["taxa_per_box",]
p.f <- p.b*params["max_size_b",]/params["max_size_f",]
choose_taxa <- sample(c(1,0),params["n_taxa",]*grid.size,replace=T,c(p.b,(1-p.b)))
choose_taxa[fb==1] <- sample(c(1,0),sum(fb),replace=T,c(p.f,(1-p.f)))
Microbes[choose_taxa==0,] <- 0
# Initialize the reproduction list
ReprodNew <- Microbes
ReprodNew[ReprodNew>0] <- 0
Colonization.reset <- ReprodNew
# Account for efficiency-specificity tradeoff by dividing Vmax_0 by the number of substrates (or monomers) targeted
# and multiplied by a specificity factor
# Leave Km unchanged, effectively increasing it by the factor that Vmax_0 is reduced
total_substrates <- colSums(ReqEnz[,,1][1:params["n_substrates",],]) + colSums(RE2)
if (params["Specif_factor",]==0) {total_substrates[total_substrates>1] <- 1} else
{total_substrates[total_substrates>1] <- total_substrates[total_substrates>1]*params["Specif_factor",]}
Vmax0 <- t(t(Vmax0)/total_substrates)
Vmax0[!is.finite(Vmax0)] <- 0
total_monomers <- colSums(Uptake_ReqEnz)
if (params["Specif_factor",]==0) {total_monomers[total_monomers>1] <- 1} else
{total_monomers[total_monomers>1] <- total_monomers[total_monomers>1]*params["Specif_factor",]}
Uptake_Vmax0 <- t(t(Uptake_Vmax0)/total_monomers)
Uptake_Vmax0[!is.finite(Uptake_Vmax0)] <- 0
# Implement Vmax-Km tradeoff as a direct correlation with slope Vmax_Km,
# and error term normally distributed with magnitude Km_error. Minimum Km constrained to Km_min
Km <- abs(rnorm(Vmax0, mean=Vmax0*params["Vmax_Km",], sd = mean(Vmax0)*params["Km_error",])+params["Vmax_Km_int",])
Km[Km < params["Km_min",]] <- params["Km_min",]
Km <- matrix(Km,nrow=params["n_substrates",],ncol=params["n_enzymes",],byrow=F,dimnames=list(sprintf("%s%03d","Sub",1:params["n_substrates",]),sprintf("%s%03d","Enz",1:params["n_enzymes",])))
Uptake_Km <- abs(rnorm(Uptake_Vmax0, mean=Uptake_Vmax0*params["Uptake_Vmax_Km",], sd = mean(Uptake_Vmax0)*params["Km_error",])+params["Uptake_Vmax_Km_int",])
Uptake_Km[Uptake_Km < params["Uptake_Km_min",]] <- params["Uptake_Km_min",]
Uptake_Km <- matrix(Uptake_Km,nrow=n_monomers,ncol=params["n_uptake",],byrow=F,dimnames=list(sprintf("%s%03d","Mon",1:n_monomers),sprintf("%s%03d","Upt",1:params["n_uptake",])))
# Expand matricies to cover the whole grid
Substrates <- expand(Substrates,grid.size)
SubInput <- expand(SubInput,grid.size)
Monomers <- expand(Monomers,grid.size)
MonInput <- expand(MonInput,grid.size)
Enzymes <- expand(Enzymes,grid.size)
Km0 <- expand(Km,grid.size)
Uptake_Km0 <- expand(Uptake_Km,grid.size)
Uptake_ReqEnz <- expand(Uptake_ReqEnz,grid.size)
MinRatios <- expand(MinRatios,grid.size)