% 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 pertubate parameter and derive solutions analytically
clear;

% Parameters
G = 1000;       % Net primary productivity
aL0 = 0.05;     % Leaf litter allocation
aR0 = 0.03;     % Root allocation
aM0 = 0.02;     % Mycorrhizal allocation
kL0 = 0.5;      % Litter turnover
kR0 = 0.3;      % Root turnover
kM0 = 3;        % Mycorrhizal turnover
kB0 = 3.75;     % Microbial turnover
rL0 = 0.9;      % Fraction of litter respired
rM0 = 1;        % Mycorrhizal respiration
rB0 = 5;        % Microbial respiration
eP = 1;         % Plant acquisition efficiency
eS0 = 0.3;      % Soil acquisition efficienty
u0 = 1.64;      % Specific mycorrhizal C uptake
V10 = 1600;     % Maximum microbial C uptake
K10 = 5000;     % Microbial half-saturation
V20 = 200;      % Maximum mycorrhizal C uptake
K20 = 5000;     % Mycorrhizal half-saturation
s = 1;          % scales NPP

%% Pertubation experiments
Trecord = zeros(1000,28);
for inde = 1:1000
    s = rand
    aL = aL0*(0.9+0.2*rand);
    aM = aM0*(0.9+0.2*rand);
    aR = aR0*(0.9+0.2*rand); 
    eS = eS0*(0.9+0.2*rand);
    kL = kL0*(0.9+0.2*rand);
    kR = kR0*(0.9+0.2*rand);
    kM = kM0*(0.9+0.2*rand);
    kB = kB0*(0.9+0.2*rand);
    rB = rB0*(0.9+0.2*rand);
    rM = rM0*(0.9+0.2*rand);
    rL = rL0*(0.9+0.2*rand);
    K1 = K10*(0.9+0.2*rand);
    K2 = K20*(0.9+0.2*rand);
    u = u0*(0.9+0.2*rand);
    V1 = V10*(0.9+0.2*rand);
    V2 = V20*(0.9+0.2*rand);

% Steady state P, L, and R respectively.     
P_star = G*s/(aL + aM + aR);
L_star = aL/kL*P_star;
R_star = aR/kR*P_star;

% Steady state M, B, and S in case of V1
M_star1 = eP*aM*P_star/(kM+rM);
B_star1 = (G*s-rL*kL*L_star-rM*M_star1-(1-eP)*aM*P_star)/rB;
S_star1 = (kB+rB)*B_star1*K1/(V1-(kB+rB)*B_star1);

% Steady state M, B, and S in case of V2
M_star2 = eP*aM*P_star/(kM+rM-eS*u);
B_star2 = (G*s-rL*kL*L_star-rM*M_star2-(1-eP)*aM*P_star - (1-eS)*u*M_star2)/rB;
S_star2 = (kB+rB)*B_star2*K1/(V1-(kB+rB)*B_star2);

% Steady state M, B, and S in case of V3
syms B M S
B_hat = solve(V1*S/(K1+S)==kB*B+rB*B,B);
M_hat = solve(eP*aM*P_star+eS*V2*S/(K2+S)==kM*M+rM*M,M);
S_hat = solve(G*s-rL*kL*L_star-rM*M_hat-(1-eP)*aM*P_star-(1-eS)*V2*S/(K2+S)==rB*B_hat,S);
S_star3 = eval(S_hat(1));
M_star3 = subs(M_hat,S,S_star3);
B_star3 = subs(B_hat,S,S_star3);

% For each uptake strategy, calculate key C flux
MgetfromP1 = eP*aM*P_star;      % Mycorrhizae get C from plant
MgetfromS1 = 0;                 % Mycorrhizae get C from soil
SdecombyB1 = V1*S_star1./(K1+S_star1);  % Soil C decomposed by bacteria
SdecombyM1 = 0;                         % Soil C decomposed by mycorrhizae
MfluxtoS1 = kM*M_star1;         % Mycorrhizae turnover C to soil

MgetfromP2 = eP*aM*P_star;
MgetfromS2 = eS*u*M_star2;
SdecombyB2 = V1*S_star2./(K1+S_star2);
SdecombyM2 = u*M_star2;
MfluxtoS2 = kM*M_star2;

MgetfromP3 = eP*aM*P_star;
MgetfromS3 = eS*V2*S_star3./(K2+S_star3);
SdecombyB3 = V1*S_star3./(K1+S_star3);
SdecombyM3 = V2*S_star3./(K2+S_star3);
MfluxtoS3 = kM*M_star3;

% Record all the data for latter analysis
Trecord(inde,:) = [P_star L_star R_star M_star1 M_star2 M_star3 ...
    B_star1 B_star2 B_star3 S_star1 S_star2 S_star3 ...
    MgetfromP1 MgetfromS1 SdecombyB1 SdecombyM1 MfluxtoS1 ...
    MgetfromP2 MgetfromS2 SdecombyB2 SdecombyM2 MfluxtoS2 ...
    MgetfromP3 MgetfromS3 SdecombyB3 SdecombyM3 MfluxtoS3 s];

end

%%
% Retrieve data from the dataset.
load Trecord1
Pstar = Trecord(:,1);
Lstar = Trecord(:,2);
Rstar = Trecord(:,3);
Mstar1 = Trecord(:,4);
Mstar2 = Trecord(:,5);
Mstar3 = Trecord(:,6);
Bstar1 = Trecord(:,7);
Bstar2 = Trecord(:,8);
Bstar3 = Trecord(:,9);
Sstar1 = Trecord(:,10);
Sstar2 = Trecord(:,11);
Sstar3 = Trecord(:,12);

MgetfromP1 = Trecord(:,13);
MgetfromS1 = Trecord(:,14);
SdecombyB1 = Trecord(:,15);
SdecombyM1 = Trecord(:,16);
MfluxtoS1 = Trecord(:,17);

MgetfromP2 = Trecord(:,18);
MgetfromS2 = Trecord(:,19);
SdecombyB2 = Trecord(:,20);
SdecombyM2 = Trecord(:,21);
MfluxtoS2 = Trecord(:,22);

MgetfromP3 = Trecord(:,23);
MgetfromS3 = Trecord(:,24);
SdecombyB3 = Trecord(:,25);
SdecombyM3 = Trecord(:,26);
MfluxtoS3 = Trecord(:,27);

ss = Trecord(:,28);


%% percent change in soil C compare to V1

xs = ss*100;

p2 = (Sstar2-Sstar1)./Sstar1*100;
p3 = (Sstar3-Sstar1)./Sstar1*100;
pp2 = p2;
pp3 = p3;

xs(p2<-80 | p3<-80)= [];
pp2(p2<-80 | p3<-80)=[];
pp3(p2<-80 | p3<-80)=[];
inde = find(p2<-80 | p3<-80);

poly2 = polyfit(xs,pp2,2);
poly3 = polyfit(xs,pp3,2);

x = [0:100];
f2 = polyval(poly2,x);
f3 = polyval(poly3,x);
figure
hold on
plot(x,f2,'k-.')
plot(x,f3,'k')
plot(xs,pp2,'ko','MarkerSize',4)
plot(xs,pp3,'ko','MarkerFaceColor',[0 0 0],'MarkerSize',4)
plot(x,f2,'k-.','Linewidth',1.5)
plot(x,f3,'k','Linewidth',1.5)
xlabel('Level of plant productivity (%)','fontsize',14)
ylabel('Percent change in soil C (%)','fontsize',14)
legend('Constant uptake','Variable uptake','Location',[0.2 0.2 0.2 0.1])
legend('boxoff')
set(gca,'Fontname','Arial')

%% proportion of Mycorrhizae C uptake  
figure
hold on

x = ss*100;
y = MgetfromS3./(MgetfromP3+MgetfromS3)*100;

x(inde) = [];
y(inde) = [];
plot(x,y,'ko','MarkerFaceColor',[0 0 0],'MarkerSize',4)
xlabel('Level of plant productivity (%)','fontsize',14)
ylabel('Percent of mycorrhizae C uptake from soil (%)','fontsize',14)
p = polyfit(x,y,2);
x1 = [0:100];
f = polyval(p,x1);
hold on
plot(x1,f,'k','linewidth',1.5)
set(gca,'Fontname','Arial')


%% percent of soil C transfferred to mycorrhizae 

figure
hold on
x=ss*100;
y = SdecombyM3./(SdecombyB3+SdecombyM3)*100;

x(inde) = [];
y(inde) = [];
plot(x,y,'ko','MarkerFaceColor',[0 0 0],'MarkerSize',4)
xlabel('Level of plant productivity (%)','fontsize',14)
ylabel('Percent of soil C transferred to mycorrhizae (%)','fontsize',14)
p = polyfit(x,y,2);
x1 = [0:100];
f = polyval(p,x1);
hold on
plot(x1,f,'k','linewidth',1.5)
set(gca,'Fontname','Arial')