Clear["Global'*"];
 errortrial1 = 
  Table[ro2 = 3400; kmrs = 70; kisu = 1.4; kmp = .5; kres = 25;
   kccc1 = 38; k23 = 16; kvp = .04; kc = 21; nc = 5; km1 = 24; nm1 = 3;
   km2 = 15; nm2 = 4.5; kv1 = 17; nv1 = 6; kv2 = 200; nv2 = 4; 
   k32 = 18; n32 = 9; a = .5;
   odec = -(a c[t]) - (
     kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)))/((f3[t]/kv2)^nv2 + 1) - (
     kmrs c[t])/((fm[t]/km2)^nm2 + 1) + (khi nut)/((c[t]/kc)^nc + 1);
  odefm = -(a fm[t]) + (kmrs c[t])/((fm[t]/km2)^nm2 + 1) - 
    kisu fm[t] - kmp fm[t] o2[t];
  odefs = kisu fm[t] - a fs[t];
  odemp = kmp fm[t] o2[t] - a mp[t];
  odeo2 = -(kmp fm[t] o2[t]) - kres fs[t] o2[t] + 
    ro2/((o2[t]/245)^9 + 1);
  odef2 = -(a f2[t]) - k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) + (
    kccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)))/((f3[t]/kv2)^nv2 + 1);
  odef3 = -(a f3[t]) + k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) - 
    kvp f3[t];
  odevp = kvp f3[t] - a vp[t];
  vars = {khi, nut, c[t], f2[t], f3[t], fm[t], fs[t], mp[t], o2[t], 
    vp[t], rmsd};
  solution = 
   NDSolve[{Derivative[1][c][t] == odec, 
     Derivative[1][f2][t] == odef2, Derivative[1][f3][t] == odef3, 
     Derivative[1][fm][t] == odefm, Derivative[1][fs][t] == odefs, 
     Derivative[1][mp][t] == odemp, Derivative[1][o2][t] == odeo2, 
     Derivative[1][vp][t] == odevp, c[0] == 30, f2[0] == 100, 
     f3[0] == 100, fm[0] == 50, fs[0] == 50, mp[0] == 0, o2[0] == 100,
      vp[0] == 100}, vars, {t, 0, 31}];
  vars1 = vars /. solution /. t -> 31;
  vars2 = Flatten[vars1];
  If[nut == 4, {fecell = 153, fef3 = 0, fem = 97, fevp = 0}, 
   "nut not 4"];
  If[nut == 7, {fecell = 250, fef3 = 100, fem = 60, fevp = 30}, 
   "nut not 7"];
  If[nut == 16, {fecell = 395, fef3 = 250, fem = 80, fevp = 0}, 
   "nut not 16"];
  If[nut == 46, {fecell = 440, fef3 = 295, fem = 65, fevp = 45}, 
   "nut not 46"];
  If[nut == 106, {fecell = 440, fef3 = 350, fem = 40, fevp = 20}, 
   "nut not 106"];
  If[nut == 1006, {fecell = 455, fef3 = 360, fem = 30, fevp = 30}, 
   "nut not 1006"];
  If[nut == 10006, {fecell = 465, fef3 = 365, fem = 25, fevp = 35}, 
   "nut not 10006"];
  rmsd = ((((Part[vars2, 3] + Part[vars2, 4] + Part[vars2, 5] + 
           Part[vars2, 6] + Part[vars2, 7] + Part[vars2, 8] + 
           Part[vars2, 10] - fecell)/
         fecell)^2 + ((Part[vars2, 5] - fef3)/
         fecell)^2 + ((Part[vars2, 6] + Part[vars2, 7] + 
           Part[vars2, 8] - fem)/
         fecell)^2 + ((Part[vars2, 10] - fevp)/fecell)^2));
  vars2, {khi, 10, 16, 1}, {nut, {4, 7, 16, 46, 106, 1006, 10006}}]
Map[Mean[{rmsd, #}] &, errortrial1, {3}]
For the first part, do you mean something like this? I get values, but I am not entirely sure what they are.