[R] Using odesolve to produce non-negative solutions

[Jeremy Goldhaber-Fiebert]

Hello,

I am using odesolve to simulate a group of people moving through time and transmitting infections to one another. 

In Matlab, there is a NonNegative option which tells the Matlab solver to keep the vector elements of the ODE solution non-negative at all times. What is the right way to do this in R?

Thanks,
Jeremy

P.S., Below is a simplified version of the code I use to try to do this, but I am not sure that it is theoretically right 

dynmodel <- function(t,y,p) 
{ 
## Initialize parameter values

	birth <- p$mybirth(t)
	death <- p$mydeath(t)
	recover <- p$myrecover
	beta <- p$mybeta
	vaxeff <- p$myvaxeff
	vaccinated <- p$myvax(t)

	vax <- vaxeff*vaccinated/100

## If the state currently has negative quantities (shouldn't have), then reset to reasonable values for computing meaningful derivatives

	for (i in 1:length(y)) {
		if (y[i]<0) {
			y[i] <- 0
		}
	}

	S <- y[1]
	I <- y[2]
	R <- y[3]
	N <- y[4]

	shat <- (birth*(1-vax)) - (death*S) - (beta*S*I/N)
	ihat <- (beta*S*I/N) - (death*I) - (recover*I)
	rhat <- (birth*(vax)) + (recover*I) - (death*R)

## Do we overshoot into negative space, if so shrink derivative to bring state to 0 
## then rescale the components that take the derivative negative

	if (shat+S<0) {
		shat_old <- shat
		shat <- -1*S
		scaled_transmission <- (shat/shat_old)*(beta*S*I/N)
		ihat <- scaled_transmission - (death*I) - (recover*I)
		
	}	
	if (ihat+I<0) {
		ihat_old <- ihat
		ihat <- -1*I
		scaled_recovery <- (ihat/ihat_old)*(recover*I)
		rhat <- scaled_recovery +(birth*(vax)) - (death*R)
	
	}	
	if (rhat+R<0) {
		rhat <- -1*R
	}	

	nhat <- shat + ihat + rhat

	if (nhat+N<0) {
		nhat <- -1*N	
	}	

## return derivatives

	list(c(shat,ihat,rhat,nhat),c(0))

}