function sol = arab

u = [0.15,...%  u(1) - TFmin    Level of TF below which PCD does not occur.
     0.5,...%   u(2) - TFmax    Level of TF above which PCD no longer depends on TF levels.
     0.17,...%  u(3) - v11      Flux representing the maximal rate of TF-dependent PCD in the absence of negative regulation.
     0.2,...%   u(4) - K11      Level of TF (relative to a maximum of 1.0) when the effect of TF on PCD is half-maximal.
     2.5,...%   u(5) - O2- min  Level of superoxide below which direct induction of PCD by high superoxide levels does not happen.
     2,...%     u(6) -          DELETED FROM MODEL
     0.11,...%  u(7) - c12      Proportionality constant relating increases in superoxide levels (above O2?min) to changes in PCD.    
     3,...%     u(8) - K13      Level of O2? when direct induction of PCD is half-maximal.
     1,...%     u(9) - c11      Proportionality constant relating increases in RNApcd to changes in PCD.  The expression including this constant accounts for NPR1-dependent negative regulation of PCD.
     1,...%     u(10) - NPR1    Gene dosage of NPR1.  Value is one for wild type, zero for npr1 mutant.
     1,...%     u(11) - K12     Level of RNApcd when NPR1-dependent negative regulation of PCD is half-maximal.
     1,...%     u(12) - t11     Delay between changes in [SA] and consequent potentiation  of TF action.
     6,...%     u(13) - K14     Level of SA when potentiation of TF action is half-maximal. 
     50,...%    u(14) - SAmax   DELETED FROM MODEL
     .35,...%   u(15) - v21     Flux representing the basal rate of SA biosynthesis from CM.
     0.01,...%  u(16) - K21     Level of CM when basal SA biosynthesis from CM is half-maximal.  This value doesn’t have a major effect in simulations presented in this paper, since it was taken as being much lower than [CM].
     1,...%     u(17) - CM      Level of CM.  Taken as large and unchanging under the conditions used in simulations.
     1,...%     u(18) - NDR1    Gene dosage of NDR1.  Value is one for wild type, zero for ndr1 mutant.
     8,...%     u(19) - c21     Proportionality constant relating effect of NDR1 gene dosage to extent of superoxide-induced increases in SA biosynthesis from CM.
     3,...%     u(20) - K22     Level of superoxide when superoxide-induced increase in SA biosynthesis from CM is half-maximal. 
     1,...%     u(21) - c22     Proportionality constant relating level of avr-R to changes in SA biosynthesis from CM.  This term accounts for superoxide-independent, NDR1-independent upregulation of this pathway. 
     0.2,...%   u(22) - K23     Level of avr-R when superoxide-independent, NDR1-independent upregulation of SA biosynthesis from CM is half-maximal.
     0.27,...%  u(23) - v22     Flux representing the basal rate of SA biosynthesis from the alternative pathway.
     0.01,...%  u(24) - K24     Level of Alt (hypothetical precursor in alternative SA biosynthesis pathway, likely to be phenylalanine) when basal SA biosynthesis from Alt is half-maximal.  This value doesn’t have a major effect in simulations presented in this paper, since it was taken as being much lower than [Alt].
     1,...%     u(25) - Alt     Level of Alt.  Taken as large and unchanging under the conditions used in simulations.
     8,...%     u(26) - c23     Proportionality constant relating effect of NDR1 gene dosage to extent of superoxide-induced increases in SA biosynthesis from Alt.
     3,...%     u(27) - K25     Level of superoxide when superoxide-induced increase in SA biosynthesis from Alt is half-maximal. 
     1,...%     u(28) - c24     Proportionality constant relating level of avr-R to changes in SA biosynthesis from Alt.  This term accounts for superoxide-independent, NDR1-independent upregulation of this pathway.
     0.2,...%   u(29) - K26     Level of avr-R when upregulation of superoxide-independent, NDR1-independent SA biosynthesis from Alt is half-maximal.
     1,...%     u(30) - c25     Proportionality constant relating level of hydrogen peroxide to changes in SA biosynthesis from Alt.  
     3,...%     u(31) - K27     Level of hydrogen peroxide when hydrogen peroxide-dependent upregulation of SA biosynthesis from Alt is half-maximal.
     1.2,...%   u(32) -         DELETED FROM MODEL Rate constant modifying rate of SA biosynthesis due to NPR1-dependent negative regulation
     1.5,...%   u(33) -         DELETED FROM MODEL Equilibrium constant governing impact of c24 at varying levels of SA.
     0.25,...%  u(34) - k21     First-order rate constant for SA degradation. 
     0.3,...%   u(35) - v31     Flux representing the basal rate of superoxide production. 
     0.1,...%   u(36) - c31     Proportionality constant relating level of avr-R to changes in superoxide production (minor).
     0.15,...%  u(37) - K31     Level of avr-R when upregulation of superoxide production is half-maximal.
     5,...%     u(38) - c32     Proportionality constant relating changes in PCD to consequent changes in superoxide production.  
     0.000025,...%   u(39) - k31 Bimolecular rate constant for uncatalyzed production of hydrogen peroxide from superoxide.
     0.10,...%  u(40) - k32     Rate constant for SOD-catalyzed production of hydrogen peroxide from superoxide
     0.04,...%  u(41) - k33     First-order degradation rate constant for destruction of superoxide by all alternative oxidants and reductants.  This is a lumped term representing all such pathways.
     0.09,...%  u(42) - k41     First-order degradation rate constant for destruction of superoxide.  This is a lumped term representing all apoplastic catalase and peroxidase activities.
     0.3,...%   u(43) - v51     Flux representing the basal rate of SOD production.  Basal specific activity of SOD inuninoculated wild type plants was defined as 1.
     1,...%     u(44) -         DELETED FROM MODEL 
     1,...%     u(45) -         DELETED FROM MODEL
     0.3,...%   u(46) - c51     Proportionality constant relating level of RNAsod to changes in SOD production.  The term containing this constant relates changes in SA to changes in apoplastic SOD specific activity. The assumption was made that transcription of the SOD gene is slower than all subsequent steps in the process of changing SOD activity.
     1,...%     u(47) - LSD1    Gene dosage of LSD1.  Value is one for wild type, zero for lsd1 mutant.
     6,...%     u(48) -         DELETED FROM MODEL
     0.1,...%   u(49) - c52     Proportionality constant relating level of avr • R to changes in SOD production. 
     1,...%     u(50) - K51     Level of avr-R when upregulation of SOD production is half-maximal.
     0.3,...%   u(51) - k51     First-order degradation rate constant for SOD.
     .23,...%   u(52) - v61     Flux representing the basal rate of RNApcd production.
     20,...%    u(53) - c61     Proportionality constant relating effect of NPR1 gene dosage to extent of SA-induced increase in RNApcd production.
     18.8,...%  u(54) - K61     Level of SA when SA-induced increase in RNApcd production is half-maximal.
     3.5,...%   u(55) - t61     Delay between changes in [SA] and production of RNApcd. 
     0.23,...%  u(56) - k61     First-order degradation rate constant for RNApcd. 
     0.8,...%   u(57) - v71     Flux representing the basal rate of RNAsod production.
     10,...%    u(58) - c71     Proportionality constant relating effect of NPR1 gene dosage to extent of SA-induced increase in RNAsod production.
     12,...%    u(59) - K71     Level of SA when SA-induced increase in RNAsod production is half-maximal.
     1,...%     u(60) - t71     Delay between changes in [SA] and production of RNAsod.
     0.8,...%   u(61) - k71     First-order degradation rate constant for RNAsod. 
     0.055,...% u(62) -         DELETED FROM MODEL Basal cellular production rate of the RNA corresponding to the putative gene that is activated by SA, where the gene product functions to negatively autoregulate SA production.  
     14,...%    u(63) -         DELETED FROM MODEL Effect of SA signaling via NPR1 on v91.  
     8,...%     u(64) -         DELETED FROM MODEL Effect of SA signaling via NPR1 on v91.  
     2,...%     u(65) -         DELETED FROM MODEL Delay between changes in [SA] and production of RNASA
     0.15,...%  u(66) -         DELETED FROM MODEL Degradation rate constant for RNASA.  See above.
     0.7,...%   u(67) - k81     Proportionality constant allowing a sigmoid function to be altered to simulate response to different avr genes.  This was set to 0.7 for avrB and to 0.35  for avrRpt2.
     2.5,...%   u(68) - t81     Beginning of sigmoid in equation for avr-R.
     1,...%     u(69) - t82     Time between beginning and midpoint of sigmoid in equation for avr-R.  
     0.15,...%  u(70) - k82     First-order degradation rate constant for avr-R.  Note: the same variable was used in the code for k92, which is the first-order degradation rate constant for TF.  The intention was to make TF levels closely parallel avr-R levels.  
     1,...%     u(71) - k91     Proportionality constant relating level of avr-R to level of TF.
     1]; %      u(72) - t91     Delay between changes in level of avr-R  and action of TF
              
tfinal = 8;   
tspan = [0 tfinal];
lags = [1, 2, 3.5]; 
yhist = [.1,...%PCD
         1,...%SA
         1,...%Superoxide
         1,...%DELETED FROM MODEL
         1,...%Hydrogen Peroxide
         1,...%Apoplastic SOD Activity
         1,...%RNApcd
         1,...%RNAsod
         1,...%RNAsa DELETED FROM MODEL
         .01,...%avr-R
         .01];...%TF
opts = [];

tic

DataFile = [];

for i=1:4
    if i == 1 % Wild Type
        sol = dde15s(@arabidopsis,lags,yhist,tspan,opts,u);
        tout = sol.x; 
        xint = 0:0.1:8;
        tout = xint;
        sol = deval(sol, xint);
        timewild = tout;
        sol(1,:) = 1.3*sol(1,:);
        wild = sol;
        DataFile = [DataFile, transpose(tout), transpose(sol)];
    elseif i == 2 % npr1 mutant
        u(10) = 0;
        sol = dde15s(@arabidopsis,lags,yhist,tspan,opts,u);
        tout = sol.x; 
        xint = 0:0.1:8;
        tout = xint;
        sol = deval(sol, xint);
        timenpr1 = tout;
        sol(1,:) = 1.3*sol(1,:);
        npr1 = sol;
        DataFile = [DataFile, transpose(tout), transpose(sol)];
    elseif i == 3 % ndr1 mutant
        u(10) = 1;
        u(18) = 0;
        sol = dde15s(@arabidopsis,lags,yhist,tspan,opts,u);
        tout = sol.x; 
        xint = 0:0.1:8;
        tout = xint;
        sol = deval(sol, xint);
        timendr1 = tout;
        sol(1,:) = 1.3*sol(1,:);
        ndr1 = sol;
        DataFile = [DataFile, transpose(tout), transpose(sol)];
    end
end
    
toc

subplot(3,3,1); plot(timewild, wild(1,:),'k'); title('PCD in DC3000 avrB');
hold on
subplot(3,3,1); plot(timenpr1, npr1(1,:),'b'); 
subplot(3,3,1); plot(timendr1, ndr1(1,:),'r'); 

subplot(3,3,2); plot(timewild, wild(2,:),'k'); title('Salicylate Levels');
hold on
subplot(3,3,2); plot(timenpr1, npr1(2,:),'b'); 
subplot(3,3,2); plot(timendr1, ndr1(2,:),'r'); 

subplot(3,3,3); plot(timewild, wild(3,:),'k'); title('Superoxide');
hold on
subplot(3,3,3); plot(timenpr1, npr1(3,:),'b'); 
subplot(3,3,3); plot(timendr1, ndr1(3,:),'r'); 

subplot(3,3,4); plot(timewild, wild(5,:),'k'); title('Hydrogen Peroxide');
hold on
subplot(3,3,4); plot(timenpr1, npr1(5,:),'b'); 
subplot(3,3,4); plot(timendr1, ndr1(5,:),'r'); 

subplot(3,3,5); plot(timewild, wild(6,:),'k'); title('Apoplastic SOD Activity');
hold on
subplot(3,3,5); plot(timenpr1, npr1(6,:),'b'); 
subplot(3,3,5); plot(timendr1, ndr1(6,:),'r'); 

subplot(3,3,6); plot(timewild, wild(7,:),'k'); title('RNApcd');
hold on
subplot(3,3,6); plot(timenpr1, npr1(7,:),'b'); 
subplot(3,3,6); plot(timendr1, ndr1(7,:),'r'); 

subplot(3,3,7); plot(timewild, wild(8,:),'k'); title('RNAsod');
hold on
subplot(3,3,7); plot(timenpr1, npr1(8,:),'b'); 
subplot(3,3,7); plot(timendr1, ndr1(8,:),'r'); 

subplot(3,3,8); plot(timewild, wild(10,:),'k'); title('avr-R');
hold on
subplot(3,3,8); plot(timenpr1, npr1(10,:),'b'); 
subplot(3,3,8); plot(timendr1, ndr1(10,:),'r'); 

subplot(3,3,9); plot(timewild, wild(11,:),'k'); title('Triggering Factors');
hold on
subplot(3,3,9); plot(timenpr1, npr1(11,:),'b'); 
subplot(3,3,9); plot(timendr1, ndr1(11,:),'r'); 

filename = '';
name = ['avrBplot-' date];
filename = sprintf(name);
saveas(gcf, filename);
saveas(gcf, filename, 'tiff');
saveas(gcf, filename, 'jpg');
saveas(gcf, filename, 'emf');
name = ['avrBplot-' date '.txt']; % Data file for plotting in other software
filename = sprintf(name);
dlmwrite(filename, DataFile, '\t');


%=============================================
function yp = arabidopsis(t,y,Z,u)
% variables
%	y(1)  = Programmed Cell Death [PCD]
%	y(2)  = Salicylate Levels [SA]
%	y(3)  = Superoxide [O2-]
%   y(4)  = DELETED FROM MODEL
%   y(5)  = Hydrogen Peroxide [H2O2]
%   y(6)  = Apoplastic SOD activity [SOD]
%	y(7)  = [RNApcd]
%	y(8)  = [RNAsod]
%	y(9)  = [RNAsa] DELETED FROM MODEL
%   y(10)  = [avr-R]
%   y(11)  = Triggering Factor [TR]

%  Time Delay Values
ylag1 = Z(:,1);
ylag2 = Z(:,2);
ylag3 = Z(:,3);
yp = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];

%  PCD
if (y(11)<u(1)) & (y(3)<=u(5))
    yp(1) = 0;
elseif ((y(3)<=u(5)) & (y(11)>=u(1)) & (y(11)<=u(2)) & ( y(11) > (u(1) + (u(2)-u(1))/(1+ylag1(2)/u(13)))))
    R_TF = u(3)*y(11)/(y(11)+u(4));
    %M_O2 = u(7)*y(3)/(y(3)+u(8)); %DELETED FROM MODEL
    N_NPR1_PCD = u(9)/(u(11)+y(7));
    yp(1) = R_TF*N_NPR1_PCD;% + M_O2;
elseif ((y(3)>u(5)) & (y(11)>=u(1)) & (y(11)<=u(2)) & ( y(11) > (u(1) + (u(2)-u(1))/(1+ylag1(2)/u(13)))))
    R_TF = u(3)*y(11)/(y(11)+u(4));
    M_O2 = u(7)*y(3)/(y(3)+u(8));
    N_NPR1_PCD = u(9)/(u(11)+y(7));
    yp(1) = R_TF*N_NPR1_PCD + 2*M_O2;
elseif ((y(11)>u(2)) | (y(2)>u(14)))
    R_TF_MAX = u(3)*u(2)/(u(2)+u(4));
    M_O2 = u(7)*y(3)/(y(3)+u(8));
    N_NPR1_PCD = u(9)/(u(11)+y(7));
    yp(1) = R_TF_MAX*N_NPR1_PCD+M_O2;
else
    yp(1) = 0;
end

% Salicylic Acid

if (t<=0.5)
    yp(2) = 0;
else
    M_CM = u(15)*u(17)/(u(17)+u(16));
    R_IO2 = 1+u(18)*u(19)*y(4)/(y(4)+u(20));
    R_avr_R = u(21)*y(10)/(y(10)+u(22));
    M_Alt =  u(23)*u(25)/(u(25)+u(24));
    R_IO2_Alt = 1+u(18)*(u(26)*y(4)/(y(4)+u(27)));
    R_avr_R_Alt = u(28)*y(10)/(y(10)+u(29));
    R_H2O2 = u(30)*y(5)/(y(5)+u(31));
    %N_NPR1 = (u(32)*u(10)/(1+y(9)/u(33)))^u(10); % DELETED FROM MODEL
    D_SA = u(34)*y(2);
    yp(2) = (M_CM*(R_IO2 + R_avr_R) + M_Alt*(R_IO2_Alt + R_avr_R_Alt + R_H2O2)) - D_SA;% *N_NPR1 % DELETED FROM MODEL
end

% Superoxide

if t<=1
    yp(3) = 0;
else
    M_avrR = u(36)*y(10)/(y(10)+u(37));
    M_PCD = u(38)*yp(1);
    M_O22 = u(39)*y(3)*y(3);
    M_SOD = u(40)*(y(6)^1.5)*(y(3)^1.5)/(1+y(3));
    D_O2 = u(41)*y(3);
    yp(3) = u(35) + M_avrR + M_PCD - M_O22 - M_SOD - D_O2;
end

%DELETED FROM MODEL
yp(4) = y(3); % DELETED FROM MODEL

%H2O2
M_O22 = u(39)*y(3)*y(3);
M_SOD = u(40)*(y(6)^1.5)*(y(3)^1.5)/(1+y(3));

yp(5) = M_O22 + M_SOD - u(42)*y(5);

%Superoxide dismutase
yp(6) = u(43) + u(46)*(y(8)-1) + u(49)*y(10)/(y(10)+u(50)) - u(51)*y(6);

%RNA_PCD SA negative regulation of PCD via NPR1
yp(7) = u(52) + u(10)*((u(53)*(ylag3(2)-1))/(u(54)+(ylag3(2)-1))) - u(56)*y(7);

%RNA_SOD
yp(8) = u(57) + u(47)*((u(58)*(ylag1(2)-1))/(u(59)+(ylag1(2)-1))) - u(61)*y(8);

%RNA_SA - DELETED FROM MODEL
yp(9) = u(62) + u(10)*((u(63)*ylag2(2))/(u(64)+ylag2(2))) - u(66)*y(9); % DELETED FROM MODEL

%avr-R gene interaction
yp(10) = 0.5*u(67)*((u(68)*exp(-u(68)*(t-u(69))))/(1+exp(-u(68)*(t-u(69))))^2)-u(70)*y(10);

%TF (Triggering factors leading to PCD)
if t<=2
    yp(11) = 0;
elseif t>2 
    yp(11) = (u(67)*((u(68)*exp(-u(68)*((t-2)-u(69))))/(1+exp(-u(68)*((t-2)-u(69))))^2)-u(70)*y(11));
end