How to fit the parameters of differential equations with known data?

This kind of problem is known in the literature as a nonlinear state-space system identification. Several algorithms have been proposed in the literature to solve these problems. I think a good starting point would be (1) and the references therein, in particular the works of L. Ljung.

As far I know, in general if you don't have a good initial estimate of the parameters $k_i$, there is no guarantee that the algorithm will converge to a global solution.

(1) Schön, Thomas B., Adrian Wills, and Brett Ninness. "System identification of nonlinear state-space models." Automatica 47.1 (2011): 39-49.

P.S. Many of these identification algorithms have been implemented in MATLAB (maybe there is an equivalent implementation in Octave). See for instance this freely available third party toolbox and this example from the system identification toolbox of MATLAB.

Update (08/23/2014):

Following this example and using the first data set I found the following solution: $$\begin{array}{cccc} k_1 &= &1.100350 &\pm& 8.3610 \\ k_2 &= &2.588940 &\pm& 9.6769 \\ k_3 &= &0.400037 &\pm& 14.8607 \\ k_4 &= &0.338143 &\pm& 7.8111 \\ k_5 &= &3.757100 &\pm& 168.1370 \\ k_6 &= &0.915842 &\pm& 11.7712 \end{array}$$

The algorithm used was trust-region reflective Newton method of nonlinear least-squares, where the cost is the sum of squares of errors between the measured and simulated outputs.

Notice that the high standard deviation must be due to the small amount of samples.

If you have access to Maple and the Global Optimization Toolbox, you could try something like this.

data:=[[0.0000,9.1300,0.0931,0.0899,0.1000,0.0000],
[30.0000,8.9300,0.1270,0.1230,0.2270,0.0049],
[60.0000,8.6000,0.1510,0.1390,0.4920,0.0153],
[90.0000,8.2800,0.1540,0.1490,0.7780,0.0249],
[120.0000,7.9700,0.1540,0.1570,1.0700,0.0329],
[150.0000,7.8600,0.1540,0.1600,1.1700,0.0348],
[180.0000,7.8100,0.1530,0.1530,1.2100,0.0404],
[210.0000,7.7700,0.1400,0.1420,1.2800,0.0432]]:
des:=[diff(y1(t),t) = -k1*y1(t) - k2*y1(t),
diff(y2(t),t) = k2*y1(t) - k3*y2(t),
diff(y3(t),t) = k1*y1(t) +k3*y2(t) - k4*y3(t),
diff(y4(t),t) = k4*y3(t) - k5*y2(t)*y4(t) + k6*y5(t),
diff(y5(t),t) = k5*y2(t)*y4(t) - k6*y5(t)]:
ics := seq((y || i)(0) = data[1, i+1], i = 1 .. 5):
res := dsolve({ics, des[]}, numeric, parameters = [k1, k2, k3, k4, k5, k6]):
timeList := [0, 30, 60, 90, 120, 150, 180, 210];
sse := proc (k1, k2, k3, k4, k5, k6) 

   res(parameters = [k1, k2, k3, k4, k5, k6]); 

   add(
   (rhs(select(has, res(timeList[i]), y1)[])-data[i, 2])^2+
   (rhs(select(has, res(timeList[i]), y2)[])-data[i, 3])^2+
   (rhs(select(has, res(timeList[i]), y3)[])-data[i, 4])^2+
   (rhs(select(has, res(timeList[i]), y4)[])-data[i, 5])^2+
   (rhs(select(has, res(timeList[i]), y5)[])-data[i, 6])^2, 
   i = 2 .. 8);

end proc;
c := GlobalOptimization:-GlobalSolve('sse'(k1, k2, k3, k4, k5, k6), k1 = 0 .. 1, k2 = 0 .. 1, k3 = 0 .. 1, k4 = 0 .. 1, k5 = 0 .. 1, k6 = 0 .. 1, timelimit = 10);
   [0.219132447080011505, [k1 = 8.52482740113834e-4, 

     k2 = 5.2683998680924474e-5, k3 = 0.0, 

     k4 = 0.05113239298267808, 

     k5 = 0.004363021255887466, k6 = 0.0]]
res(parameters = c[2]):
p:=Array(1..5):
for n from 1 to 5 do
   p[n]:=plots:-display(plots:-odeplot(res, [t, (y || n)(t)], t = 0 .. 210),plots:-pointplot([seq([data[i, 1], data[i, n+1]], i = 1 .. 8)]));
end do;
plots:-display(p)

The code is a parameterized numeric solution of the differential equations (I shifted all the data down by 30s, so I could get the ICs at t=0) followed by a global optimization of the least-squares problem.

I haven't tested the code for correctness, and if you let the optimizer run for longer or change the optimization bounds, you may get better results.