How to fit data with numerical solution of system of parametric ODE?

Question

I need to find the parameters (k1,k2,k3,k4,k1r,k2r,k3r,k4r) that fit my data (list of [Intensity, time]) using the function b(t), solution of the following system of parametric ODE:

{a'[t] == -k1*a[t] - k2*a[t] + k1r*g[t] + k2r*p[t],
 g'[t] == -k1r*g[t] - k3*g[t] + k1*a[t] + k3r*b[t],
 p'[t] == -k4*p[t] - k2r*p[t] + k4r*c[t] + k2*a[t],
 b'[t] == -k3r*b[t] + k3*g[t],
 c'[t] == k4*p[t] - k4r*c[t]}

Initial conditions:

a[0] == 1,
c[0] == 0,
b[0] == 0,
p[0] == 0,
g[0] == 0,

The system of ODE has no analytical solution, so I tried to use NDSolve, but it is of course not working because the ODEs are parametric. I guess ParametricNDSolve would help me, but unfortunately it is not implemented in the version of Mathematica I am using (Mathematica 8). Do you have any idea about how to address this problem, avoiding the ParametricNDSolve function?