Module in numerical model for NonlinearModelFit is slow and leaks memory

Question

I am trying to make a NonlinearModelFit where the model is the result of a numerical calculation. My code works, but it is very slow and also seems to have a memory leak. I'd appreciate advice about how to improve the code, make it faster, and stop leaking.

A simplified version of my model is:

$HistoryLength=0;
model[gnIN_?NumberQ, gpIN_?NumberQ] :=
   Module[{rep, eqs, sol, t},
    rep = {gn -> gnIN, gp -> gpIN, n0 -> 10, p0 -> 0.01};
    eqs = {
     n'[t] == -gn (n0 + n[t]) f[t],
     p'[t] == -gp ((p0 + p[t])) (1 - f[t]),
     n[t] == p[t] +f[t],
     n[0] == 10, p[0] == 10, f[0] == 0};
    sol = NDSolve[eqs /. rep, {n, p, f}, {t, 0, 10}];
    First[(n[#] + p[#])/(n[0] + p[0]) /. sol] &]

[Edit: maximum time in NDSolve changed from 50 to 10. See comment below.] There are three coupled equations (one of them not a differential equation, so it could be removed easily). The output that I want for my model, however, is not n, p, or f, but rather (n+p)/(n[0]+p[0]). I have used a pure function to pass this out of the model, but I don't know if there's a better way.

We can make some fake data by

tmp = model[1,.1];
data = Table[{i, tmp[i]}, {i, 0, 50}];

I then call

MemoryInUse[]/10^6.
nlm=NonlinearModelFit[data, {model[gn,gp][t], gn>=0 && gp>=0},{gn,gp},t];
MemoryInUse[]/10^6.

This runs successfully, giving best fit parameters {gn -> 1.04935, gp -> 0.0998232}, but it takes several minutes -- 376 seconds [Edit: revised to 68 seconds], from AbsoluteTiming on my 3.2GHz i5 iMac running MMA 9. (In this case, I could clearly give it excellent initial guesses, but that's a bit harder in my real problem. As an aside, giving it initial guesses of {2,0.2} causes it to take 450 seconds and return {1.48,0.098} as best fits, badly missing on gn.)

On this time that I'm running, the MemoryInUse commands return 92.9 and 101.1. Running it again will cause steady increase in memory usage, which is why I believe there is a memory leak.

Running Names["*$*"] returns a list with many variables like n$10849.

Names["n*$*"]//Length
Names["p*$*"]//Length

returns 289 and 336. I assume that all these n and p variables are being kept for some reason related to the return of the Module. But if I run the NonlinearModelFit again, the total memory continues to rise, but the number of n$* and p$* variables does not continue to increase -- it seems to stay very close to 300, indicating that some garbage collection is occurring. Also, the total memory used by those n$* and p$* variables, using symbolMemoryUsage from here, does not explain the increase in memory usage of the kernel. That is,

Total[symbolMemoryUsage /@ Names["*n$*"]]
Total[symbolMemoryUsage /@ Names["*p$*"]]

returns 13872 and 16176, which is much less than the 8 MB increase in memory use on running the NonlinearModelFit. In my actual case, the increased memory use can reach over 5 GB (and grind my computer to a halt) before NonlinearModelFit finishes. In this simple case, the leak is smaller, but it is still there. Running

MemoryInUse[]/10^6.
AbsoluteTiming[Do[model[1, .1], {i, 1, 1000}];]
MemoryInUse[]/10^6.

returns 104.161, 4.34, and 104.261. So the memory increase is about 100 KB with each 1000 evaluations of model. When I add a simple counter into model to count the number of times it is evaluated during the NonlinearModelFit, the answer is about 85000, indicating that the 8 MB increase in memory use is mostly due to a leak from model. Similarly, the time to evaluate model 85k times should be about 370 seconds, which is close to the total evaluation time. The slowness of my code and the memory leak are all in my model.

I understand that there's a limit to how fast a model with NDSolve could possibly be, but I suspect it can be better than this. Previous discussions of Block vs. Module have left me concerned that Block will do things I don't understand, but I haven't tried it. Any advice for coding practice (in particular, better ways to return the function from model so it's in the correct format to put into NonlinearModelFit), speed, or memory will be greatly appreciated.

EDIT: In case people have been dissuaded from working on this problem by the running time of my example code, I have changed the maximum time of NDSolve in model from 50 to 10, which makes the NonlinearModelFit run in 68 seconds, rather than 376. The memory leak is still there and of the same magnitude as in the original question. Any comments or advice welcome!