% 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 contain function to be called by matlab ode solver
function ydot = myc_ode(t,y,pars)

% Assign the values to each variable
P = y(1);   % Above ground plant C biomass
L = y(2);   % Leaf litter
R = y(3);   % Dead roots
M = y(4);   % Mycorrhizal C biomass
B = y(5);   % Bacterial C biomass
S = y(6);   % Soil C storage

% Parameters 
G = pars(1);    % Net primary productivity
aL = pars(2);   % Leaf litter allocation
aR = pars(3);   % Root allocation
aM = pars(4);   % Mycorrhizal allocation
kL = pars(5);   % Litter turnover
kR = pars(6);   % Root turnover
kM = pars(7);   % Mycorrhizal turnover
kB = pars(8);   % Microbial turnover
rL = pars(9);   % Fraction of litter respired
rM = pars(10);  % Mycorrhizal respiration
rB = pars(11);  % Microbial respiration
eP = pars(12);  % Plant acquisition efficiency
eS = pars(13);  % Soil acquisition efficienty
u = pars(14);   % Specific mycorrhizal C uptake
V1 = pars(15);  % Maximum microbial C uptake
K1 = pars(16);  % Microbial half-saturation
V2 = pars(17);  % Maximum mycorrhizal C uptake
K2 = pars(18);  % Mycorrhizal half-saturation
s = pars(19);   % scales NPP
flag = pars(20);% Indicate which C uptake strategy is used

if flag == 1
    VM = 0;     % No C uptake from soil, corresponding to V1 in paper
elseif flag == 2
    VM = u*M;   % corresponding to V2 in paper
else
    VM = V2*S/(K2+S);   % corresponding to V3 in paper
end

% ODE equations for P, L, R, M, B, and S, respectively.
rhs1 = G*s - (aL+aR+aM)*P;  % P
rhs2 = aL*P - kL*L;         % L
rhs3 = aR*P - kR*R;         % R
rhs4 = eP*aM*P - kM*M - rM*M + eS*VM;   % M
rhs5 = V1*S/(K1+S) - kB*B -rB*B;        % B
rhs6 = (1-rL)*kL*L + kR*R + kM*M + kB*B - VM - V1*S/(K1+S); % S

ydot = [rhs1; rhs2; rhs3; rhs4; rhs5; rhs6];