{ diff(x1, t) = (1*k22*k1*x1*(1 - k2*x1) + (-1)*k22*(k3*x2 + k14*x3 + k15*x4)*x1)/k22,
  diff(x2, t) = (1*k22*(k16 + k20*x2*x1*x1/(k21 + x1*x1)) + (-1)*k22*((k4*x1 + k5*x3)*x2 + k8*x2))/k22,
  diff(x3, t) = (1*k22*k17 + (-1)*k22*((k9*x4 + k6*x2 - k7*x1)*x3 - k11*x3))/k22,
  diff(x4, t) = (1*k22*(k10*x3*x1 + k19*x2*x1) + (-1)*k22*(k12*x4*x1 + k18*x2*x4*x4 + k13*x4))/k22 }