Fitting NDSolve solution to data

Question

I have a system of DDEs with 2 free parameters I solved with NDSolve, as you can see below

    n = 1;
system[l_?NumericQ, 
  eps_] := {m'[t] == 2/n*((1 - m[t])/2*lp[t] - (1 + m[t])/2*lm[t]),
  lp'[t] == 
   1/n *((1 + m[t])/2*((1 - m[t])/2 - lm[t])) - 
    1/n *((1 + m[t - l])/2*((1 - m[t - l])/2 - lm[t - l])),
  lm'[t] == 
   1/n *((1 - m[t])/2*((1 + m[t])/2 - lp[t])) - 
    1/n *((1 - m[t - l])/2*((1 + m[t - l])/2 - lp[t - l])),
  m[t /; t <= l] == eps*(2*E^(-t/(2 n)) - 1), 
  lp[t /; t <= l] == (1 - eps)/2*(1 - E^(-t/(2 n))), 
  lm[t /; t <= l] == (1 - eps)/2*(1 - E^(-t/(2 n)))}
sol[l_?NumericQ, eps_] := 
  NDSolve[system[l, eps], {m, lp, lm}, {t, 0, 300}];

This method works and gives me the chance of exploring the solutions of my system varying the two parameters l and eps.

Now I would like to fit a dataset I imported using the solution provided by NDSolve as a model, but I just can't define the model correctly to make it work with FindFit or similars.

I tried with

model[t_, a_, b_, c_] := m[t + a] /. sol[b, c]

and variations (tried Evaluate, adding ?NumericQ and so on) but it always gives me the error

ReplaceAll::reps: {sol[b,c]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

I think the problem is Mathematica attempts at evaluating m[t] with sol, which contains no numbers, just parameters, and it fails at it. How can I fix it to make it suitable for fitting?

And please excuse me if I couldn't find the answer among similar questions posted, but I find Mathematica very criptic (yet useful!).

Thanks in advance

EDIT: Here's the data I want to fit (orange) and the solution of my system I'd like to use (I have manually found it with Manipulate).

Dataset to fit:

t, m
    {{0, 0.007}, {1, -0.011}, {2, -0.017}, {3, -0.036}, {4, -0.044}, {5, \
-0.045}, {6, -0.043}, {7, -0.053}, {8, -0.059}, {9, -0.06}, {10, \
-0.061}, {11, -0.077}, {12, -0.084}, {13, -0.055}, {14, -0.054}, {15, \
-0.057}, {16, -0.061}, {17, -0.058}, {18, -0.041}, {19, -0.024}, {20, \
-0.018}, {21, 0.006}, {22, 0.032}, {23, 0.017}, {24, -0.002}, {25, 
  0.02}, {26, 0.042}, {27, 0.067}, {28, 0.062}, {29, 0.062}, {30, 
  0.053}, {31, 0.052}, {32, 0.06}, {33, 0.068}, {34, 0.088}, {35, 
  0.065}, {36, 0.052}, {37, 0.074}, {38, 0.072}, {39, 0.075}, {40, 
  0.087}, {41, 0.074}, {42, 0.046}, {43, 0.018}, {44, 0.014}, {45, 
  0.015}, {46, 0.005}, {47, 0.007}, {48, 0.019}, {49, 0.014}, {50, 
  0.004}, {51, -0.002}, {52, 
  0.001}, {53, -0.006}, {54, -0.021}, {55, -0.024}, {56, -0.018}, \
{57, -0.026}, {58, -0.03}, {59, -0.047}, {60, -0.047}, {61, -0.063}, \
{62, -0.081}}

UPDATE: I hope this might be useful to somebody. I've found a way of passing my model to FindFit or NonLinearModelFit:

n = 1;
model[l_?NumberQ, 
  eps_?NumberQ] := (model[a, b] = 
   Module[{m, lp, lm, t}, 
    First[m /. 
      NDSolve[{m'[t] == 2/n*((1 - m[t])/2*lp[t] - (1 + m[t])/2*lm[t]),
        lp'[t] == 
         1/n *((1 + m[t])/2*((1 - m[t])/2 - lm[t])) - 
          1/n *((1 + m[t - l])/2*((1 - m[t - l])/2 - lm[t - l])),
        lm'[t] == 
         1/n *((1 - m[t])/2*((1 + m[t])/2 - lp[t])) - 
          1/n *((1 - m[t - l])/2*((1 + m[t - l])/2 - lp[t - l])),
        m[t /; t <= l] == eps*(2*E^(-t/(2 n)) - 1), 
        lp[t /; t <= l] == (1 - eps)/2*(1 - E^(-t/(2 n))), 
        lm[t /; t <= l] == (1 - eps)/2*(1 - E^(-t/(2 n)))}, {m, lp, 
        lm}, {t, 0, 300}]]])
nlm = FindFit[
  death, {model[a, b][
    t + c], {11 < a < 12}, {0.0001 < b < 0.01}, {150 < c < 180}}, {{a,
     11.2}, {b, 0.001}, {c, 160}}, t]

If I run this, my code works, but the answer isn't quite the one I expected (although I set the initial values of the parameters with the one I found manually with Manipulate).

Does anybody know how I can improve FindFit or NonLinearModelFit? The documentation is not explanatory at all.