I have a group of equations that I was solving pretty consistently with NDsolve. I added a couple more to add complexity to this model, and it all went haywire. I've stared at the code for a couple hours now and I have no idea where it went wrong. any ideas?

nut = 46; khi = 410; ko2 = 52; ro2 = 100; kmrs = 900; kisu = 520; kmp \
= 0.18; kros = .014; kres = 150; kccc1 = 1500; k23 = 7.3; kvp = \
0.020; kc = 27; nc = 5; kfc = 190; nfc = 9; km = 80; nm = 6; km2 = \
200; nm2 = 4; kv1 = 1.2; nv1 = 3; kvs = 190; nvs = 8; k32 = 2; n32 = \
10; k3s = 210; n3s = 10; volc = .65; volm = .1; volv = .25; vhi = \
khi; kmhi = 14; kmmrs = 12; vmrs = kmrs; vkccc1 = kccc1; kmccc1 = \
5.5; kmisu = 610; kmres = 9;
vrim = 90;
kmrim = 6;
kmd = 52;
kmmd = 300;
a = .2; odec = -(a c[t]) - (
  vkccc1 volv c[
    t] (1 - 1/((c[t]/kv1)^nv1 + 1)) (1 - 1/((fs[t]/kvs)^nvs + 1)))/(
  volc (c[t] + kmccc1)) + (
  nut vhi)/((kmhi + nut) ((c[t]/kc)^nc + 1) ((fs[t]/kfc)^nfc + 1)) - (
  vmrs volm c[t])/(
  volc (c[t] + kmmrs) ((c[t]/km)^nm + 1) ((fs[t]/km2)^nm2 + 1)) - (
  vrim volm c[t])/(
  volc (c[t] + kmrim) ((c[t]/km3)^nm3 + 1) ((fs[t]/km4)^nm4 + 
     1)); oderos = 
 kmp fm[t] o2[t] - a ros[t]; odefm = -(a fm[t]) + (
  vmrs c[t])/((c[t] + kmmrs) ((c[t]/km)^nm + 1) ((fs[t]/km2)^nm2 + 
     1)) - kisu fm[t] fe[t] - kmp fm[t] o2[t];

odefmb = (vrim volm c[t])/(
   volc (c[t] + kmrim) ((c[t]/km3)^nm3 + 1) ((fs[t]/km4)^nm4 + 1)) - (
   kmd fmb[t])/(kmmd + fmb[t]) - a fmb[t];
odefe = (kmd fmb[t])/(kmmd + fmb[t]) - kisu fm[t] fe[t] - a fe[t];
odefs = kisu fm[t] fe[t] - a fs[t];
odemp = kmp fm[t] o2[t] - a mp[t]; odeo2 = -(kmp fm[t] o2[t]) - (
  kres fs[t] o2[t])/(kmres + o2[t]) + 
  ko2 (ro2 - o2[t]); odef2 = -(a f2[t]) - 
  k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) (1 - 
     1/((fs[t]/k3s)^n3s + 1)) + (
  vkccc1 c[t] (1 - 1/((c[t]/kv1)^nv1 + 1)) (1 - 
     1/((fs[t]/kvs)^nvs + 1)))/(c[t] + kmccc1); odef3 = -(a f3[t]) + 
  k23 f2[t] (1 - 1/((c[t]/k32)^n32 + 1)) (1 - 
     1/((fs[t]/k3s)^n3s + 1)) - kvp f3[t]; odevp = 
 kvp f3[t] - a vp[t]; vars = {c[t], f2[t], f3[t], vp[t], fm[t], 
  fmb[t], fe[t], fs[t], mp[t], o2[t], ros[t]}; solution = 
 NDSolve[{Derivative[1][c][t] == odec, Derivative[1][f2][t] == odef2, 
   Derivative[1][f3][t] == odef3, Derivative[1][fm][t] == odefm, 
   fmb'[t] == odefmb, fe'[t] == odefe, Derivative[1][fs][t] == odefs, 
   Derivative[1][mp][t] == odemp, Derivative[1][o2][t] == odeo2, 
   Derivative[1][vp][t] == odevp, Derivative[1][ros][t] == oderos, 
   c[0] == 31.3972, f2[0] == 14.1354, f3[0] == 1.53991, 
   fm[0] == 624.774, fmb[0] == 1, fe[0] == 100, fs[0] == 81.3915, 
   mp[0] == 7673.09, o2[0] == 5.02413, vp[0] == 8.44881, 
   ros[0] == 7673.09}, vars, {t, 0, 31}]; vars; vars1 = 
 vars /. solution /. t -> 31