When using NDsolve, how to determine the positions of steady states?

[Update notice: NDSolve code below used to work without the Method option, but somewhere around V10.4, the option became necessary. Also see this.]

Based on the OP's original "set-up":

des = {x[1]'[t], x[2]'[t], x[3]'[t], x[4]'[t], x[5]'[t],
    x[6]'[t], x[7]'[t], x[8]'[t], 
    x[9]'[t]} == {-k[1] x[1][t] x[3][t] + k[2] x[5][t] + 
     k[3] x[5][t] - k[7] x[1][t] x[7][t] + k[8](*[t]*) x[8][t],
    -k[4] x[2][t] x[4][t] + k[5] x[6][t] + k[6] x[6][t] - 
     k[9] x[2][t] x[7][t] + k[10] x[9][t],
    -k[1] x[1][t] x[3][t] + k[2](*[t]*) x[5][t] + 
     k[6] x[6][t], -k[4] x[2][t] x[4][t] + k[3] x[5][t] + k[5] x[6][t],
    k[1] x[1][t] x[3][t] - k[2] x[5][t] - k[3] x[5][t],
    k[4] x[2][t] x[4][t] - k[5] x[6][t] - k[6] x[6][t],
    -k[7] x[1][t] x[7][t] - k[9] x[2][t] x[7][t] + k[8] x[8][t] + 
     k[10] x[9][t], k[7] x[1][t] x[7][t] - k[8] x[8][t],
    k[9] x[2][t] x[7][t] - k[10] x[9][t]};

init = {T[1], T[2], T[3], 0, 0, 0, T[4], 0, 0};

Use the Norm of the derivative (with WhenEvent) to determine when the system is close to a steady-state (same idea as belisarius, who answered first).

vars = Array[x, 9];
dvars = Thread[Derivative[1][vars]];

Block[{k, T, ssthreshold},
  k[n_] := k[n] = (SeedRandom[n]; RandomReal[]);
  T[n_] := T[n] = (SeedRandom[n + 10]; RandomReal[]);
  ssthreshold = 1.*^-4;
  (* Print[des]; *) (* to see the ODE *)
  {sol} = 
   NDSolve[{des, Through[vars[0]] == init, 
     With[{df = Through[dvars[t]]}, 
      WhenEvent[Norm[df] < ssthreshold, x[7][t] -> x[7][t] + 0.1]]}, 
    vars, {t, 0, 200}, MaxSteps -> 100000,
    Method -> {"EquationSimplification" -> "Residual"}  (* needed as of V10.4 or so *)
    ]
  ];

Plot @@ {Through[vars[t]] /. sol, Flatten@{t, x[1]["Domain"] /. sol}, 
  PlotLegends -> Automatic}

Because the equations are nonlinear, there is no assurance that equilibrium solutions even exist for a given set of k and initial conditions. However, if they do exist, an alternative, and perhaps more informative, approach is to compute them directly:

xs = Solve[Thread[Table[0, {i, Length[des[[1]]]}] == des[[2]] /. x[n_][t] -> x[n]],
     Table[x[i], {i, Length[des[[1]]]}]]

Although Solve warns

Solve::svars: Equations may not give solutions for all "solve" variables. >>

playing a bit with Reduce strongly suggests that Solve gives all the equilibrium solutions here:

{* {{x[3] -> ((k[2] + k[3]) x[5])/(k[1] x[1]), 
      x[4] -> (k[3] (k[5] + k[6]) x[5])/(k[4] k[6] x[2]), 
      x[6] -> (k[3] x[5])/k[6], x[8] -> (k[7] x[1] x[7])/k[8], 
      x[9] -> (k[9] x[2] x[7])/k[10]},
   {x[1] -> 0, x[4] -> 0, x[5] -> 0, x[6] -> 0, x[8] -> 0, x[9] -> (k[9] x[2] x[7])/k[10]},
   {x[1] -> 0, x[2] -> 0, x[5] -> 0, x[6] -> 0, x[8] -> 0, x[9] -> 0},
   {x[2] -> 0, x[3] -> 0, x[5] -> 0, x[6] -> 0, x[8] -> (k[7] x[1] x[7])/k[8], 
      x[9] -> 0}} *}

(For each of the four solutions, any x that do not appear can take any value.) With the equilibrium solutions known, one can linearize the equations des about them, Laplace transform the resulting linear equations, and Solve them to obtain the Eigensystems, completely determining the perturbed solutions in closed form.

Addendum

In answer to the question posed below in a comment, the right side of des can be linearized about an equilibrium solution without difficulty. (The left side already is linear.)

desl = Coefficient[des[[2]] /. x[n_][t] -> x[n] + e dx[n], e]

(* {dx[5] k[2] + dx[5] k[3] + dx[8] k[8] - dx[3] k[1] x[1] - 
      dx[7] k[7] x[1] - dx[1] k[1] x[3] - dx[1] k[7] x[7], 
    dx[6] k[5] + dx[6] k[6] + dx[9] k[10] - dx[4] k[4] x[2] - 
      dx[7] k[9] x[2] - dx[2] k[4] x[4] - dx[2] k[9] x[7], 
    dx[5] k[2] + dx[6] k[6] - dx[3] k[1] x[1] - dx[1] k[1] x[3], 
    dx[5] k[3] + dx[6] k[5] - dx[4] k[4] x[2] - dx[2] k[4] x[4], 
    -dx[5] k[2] - dx[5] k[3] + dx[3] k[1] x[1] + dx[1] k[1] x[3], 
    -dx[6] k[5] - dx[6] k[6] + dx[4] k[4] x[2] + dx[2] k[4] x[4], 
    dx[8] k[8] + dx[9] k[10] - dx[7] k[7] x[1] - dx[7] k[9] x[2] - dx[1] k[7] x[7] - 
      dx[2] k[9] x[7], 
    -dx[8] k[8] + dx[7] k[7] x[1] + dx[1] k[7] x[7], 
    -dx[9] k[10] + dx[7] k[9] x[2] + dx[2] k[9] x[7]} *)

At this point, equilibrium values of x are substituted (via Rules) into the last result and the Eigenvaluess computed. These Eigenvaluess are the growth rates of the perturbations dx. (Complex Eigenvaluess indicate oscillatory solutions.) For instance, if all equilibrium values are 0, then

Eigenvalues[(CoefficientArrays[desl, Array[dx, 9]][[2]] // Normal) /. 
  Thread[Array[x, 9] -> 0]]
(* {0, 0, 0, 0, 0, -k[2] - k[3], -k[5] - k[6], -k[8], -k[10]} *)

Hence, if any of k[2] + k[3], k[5] + k[6], k[8], or k[10] are negative, then perturbations grow exponentially. More generally, one would substitute one of the four equilibrium solutions xs derived earlier into the perturbed equations.

Another advantage of the procedure described here is this: In the event that an equilibrium admits exponentially growing perturbations, then attempts to solve for that equilibrium using NDSolve are likely to be unstable too.

An example perturbing a simple RL circuit when it's about to reach the steady state.

R = 1;
L = 1;
V = 1;
SeedRandom[43];
sol = NDSolve[{R i[t] + L i'[t] == V, i[0] == 0,
              WhenEvent[{i'[t] ==  .01}, i[t] -> i[t] + RandomReal[{-.5, .5}]],
              WhenEvent[{i'[t] == -.01}, i[t] -> i[t] + RandomReal[{-.5, .5}]]}, 
              i, {t, 0, 10}]
Plot[(i /. sol[[1]])[t], {t, 0, 10}, PlotRange -> All]