I can only propose a crude starting point based on a more modest work I did. My goal was to determine the bleaching coefficient of a fluorophore when submitted to varying illumination intensity over time. I had the theoretical law for the bleaching over time as a function of intensity (depending on unknown parameters), the amount of fluorophore over time, and the starting intensity and intensity at the end of the experiment. The goal was to retrieve the parameters for the bleaching law and the whole intensity over time.
Anyway. What I did was to construct (using MATLAB) a function to minimize. In your case, given that you know $x(t)$ experimentally, I suggest to build a function $O(x, y)$ such that:
$$ O(x, y) = \sum_t |x(t) - y(t)|^2 $$
with $x(t)$ being your experimental data and $y(t)$ being the solution of your differential system for some parameters $k, M, a$. Put this function as the objective in a simplex optimizer that will find its minimum, allowing $k, M, a$ to vary.
This is done through numerical solution. So you have a function to minimize, say $F(k,a,M)$ that returns the value of $O(x,y)$ (which is scalar btw) for the current numerical resolution with the given $k,a,M$.
At each step, the whole system will be solved again, and the solution will be matched against the experimental data until they fit.
I have found this to be working quite well in my case, as long as you are sure you have the correct model to describe the experimental data. Otherwise, expect some seriously irrelevant results.