ODE solver problems with 'NonNegative' option

I am implementing an MCMC type fit on the parameters in an ODE model. The model represents biochemical species interactions, therefore it makes little sense to let the values go negative. However, when I implement the 'NonNegative' constraint, the ODE solver frequently fails. When I look at the failed parameter sets, the values get extremely small (i.e. 1e-47). First off, I am modeling molecular number, so a value less than 1 should simply be 0. Does anyone know a way I can designate integer steps for the ODE solver without moving to a stochastic set of equations or writing my own ODE solver? Or at least that the value be approximated as 0 below a certain tolerance? I have tried setting the absolute and relative tolerances quite high (i.e., abstol = 1) but I still get this error. Alternatively, I am unsure of the structure of my ODE model is suited to the solver I am using (ODE 15s). Per Matlab documentation: 'This option is not available for ode23s, ode15i, or for implicit solvers (ode15s, ode23t, ode23tb) applied to problems with a mass matrix.' Below I have pasted the function definition portion of my script. % direct activation model

g = @(t,x) [-deltaI*x(1) - kf*x(2)*x(1) + kb*x(3);                                                
             K - deltaN*x(2) - kf*x(2)*x(1) + kb*x(3);                                          
             -deltaN*x(3) + kf*x(2)*x(1) - kb*x(3);                                             
             alphaZ + (Vm1*x(3)^n1)/(Kp1 + x(3)^n1) - deltaZ*x(4);                              
             alphaL + (Vm2*x(4)^n2)/(Kp2 + x(4)^n2) - deltaL*x(5);];