Team:Virginia/Model/Files/denitIntermediates.nb

ClearAll["Global`*"]

rates = {R1[t] == 
    Subscript[r, COD]*X[t]*
     Subscript[S, S][t]/(Subscript[k, S] + Subscript[S, S][t])*
     Subscript[S, Mox][t]/(Subscript[k, Mox] + Subscript[S, Mox][t]), 
   R2[t] == 
    Subscript[r, NO3]*X[t]*
     Subscript[S, NO3][t]/(Subscript[k, NO3] + Subscript[S, NO3][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred1] + Subscript[S, Mred][t]), 
   R3[t] == 
    Subscript[r, NO2]*X[t]*
     Subscript[S, NO2][t]/(Subscript[k, NO2] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred2] + Subscript[S, Mred][t]), 
   R4[t] == 
    Subscript[r, NO]*X[t]*
     Subscript[S, NO][t]/(Subscript[k, NO] + Subscript[S, NO][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred3] + Subscript[S, Mred][t]), 
   R5[t] == 
    Subscript[r, N2O]*X[t]*
     Subscript[S, N2O][t]/(Subscript[k, N2O] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred4] + Subscript[S, Mred][t])};

conc = {Subscript[S, NO3]'[t] == -R2[t], 
   Subscript[S, NO2]'[t] == R2[t] - R3[t], 
   Subscript[S, NO]'[t] == R3[t] - R4[t], 
   Subscript[S, N2O]'[t] == 0.5 R4[t] - R5[t], 
   Subscript[S, N2]'[t] == R5[t], Subscript[S, S]'[t] == -R1[t], 
   Subscript[S, Mox]'[t] == -(1 - Subscript[Y, H]) R1[t] + R2[t] + 
     0.5 R3[t] + 0.5 R4[t] + R5[t], X'[t] == Subscript[Y, H]*R1[t], 
   Subscript[S, Mox][t] + Subscript[S, Mred][t] == Subscript[c, tot]};

(*Case 1:*)
Subscript[r, COD] = 0.064; Subscript[r, NO3] = 0.045; 
Subscript[r, NO2] = 0.059; Subscript[r, NO] = 0.56; 
Subscript[r, N2O] = 0.23;

Subscript[k, S] = 0.1; Subscript[k, NO3] = 0.018; 
Subscript[k, NO2] = 0.0041;
Subscript[k, NO] = 0.000011;
Subscript[k, N2O] = 0.0025;

Subscript[k, Mox] = 0.0001;
Subscript[k, Mred1] = 0.0015;
Subscript[k, Mred2] = 0.00058;
Subscript[k, Mred3] = 0.000010;
Subscript[k, Mred4] = 0.00024;
Subscript[Y, H] = 0.6;
Subscript[c, tot] = 0.01;

initc = {WhenEvent[t == 0.5, Subscript[S, NO2][t] -> 5/14], 
   Subscript[S, NO3][0] == 1, Subscript[S, NO2][0] == 0, 
   Subscript[S, NO][0] == 0, Subscript[S, N2O][0] == 0, 
   Subscript[S, N2][0] == 0, Subscript[S, S][0] == 170/14, 
   Subscript[S, Mox][0] == 0.005, X[0] == 70};

Subscript[t, 0] = 1.5;

{SNO3, SNO2, SNO, SN2O, SN2, COD, SMox, SMred, SX} = 
     NDSolveValue[
   Join[rates, conc, initc], {Subscript[S, NO3], Subscript[S, NO2], 
         Subscript[S, NO], Subscript[S, N2O], Subscript[S, N2], 
    Subscript[S, S], 
         Subscript[S, Mox], Subscript[S, Mred], X}, {t, 0., 
    Subscript[t, 0]} ,
       Method -> {"TimeIntegration" -> "StateSpace"}]; 

Plot[{SNO3[t], SNO2[t], SN2O[t]}, {t, 0., Subscript[t, 0]}, 
 PlotRange -> {0, 3}, 
   PlotStyle -> {{Red, Thick}, {Blue, Thick}, {Brown, Dashed, 
    Thick}}]
Plot[SX[t], {t, 0, Subscript[t, 0]}, PlotRange -> Full, 
 PlotStyle -> {Red, Thick}]
Plot[{SMox[t], SMred[t]}, {t, 0, Subscript[t, 0]}]

(* Nitrification addition and parameter sensitivity *)
\
ClearAll["Global`*"]

rates = {R1[t] == 
    Subscript[r, COD]*X[t]*
     Subscript[S, S][t]/(Subscript[k, S] + Subscript[S, S][t])*
     Subscript[S, Mox][t]/(Subscript[k, Mox] + Subscript[S, Mox][t]), 
   R2[t] == 
    Subscript[r, NO3]*X[t]*
     Subscript[S, NO3][t]/(Subscript[k, NO3] + Subscript[S, NO3][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred1] + Subscript[S, Mred][t]), 
   R3[t] == 
    Subscript[r, NO2]*X[t]*
     Subscript[S, NO2][t]/(Subscript[k, NO2] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred2] + Subscript[S, Mred][t]), 
   R4[t] == 
    Subscript[r, NO]*X[t]*
     Subscript[S, NO][t]/(Subscript[k, NO] + Subscript[S, NO][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred3] + Subscript[S, Mred][t]), 
   R5[t] == 
    Subscript[r, N2O]*X[t]*
     Subscript[S, N2O][t]/(Subscript[k, N2O] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred4] + Subscript[S, Mred][t]),
   RN1[t] == 
    Subscript[rN, NH3]*X[t]*
     Subscript[S, NH3][t]/(Subscript[kN, NH3] + Subscript[S, NH3][t])*
     Subscript[S, Mox][
       t]/(Subscript[kN, Mox1] + Subscript[S, Mox][t])};

conc = {Subscript[S, NH3]'[t] == -RN1[t], 
   Subscript[S, NO3]'[t] == -R2[t], 
   Subscript[S, NO2]'[t] == R2[t] - R3[t], 
   Subscript[S, NO]'[t] == R3[t] - R4[t] + RN1[t], 
   Subscript[S, N2O]'[t] == 0.5 R4[t] - R5[t], 
   Subscript[S, N2]'[t] == R5[t], Subscript[S, S]'[t] == -R1[t], 
   Subscript[S, Mox]'[t] == -(1 - Subscript[Y, H]) R1[t] + R2[t] + 
     0.5 R3[t] + 0.5 R4[t] + R5[t] - RN1[t], 
   X'[t] == Subscript[Y, H]*R1[t], 
   Subscript[S, Mox][t] + Subscript[S, Mred][t] == Subscript[c, tot]};

Subscript[r, COD] = 0.064; Subscript[r, NO3] = 0.045; 
Subscript[r, NO2] = 0.059; Subscript[r, NO] = 0.56; 
Subscript[r, N2O] = 0.23;

Subscript[k, S] = 0.1; Subscript[k, NO3] = 0.018; 
Subscript[k, NO2] = 0.0041; Subscript[k, NO] = 0.000011; 
Subscript[k, N2O] = 0.0025;

Subscript[k, Mox] = 0.0001;
Subscript[k, Mred1] = 0.0015;
Subscript[k, Mred2] = 0.00058;
Subscript[k, Mred3] = 0.000010;
Subscript[k, Mred4] = 0.00024;
Subscript[Y, H] = 0.6;
Subscript[c, tot] = 0.01;

## ## ## CASE 1 + BIOREACTOR - 
 pick only one of these, comment out the others 
initc = {Subscript[S, NH3][0] == 480.3/14, 
   WhenEvent[t == 0.5, Subscript[S, NO2][t] -> 5/14], 
   Subscript[S, NO3][0] == 1, Subscript[S, NO2][0] == 0, 
   Subscript[S, NO][0] == 0, Subscript[S, N2O][0] == 0, 
   Subscript[S, N2][0] == 0, Subscript[S, S][0] == 170/14, 
   Subscript[S, Mox][0] == 0.005, X[0] == 70};
## ## ## Simulation 1 : no intermediates 
initc = {Subscript[S, NH3][0] == 480.3/14, 
   Subscript[S, NO3][0] == 0.0001, Subscript[S, NO2][0] == 0.0001, 
   Subscript[S, NO][0] == 0.0001, Subscript[S, N2O][0] == 0.0001, 
   Subscript[S, N2][0] == 0.0001, Subscript[S, S][0] == 170/14, 
   Subscript[S, Mox][0] == 0.005, X[0] == 70};
## ## ## Simulation 2 : 
   initc = {Subscript[S, NH3][0] == 480.3/14, 
   WhenEvent[t == 0.5, Subscript[S, NO2][t] -> 5/14], 
   Subscript[S, NO3][0] == 1, Subscript[S, NO2][0] == 0, 
   Subscript[S, NO][0] == 0, Subscript[S, N2O][0] == 0, 
   Subscript[S, N2][0] == 0, Subscript[S, S][0] == 170/14, 
   Subscript[S, Mox][0] == 0.005, X[0] == 70};
## ## ## Final simulation
initc = {Subscript[S, NH3][0] == 480.3/14, 
   WhenEvent[t == 0.5, Subscript[S, NO2][t] -> 5/14], 
   Subscript[S, NO3][0] == 1, Subscript[S, NO2][0] == 0, 
   Subscript[S, NO][0] == 0, Subscript[S, N2O][0] == 0, 
   Subscript[S, N2][0] == 0, Subscript[S, S][0] == 170/14, 
   Subscript[S, Mox][0] == 0.005, X[0] == 70};

Subscript[t, 0] = 1.5;

params = {Subscript[rN, NH3], Subscript[kN, NH3], Subscript[kN, Mox1]};

sol = ParametricNDSolveValue[
  Join[rates, conc, initc], {Subscript[S, NH3], Subscript[S, NO3], 
   Subscript[S, NO2], 
        Subscript[S, NO], Subscript[S, N2O], Subscript[S, N2], 
   Subscript[S, S], 
        Subscript[S, Mox], Subscript[S, Mred], X}, {t, 0., 
   Subscript[t, 0]}, params, 
      Method -> {"TimeIntegration" -> "StateSpace"}]

abscissae = Range[0., 1.5, 0.1];
ordinates = {Function[x, x/14]@{480.3, 440.9, 427.7, 388.1, 372, 
     357.2, 336.1, 347, 342.9, 342.6, 327.8, 311.6, 301.5, 298.6, 
     291.5, 281.9}}; 

g1 = Plot[{sol[0.7681335, 0.010363, 0.0341977][[1]][t], 
    sol[0.7681335, 0.010363, 0.0441977][[1]][t], 
    sol[0.7681335, 0.010363, 0.0241977][[1]][t], 
    sol[0.7681335, 0.010363, 0.0141977][[1]][t], 
    sol[0.7681335, 0.010363, 0.0541977][[1]][t]}, {t, 0, 1.5}, 
   PlotRange -> {0, 36}, 
   PlotLabels -> {"Best Fit", "+0.01", "-0.01", "-0.02", "+0.02"}, 
   AxesLabel -> {t, 
     "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NH\), \(3\)]]\)"}, 
   PlotLabel -> "\!\(\*SubscriptBox[\(k\), \(Mox, 1\)]\) Variation"];
g2 = Plot[{sol[0.7681335, 0.010363, 0.0341977][[1]][t], 
    sol[0.7681335, 1.100363, 0.0341977][[1]][t], 
    sol[0.7681335, 10.100363, 0.0341977][[1]][t], 
    sol[0.7681335, 100.110363, 0.0341977][[1]][t], 
    sol[0.7681335, 0, 0.0341977][[1]][t]}, {t, 0, 1.5}, 
   PlotRange -> {0, 36}, 
   PlotLabels -> {"Best Fit", "+1.0", "+10.0", "+100.0", "kN=0"}, 
   AxesLabel -> {t, 
     "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NH\), \(3\)]]\)"}, 
   PlotLabel -> 
    "\!\(\*SubscriptBox[\(k\), SubscriptBox[\(NH\), \(3\)]]\) \
Variation"];
g3 = Plot[{sol[0.7681335, 0.010363, 0.0341977][[1]][t], 
    sol[0.5681335, 0.010363, 0.0341977][[1]][t], 
    sol[0.3681335, 0.010363, 0.0341977][[1]][t], 
    sol[0.9681335, 0.010363, 0.0341977][[1]][t], 
    sol[1.1681335, 0.010363, 0.0341977][[1]][t], 
    sol[1.7681335, 0.010363, 0.0341977][[1]][t]}, {t, 0, 1.5}, 
   PlotRange -> {0, 36}, 
   PlotLabels -> {"Best Fit", "-0.2", "-0.4", "+0.2", "+0.4", "+1.0"},
    AxesLabel -> {t, 
     "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NH\), \(3\)]]\)"}, 
   PlotLabel -> 
    "\!\(\*SubscriptBox[\(r\), \(\*SubscriptBox[\(NH\), \(3\)], \
max\)]\) Variation"];
g4 = ListPlot[
   Table[{abscissae[[k]], ordinates[[1]][[k]]}, {k, 1, 16}], 
   PlotRange -> All];
{Show[g1, g4],
 Show[g2, g4],
 Show[g3, g4]}
Plot[{sol[0.7681335, 0.010363, 0.0341977][[6]][t], 
  sol[0, 0, 0][[6]][t]}, {t, 0, 1.5}, PlotRange -> All, 
 AxesLabel -> {t, 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(N\), \(2\)]]\)"}, 
 PlotStyle -> {Thick, Thick}, PlotLabel -> "Simulation 2"]


(* Parameter estimation *)
ClearAll["Global`*"];

abscissae = Range[0., 1.5, 0.1];
ordinates = {Function[x, x/14]@{480.3, 440.9, 427.7, 388.1, 372, 
     357.2, 336.1, 347, 342.9, 342.6, 327.8, 311.6, 301.5, 298.6, 
     291.5, 281.9}}; 

rates = {R1[t] == 
    Subscript[r, COD]*X[t]*
     Subscript[S, S][t]/(Subscript[k, S] + Subscript[S, S][t])*
     Subscript[S, Mox][t]/(Subscript[k, Mox] + Subscript[S, Mox][t]), 
   R2[t] == 
    Subscript[r, NO3]*X[t]*
     Subscript[S, NO3][t]/(Subscript[k, NO3] + Subscript[S, NO3][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred1] + Subscript[S, Mred][t]), 
   R3[t] == 
    Subscript[r, NO2]*X[t]*
     Subscript[S, NO2][t]/(Subscript[k, NO2] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred2] + Subscript[S, Mred][t]), 
   R4[t] == 
    Subscript[r, NO]*X[t]*
     Subscript[S, NO][t]/(Subscript[k, NO] + Subscript[S, NO][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred3] + Subscript[S, Mred][t]), 
   R5[t] == 
    Subscript[r, N2O]*X[t]*
     Subscript[S, N2O][t]/(Subscript[k, N2O] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred4] + Subscript[S, Mred][t]),
   RN1[t] == 
    Subscript[rN, NH3]*X[t]*
     Subscript[S, NH3][t]/(Subscript[kN, NH3] + Subscript[S, NH3][t])*
     Subscript[S, Mox][
       t]/(Subscript[kN, Mox1] + Subscript[S, Mox][t])};

conc = {Subscript[S, NH3]'[t] == -RN1[t], 
   Subscript[S, NO3]'[t] == -R2[t], 
   Subscript[S, NO2]'[t] == R2[t] - R3[t], 
   Subscript[S, NO]'[t] == R3[t] - R4[t] + RN1[t], 
   Subscript[S, N2O]'[t] == 0.5 R4[t] - R5[t], 
   Subscript[S, N2]'[t] == R5[t], Subscript[S, S]'[t] == -R1[t], 
   Subscript[S, Mox]'[t] == -(1 - Subscript[Y, H]) R1[t] + R2[t] + 
     0.5 R3[t] + 0.5 R4[t] + R5[t] - RN1[t], 
   X'[t] == Subscript[Y, H]*R1[t], 
   Subscript[S, Mox][t] + Subscript[S, Mred][t] == Subscript[c, tot]};

Subscript[r, COD] = 0.064; Subscript[r, NO3] = 0.045; 
Subscript[r, NO2] = 0.059; Subscript[r, NO] = 0.56; 
Subscript[r, N2O] = 0.23;

Subscript[k, S] = 0.1; Subscript[k, NO3] = 0.018; 
Subscript[k, NO2] = 0.0041; Subscript[k, NO] = 0.000011; 
Subscript[k, N2O] = 0.0025;

Subscript[k, Mox] = 0.0001;
Subscript[k, Mred1] = 0.0015;
Subscript[k, Mred2] = 0.00058;
Subscript[k, Mred3] = 0.000010;
Subscript[k, Mred4] = 0.00024;
Subscript[Y, H] = 0.6;
Subscript[c, tot] = 0.01;

initc = {Subscript[S, NH3][0] == 490/14, 
   WhenEvent[t == 0.5, Subscript[S, NO2][t] -> 5/14], 
   Subscript[S, NO3][0] == 1, Subscript[S, NO2][0] == 0, 
   Subscript[S, NO][0] == 0, Subscript[S, N2O][0] == 0, 
   Subscript[S, N2][0] == 0, Subscript[S, S][0] == 170/14, 
   Subscript[S, Mox][0] == 0.005, X[0] == 70};

Subscript[t, 0] = 1.5;

params = {Subscript[rN, NH3], Subscript[kN, NH3], Subscript[kN, Mox1]};

sol = ParametricNDSolveValue[
  Join[rates, conc, initc], {Subscript[S, NH3], Subscript[S, NO3], 
   Subscript[S, NO2], 
        Subscript[S, NO], Subscript[S, N2O], Subscript[S, N2], 
   Subscript[S, S], 
        Subscript[S, Mox], Subscript[S, Mred], X}, {t, 0., 
   Subscript[t, 0]}, params, 
      Method -> {"TimeIntegration" -> "StateSpace"}]

transformedData = {ConstantArray[Range@Length[ordinates], 
      Length[abscissae]] // Transpose, 
    ConstantArray[abscissae, Length[ordinates]], ordinates}~
   Flatten~{{2, 3}, {1}};

model[r_, kNH3_, kM1_][i_, t_] := 
  Through[sol[r, kNH3, kM1][t], List][[i]] /; 
   And @@ NumericQ /@ {r, kNH3, kM1, i, t};

fit = NonlinearModelFit[
   transformedData, {model[r, kNH3, kM1][i, t], {0. < r < 1., 
     0. < kNH3 < 1, 0. < kM1 < 0.1}}, {{r, 0.3}, {kNH3,  
     0.05}, {kM1,  0.05}}, {i, t}, 
   Method -> {"NMinimize", 
     Method -> {"SimulatedAnnealing", "PerturbationScale" -> 3}}, 
   EvaluationMonitor :> 
    Print["r=", r, " kNH3=", kNH3, " kM1=", kM1]];
fit["BestFitParameters"]
fit["ParameterConfidenceIntervalTable"]


(* Bioreactor Simulation *)
ClearAll["Global`*"]

rates = {R1[t] == 
    Subscript[r, COD]*X[t]*
     Subscript[S, S][t]/(Subscript[k, S] + Subscript[S, S][t])*
     Subscript[S, Mox][t]/(Subscript[k, Mox] + Subscript[S, Mox][t]), 
   R2[t] == 
    Subscript[r, NO3]*X[t]*
     Subscript[S, NO3][t]/(Subscript[k, NO3] + Subscript[S, NO3][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred1] + Subscript[S, Mred][t]), 
   R3[t] == 
    Subscript[r, NO2]*X[t]*
     Subscript[S, NO2][t]/(Subscript[k, NO2] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred2] + Subscript[S, Mred][t]), 
   R4[t] == 
    Subscript[r, NO]*X[t]*
     Subscript[S, NO][t]/(Subscript[k, NO] + Subscript[S, NO][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred3] + Subscript[S, Mred][t]), 
   R5[t] == 
    Subscript[r, N2O]*X[t]*
     Subscript[S, N2O][t]/(Subscript[k, N2O] + Subscript[S, NO2][t])*
     Subscript[S, Mred][
       t]/(Subscript[k, Mred4] + Subscript[S, Mred][t]),
   RN1[t] == 
    Subscript[rN, NH3]*X[t]*
     Subscript[S, NH3][t]/(Subscript[kN, NH3] + Subscript[S, NH3][t])*
     Subscript[S, Mox][
       t]/(Subscript[kN, Mox1] + Subscript[S, Mox][t])};

conc = {Subscript[S, NH3]'[t] == -RN1[t], 
   Subscript[S, NO3]'[t] == -R2[t], 
   Subscript[S, NO2]'[t] == R2[t] - R3[t], 
   Subscript[S, NO]'[t] == R3[t] - R4[t] + RN1[t], 
   Subscript[S, N2O]'[t] == 0.5 R4[t] - R5[t], 
   Subscript[S, N2]'[t] == R5[t], Subscript[S, S]'[t] == -R1[t], 
   Subscript[S, Mox]'[t] == -(1 - Subscript[Y, H]) R1[t] + R2[t] + 
     0.5 R3[t] + 0.5 R4[t] + R5[t] - RN1[t], 
   X'[t] == Subscript[Y, H]*R1[t], 
   Subscript[S, Mox][t] + Subscript[S, Mred][t] == Subscript[c, tot]};

Subscript[r, COD] = 0.064; Subscript[r, NO3] = 0.045; 
Subscript[r, NO2] = 0.059; Subscript[r, NO] = 0.56; 
Subscript[r, N2O] = 0.23;

Subscript[k, S] = 0.1; Subscript[k, NO3] = 0.018; 
Subscript[k, NO2] = 0.0041; Subscript[k, NO] = 0.000011; 
Subscript[k, N2O] = 0.0025;

Subscript[k, Mox] = 0.0001;
Subscript[k, Mred1] = 0.0015;
Subscript[k, Mred2] = 0.00058;
Subscript[k, Mred3] = 0.000010;
Subscript[k, Mred4] = 0.00024;
Subscript[Y, H] = 0.6;
Subscript[c, tot] = 0.01;

initc = {Subscript[S, NH3][0] == 10, Subscript[S, NO3][0] == 8, 
   WhenEvent[{t == 1.5, t == 3, t == 4.5, t == 6, t == 7.5, 
     t == 9}, {Subscript[S, NH3][t] -> Subscript[S, NH3][t] + 10, 
     Subscript[S, NO3][t] -> Subscript[S, NO3][t] + 2}], 
   Subscript[S, NO2][0] == 0.001, Subscript[S, NO][0] == 0.001, 
   Subscript[S, N2O][0] == 0.001, Subscript[S, N2][0] == 0.001, 
   Subscript[S, S][0] == 170/14, Subscript[S, Mox][0] == 0.005, 
   X[0] == 70};

Subscript[t, 0] = 9;

params = {Subscript[rN, NH3], Subscript[kN, NH3], Subscript[kN, Mox1]};

sol = ParametricNDSolveValue[
  Join[rates, conc, initc], {Subscript[S, NH3], Subscript[S, NO3], 
   Subscript[S, NO2], 
        Subscript[S, NO], Subscript[S, N2O], Subscript[S, N2], 
   Subscript[S, S], 
        Subscript[S, Mox], Subscript[S, Mred], X}, {t, 0., 
   Subscript[t, 0]}, params, 
      Method -> {"TimeIntegration" -> "StateSpace"}]


Plot[{sol[0.7681335, 0.010363, 0.0341977][[1]][t], 
  sol[0.7681335, 0.010363, 0.0341977][[2]][t], 
  sol[0.7681335, 0.010363, 0.0341977][[6]][t]}, {t, 0, 
  Subscript[t, 0]}, PlotRange -> All, AxesLabel -> {t, "S"}, 
 PlotStyle -> {{Blue}, {Red}, {Black}}, 
 PlotLabel -> "Bioreactor Simulation", 
 PlotLabels -> {"\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NH\), \
\(3\)]]\)", 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NO\), \(3\)]]\)", 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(N\), \(2\)]]\)"}]
Plot[{sol[0.7681335, 0.010363, 0.0341977][[1]][t], 
  sol[0.7681335, 0.010363, 0.0341977][[2]][t], 
  sol[0.7681335, 0.010363, 0.0341977][[6]][t]}, {t, 0, 1}, 
 PlotRange -> All, AxesLabel -> {t, "S"}, 
 PlotStyle -> {{Blue}, {Red}, {Black}}, 
 PlotLabel -> "Bioreactor Simulation", 
 PlotLabels -> {"\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NH\), \
\(3\)]]\)", 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NO\), \(3\)]]\)", 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(N\), \(2\)]]\)"}]
g1 = Plot[{sol[0.7681335, 0.010363, 0.0341977][[6]][t]}, {t, 0, 
    Subscript[t, 0]}, PlotRange -> All, AxesLabel -> {t, "S"}, 
   PlotStyle -> {Magenta, Thick}, 
   PlotLabel -> "Bioreactor Simulation", 
   PlotLabels -> {"\!\(\*SubscriptBox[\(S\), SubscriptBox[\(N\), \
\(2\)]]\)"}];

initc2 = {Subscript[S, NH3][0] == 10, Subscript[S, NO3][0] == 8, 
   WhenEvent[{t == 0.80, t == 2*0.80, t == 3*0.80, t == 4*0.80, 
     t == 5*0.80, 
     t == 6*0.80}, {Subscript[S, NH3][t] -> Subscript[S, NH3][t] + 10,
      Subscript[S, NO3][t] -> Subscript[S, NO3][t] + 2}], 
   Subscript[S, NO2][0] == 0.001, Subscript[S, NO][0] == 0.001, 
   Subscript[S, N2O][0] == 0.001, Subscript[S, N2][0] == 0.001, 
   Subscript[S, S][0] == 170/14, Subscript[S, Mox][0] == 0.005, 
   X[0] == 70};

sol = ParametricNDSolveValue[
  Join[rates, conc, initc2], {Subscript[S, NH3], Subscript[S, NO3], 
   Subscript[S, NO2], 
        Subscript[S, NO], Subscript[S, N2O], Subscript[S, N2], 
   Subscript[S, S], 
        Subscript[S, Mox], Subscript[S, Mred], X}, {t, 0., 
   Subscript[t, 0]}, params, 
      Method -> {"TimeIntegration" -> "StateSpace"}]

Plot[{sol[0.7681335, 0.010363, 0.0341977][[1]][t], 
  sol[0.7681335, 0.010363, 0.0341977][[2]][t], 
  sol[0.7681335, 0.010363, 0.0341977][[6]][t]}, {t, 0, 
  Subscript[t, 0]}, PlotRange -> All, AxesLabel -> {t, "S"}, 
 PlotStyle -> {{Blue}, {Red}, {Black}}, 
 PlotLabel -> "Bioreactor Simulation", 
 PlotLabels -> {"\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NH\), \
\(3\)]]\)", 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(NO\), \(3\)]]\)", 
   "\!\(\*SubscriptBox[\(S\), SubscriptBox[\(N\), \(2\)]]\)"}]
g2 = Plot[{sol[0.7681335, 0.010363, 0.0341977][[6]][t]}, {t, 0, 
    Subscript[t, 0]}, PlotRange -> All, AxesLabel -> {t, "S"}, 
   PlotStyle -> {Purple, Thick}, PlotLabel -> "Bioreactor Simulation",
    PlotLabels -> {"\!\(\*SubscriptBox[\(S\), SubscriptBox[\(N\), \(2\
\)]]\)optimized"}];
Show[g1, g2]