avoid negative values when resolving a ODE

I am trying to model the behavior of a made-up networks of 5 genes, but I have the problem that I get negative values, which it has not sense biologically speaking.

Is there a way to limit the values to zero?

I managed to do it when I represent the graph, but I don't know how to use the ifelse in the main equation.

Thank you very much-1

###################################################
###preliminaries
###################################################

library(deSolve)
library(ggplot2)
library(reshape2)

###################################################
### Initial values
###################################################

values <- c(A = 1,
            B = 1,
            D = 1,
            E = 20,
            R = 1)

###################################################
### Set of constants
###################################################

constants <- c(a = 1.2,
               b = 0.5,
               c = 1.2,
               d = 1.5,
               e = 0.3,
               f = 0.5,
               g = 1.5,
               h = 0.9,
               i = 1.3,
               j = 1.3,
               m = 0.8,
               n = 0.6,
               q = 1,
               t = 0.0075,
               u = 0.0009,
               Pa = 100,
               Pb = 0.05,
               Pd = 0.1,
               Pe = 10)

###################################################
### differential equations
###################################################

Dynamic_Model<-function(t, values, constants) {
  with(as.list(c(values, constants)),{
    
    dA <- Pa + a*D - j*A - R
    dB <- Pb + b*A + e*E - m*B 
    dD <- Pd + d*B + f*E - g*A - n*D
    dE <- Pe - h*B + i*E - q*E
    dR <- t*A*B - u*D*E 
    
    list(c(dA, dB, dD, dE, dR))
  })   
}

###################################################
### time
###################################################

times <- seq(0, 200, by = 0.01)

###################################################
### print ## Ploting
###################################################

out <- ode(y = values, times = times, func = Dynamic_Model, parms = constantser)

out2 <- ifelse(out<0, 0, out)

out.df = as.data.frame(out2) 

out.m = melt(out.df, id.vars='time') 
p <- ggplot(out.m, aes(time, value, color = variable)) + geom_point(size=0.5) + ggtitle("Dynamic Model")