Estimating parameters for pair of differential equations efficiently

I'm trying to estimate the parameters for the following pair of differential equations

meq = gamma*v[t]*m[t] - mu*m[t];
veq = v[t]*(k - epsilon - mu) - v[t]^2*k + m[t]*l - m[t]^2*l - 
   m[t]*v[t]*(k + l + gamma);
equations = {m'[t] == meq, v'[t] == veq};

where gamma, k, l and epsilon are parameters to be estimated, and mu is a known parameter with value 4.58*10^(-5). Also I have initial conditions given by

initials = {m[0] == 2.03, v[0] == 3.09};
mu = 4.58*10^(-5)  

and the data for the function m is the following

data = {{60, 2.13597}, {300, 2.27247}, {390, 2.29472}, {420, 2.40096}, {630, 
  2.59312}, {660, 2.84918}, {780, 2.93677}, {1020, 3.02945}, {1110, 
  3.04794}, {1140, 3.05796}, {1140, 3.08739}, {1380, 3.21218}, {1380, 
  3.2873}, {1500, 3.3347}, {1680, 3.44467}, {1710, 3.47574}, {1710, 
  3.48421}, {1830, 3.50433}}

I defined a parametric solution for the pair of equations with

pfun = ParametricNDSolveValue[Join[equations, initials], m, {t, 0, 2000}, 

{epsilon, k, l, gamma}]

And then I tried to use FindFit as

fit = FindFit[data, pfun[epsilon, k, l, gamma][t], {epsilon, k, l, gamma}, 
t]

However this method seems to be super inefficient as it seems to take unreasonable amount of time to complete the fit. Is there a way maybe to form the parametric function in a way that the FindFit works faster or use the FindFit more efficiently? Of course if you have another faster method in mind it's very welcome.

I also tried to solve this with some different programs than Mathematica and they seems to be working much faster but in those I have no idea what algorithms they use and also in them I have serious data issues.

Answer

Here you need to use a different model. In this model, too, not everything is smooth, but the result is clear and fast

meq = gamma*v[t]*m[t] - mu*m[t];
veq = v[t]*(k - epsilon - mu) - v[t]^2*k + m[t]*l - m[t]^2*l - 
   m[t]*v[t]*(k + l + gamma);
equations = {m'[t] == meq, v'[t] == veq};
initials = {m[0] == 2.03, v[0] == 3.09};
mu = 4.58*10^(-5);
data = {{60, 2.13597}, {300, 2.27247}, {390, 2.29472}, {420, 
    2.40096}, {630, 2.59312}, {660, 2.84918}, {780, 2.93677}, {1020, 
    3.02945}, {1110, 3.04794}, {1140, 3.05796}, {1140, 
    3.08739}, {1380, 3.21218}, {1380, 3.2873}, {1500, 3.3347}, {1680, 

    3.44467}, {1710, 3.47574}, {1710, 3.48421}, {1830, 3.50433}};


 model[epsilon_?NumberQ, k_?NumberQ, l_?NumberQ, 
   gamma_?NumberQ] := (model[epsilon, k, l, gamma] = 
    First[m /. 
      NDSolve[{v'[t] == 
         v[t]*(k - epsilon - mu) - v[t]^2*k + m[t]*l - m[t]^2*l - 
          m[t]*v[t]*(k + l + gamma), 
        m'[t] == gamma*v[t]*m[t] - mu*m[t], m[0] == 2.03, 

        v[0] == 3.09}, {m, v}, {t, 2000}]]);

 fit = 
 FindFit[data, 
  model[epsilon, k, l, gamma][
   t], {{epsilon, -4.8}, {k, 4.9}, {l, 0.}, {gamma, .1}}, t, 
  PrecisionGoal -> 4, AccuracyGoal -> 4]

(*Out[]= {epsilon -> -4.78195, k -> 4.89863, l -> -0.00336569, 
 gamma -> 0.0769293}*)


 Show[
 Plot[model[epsilon, k, l, gamma][t] /. fit, {t, 0, 2000}], 
 ListPlot[data, PlotStyle -> Orange]]