How does one use a fitting routine such as NonlinearModelFit with a ParametricFunction? Example below:

First let's generate some data:

 sample[t_] = (0.002 + 101 t - 461000 t^2 + 2.218 10^9 t^3 - 
 3.64 10^12 t^4 + 3.17 10^15 t^5 ) Exp[-8653 t];
 data = Table[{t, sample[t] + 
 RandomVariate[NormalDistribution[0, 0.00001]]}, 
 {t, 0, 0.002, 0.000004}];
 ListPlot[data]

[![enter image description here][1]][1]

Now lets define a model:

 rateeqs = {a'[t] == -k1a a[t] - k12 a[t] c70gs + k21 b[t] c60gs, 
 b'[t] == -k1b b[t] - k21 b[t] c60gs + k12 a[t] c70gs, 
 a[0] == a0, b[0] == b0};
 c60gs = c70gs = 5;
 maxTime = 0.0025;
 e60 = 19060;
 e70 = 948;
 fitFunc[t_] = e60 a[t] + e70 b[t];
 params = {k1a, k1b, k12, k21, a0, b0};
 initGuesses = {8000, 100, 4500, 2000, 5. 10^-8, 8 10^-7};

Now we generate a parametric solution to the model:

 solution[t_] = ParametricNDSolveValue[rateeqs, fitFunc[t], {t, 0, maxTime}, 
 params, WorkingPrecision->30]

And we verify that the solution works:

 Show[ListPlot[data, PlotStyle -> Red, PlotRange -> Full], 
 Plot[(solution[t] @@ initGuesses), {t, 0, maxTime}]]

[![enter image description here][2]][2]

It would seem that I should be able to fit this model to the sample data starting from the initGuesses values for the parameters using something like this:

 modelFit = NonlinearModelFit[data, solution[t] @@ params, 
 {params, initGuesses}\[Transpose], t];

But I get errors:

[![enter image description here][3]][3]

I'm sure it's a syntax issue, but I have tried so many different possibilities that the screen is swimming in front of my eyes. Any help is appreciated!

 [1]: https://i.sstatic.net/itIOP.png
 [2]: https://i.sstatic.net/YRGTZ.png
 [3]: https://i.sstatic.net/T4oPK.png