Parametric plot refuses to assign unique plot styles

Question

I want to plot the motion on n bodies. To do this I have set up three variables that return a list of variables to be evaluated by an InterpolationFunction through s = NDSolve[...].

That is

xm = {x1[t], x2[t], ..., xn[t]}
ym = {y1[t], y2[t], ..., yn[t]}
zm = {z1[t], z2[t], ..., zn[t]}

I then visualize it with:

ParametricPlot3D[Flatten[Thread[Transpose[{xm, ym, zm}] /. s], 1], {t, t0, tmax}]

Where the functions return a n x 3 list of interpolating functions

Unfortunately, when I try to use PlotStyle to differentiate between the plots, all n plots take on the same style.

The entire current code is a follows:

n = 3;
m = {1, 0.005, 0.01};
positions = {{0, 0, 0}, {4, 0, 0.2}, {-2, 0, 0}};
velocities = {{0, 0.5 m[[2]] - 0.6 m[[3]], 
    0.025 m[[2]]}, {0, -0.5, -0.025}, {0, 0.6, 0}};
time = 100;
(*Here Endeth user input*)

colours = 
  Table[RGBColor[RandomReal[], RandomReal[], RandomReal[]], {i, n}];
positions = Transpose[positions];
velocities = Transpose[velocities];
xt = ToExpression[Table["x" <> ToString[i], {i, n}]];
xm = Through[xt[t]];
xz = Through[xt[0]];
yt = ToExpression[Table["y" <> ToString[i], {i, n}]];
ym = Through[yt[t]];
yz = Through[yt[0]];
zt = ToExpression[Table["z" <> ToString[i], {i, n}]];
zm = Through[zt[t]];
zz = Through[zt[0]];
rm = Flatten[
    Table[If[i != j, 
      Sqrt[(xm[[j]] - xm[[i]])^2 + (ym[[j]] - ym[[i]])^2 + (zm[[j]] - 
           zm[[i]])^2]], {i, n}, {j, n}]] /. Null -> Sequence[];
xf1 = (Flatten[
      Table[If[i != j, (xm[[j]] - xm[[i]]) m[[j]]], {i, n}, {j, 
        n}]] /. Null -> Sequence[])/rm^3;
xf = Thread[
   D[D[xm, t], t] == 
    Table[Total[Take[xf1, {(n - 1) i - n + 2, (n - 1) i}]], {i, n}]];
yf1 = (Flatten[
      Table[If[i != j, (ym[[j]] - ym[[i]]) m[[j]]], {i, n}, {j, 
        n}]] /. Null -> Sequence[])/rm^3;
yf = Thread[
   D[D[ym, t], t] == 
    Table[Total[Take[yf1, {(n - 1) i - n + 2, (n - 1) i}]], {i, n}]];
zf1 = (Flatten[
      Table[If[i != j, (zm[[j]] - zm[[i]]) m[[j]]], {i, n}, {j, 
        n}]] /. Null -> Sequence[])/rm^3;
zf = Thread[
   D[D[zm, t], t] == 
    Table[Total[Take[zf1, {(n - 1) i - n + 2, (n - 1) i}]], {i, n}]];
pos = {Thread[xz == positions[[1]]], Thread[yz == positions[[2]]], 
   Thread[zz == positions[[3]]]};
vel = {Thread[D[xm, t] == velocities[[1]]] /. t -> 0, 
   Thread[D[ym, t] == velocities[[2]]] /. t -> 0, 
   Thread[D[zm, t] == velocities[[3]]] /. t -> 0};
s = NDSolve[Flatten[{xf, yf, zf,
     vel, pos}], Flatten[{xt, yt, zt}], {t, 0, time}, 
   AccuracyGoal -> 12, PrecisionGoal -> 12, MaxSteps -> \[Infinity]];

Manipulate[
 Show[ParametricPlot3D[
   Flatten[Thread[Transpose[{xm, ym, zm}] /. s], 1], {t, t0, a}, 
   PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}}, ImageSize -> Large, 
   PlotStyle -> colours], 
  ParametricPlot3D[
   Flatten[Thread[Transpose[{xm, ym, zm}] /. s], 1], {t, a - 0.05, a},
    PlotStyle -> 
    Transpose[{ConstantArray[Thickness[0.02], n], colours}]]], {t0, 0,
   time - 1, 1}, {{a, 1, "t"}, t0, time, 1}]