% Supplementary source codes for 
% "Decomposition by ectomycorrhizal fungi alters
%   soil carbon storage in a simulation model"
% by Moore J.A.M., Post W.M., Jiang J., and Classen A.T.

% This file solve the equations by ode solver 

% Time span
tspan = [0 1000]; 

% Initial C pool
y0 = [10000, 1000, 1000, 50, 100, 5000];

% Parameters
G = 1000;       % Net primary productivity
aL = 0.05;      % Leaf litter allocation
aR = 0.03;      % Root allocation
aM = 0.02;      % Mycorrhizal allocation
kL = 0.5;       % Litter turnover
kR = 0.3;       % Root turnover
kM = 3;         % Mycorrhizal turnover
kB = 3.75;      % Microbial turnover
rL = 0.9;       % Fraction of litter respired
rM = 1;         % Mycorrhizal respiration
rB = 5;         % Microbial respiration
eP = 1;         % Plant acquisition efficiency
eS = 0.3;       % Soil acquisition efficienty
u = 1.64;       % Specific mycorrhizal C uptake
V1 = 1600;      % Maximum microbial C uptake
K1 = 5000;      % Microbial half-saturation
V2 = 200;       % Maximum mycorrhizal C uptake
K2 = 5000;      % Mycorrhizal half-saturation
s = 1;          % scales NPP (any number between 0 and 1)
flag = 2;       % Indicate which C uptake strategy is used (1,2,or 3)
% Another version of parameters can be K1=5000, K2=2000, and u=2.29
pars = [G aL aR aM kL kR kM kB rL rM rB eP eS u V1 K1 V2 K2 s flag];

% Solve the equation by matlab ode45 solver, need to call "myc_ode.m"
[t,y] = ode45(@(t,y) myc_ode(t,y,pars),tspan,y0);


% Plot results
figure
subplot(3,2,1)
plot(t,y(:,1))
ylabel('P')

subplot(3,2,2)
plot(t,y(:,6))
ylabel('S')

subplot(3,2,3)
plot(t,y(:,2))
ylabel('L')

subplot(3,2,4)
plot(t,y(:,3))
ylabel('R')

subplot(3,2,5)
plot(t,y(:,4))
ylabel('M')

subplot(3,2,6)
plot(t,y(:,5))
ylabel('B')