Have you tried implementing the PositiveDomain() in DiffEqCallbacks.jl, it was suggested to me for a similar problem, and it appears faster than the “manual” callback-integrator.
However, I noted that my system is indeed keep pushing to negative, and while the u stays close to 0, there a small negative perturbations that cause other u to exceed other physical limits (e.g., exceed mass conservation). I have not tried yet the combinations mentioned earlier by @Dan and @stevengj.