r - Preventing a Gillespie SSA Stochastic Model From Running Negative -

i have produce stochastic model of infection (parasitic worm), using gillespie ssa. model used "gillespiessa"package (https://cran.r-project.org/web/packages/gillespiessa/index.html). in short code models population of discrete compartments. movement between compartments dependent on user defined rate equations. ssa algorithm acts calculate number of events produced each rate equation given timestep (tau) , updates population accordingly, process repeats given time point. problem is, number of events assumed poisson distributed (poisson(rate[i]*tau)), produces error when rate negative, including when population numbers become negative.

# parameter values  sir.parms <- c(deltahinfinity=0.00299, chi=0.00586, deltah0=0.0854, ah=0.5,                muh=0.02, sigmaw=0.1, sigmam =0.8, sigmal=104, phi=1.15, f = 0.6674,                deltavo=0.0166, cvo=0.0205, alphavo=0.5968, beta=52, mbeta=7300 ,muv=52, g=0.0096, n=100) # inital population values sir.x0 <- c(w=20,m=10,l=0.02) # rate equations sir.a <- c("((deltah0+deltahinfinity*chi*mbeta*l)/(1+chi*mbeta*l))*mbeta*l*n"            ,"sigmaw*w*n", "muh*w*n", "((1/2)*phi*f)*w*n", "sigmam*m*n", "muh*m*n",            "(deltavo/(1+cvo*m))*beta*m*n", "sigmal*l*n", "muv*l*n", "alphavo*m*l*n", "(ah/g)*l*n") # population change sir.nu <- matrix(c(+0.01,0,0,                    -0.01,0,0,                    -0.01,0,0,                    0,+0.01,0,                    0,-0.01,0,                    0,-0.01,0,                    0,0,+0.01/230,                    0,0,-0.01/230,                    0,0,-0.01/230,                    0,0,-0.01/230,                    0,0,-0.01/32),nrow=3,ncol=11,byrow=false) runs <- 10 set.seed(1)  # data frame of output sir.out <- data.frame(time=numeric(),w=numeric(),m=numeric(),l=numeric()) # multiple runs , combining data , ssa methods  for(i in 1:runs){   sim <- ssa(sir.x0,sir.a,sir.nu,sir.parms, method="etl", tau=1/12, tf=140, simname="sir")   sim.out <- data.frame(time=sim$data[,1],w=sim$data[,2],m=sim$data[,3],l=sim$data[,4])    sim.out$run <-   sir.out <- rbind(sir.out,sim.out) } 

thus, rates computed , model updates population values each time step, data store in data frame, attached previous runs. however, when levels of population low events can occur such number of events occurs reducing population greater number in compartment. 1 method make time step small, increases length of simulation long.

my question there way augment code data created/ calculated @ each time step values of population numbers negative converted 0?

i have tried working on problem, seem able come methods alter values once simulation complete, negative values still causing issues in runs themselves. e.g.
if (sir.out$l < 0) sir.out$l == 0

any appreciated

i believe problem method set ("etl") in ssa function. etl method produce negative numbers. can try "otl" method, based on efficient step size selection tau-leaping simulation method- in there few more parameters can tweak, basic command is:

ssa(sir.x0,sir.a,sir.nu,sir.parms, method="otl", tf=140, simname="sir") 

or direct method, not produce negative number whatsoever:

ssa(sir.x0,sir.a,sir.nu,sir.parms, method="d", tf=140, simname="sir")