Error in R: Error in eval(substitute(expr), data, enclos = parent.frame()) : attempt to apply non-function

Question

I have been stuck on this code for months. Can't figure out what's wrong with the code here. Any help is appreciated. The code stops working after "time_values" line. I am new to R so please excuse the code.
I dont understand the error. I have a duplicate code with just lesser parameters, but that works fine, so no clue what is happening

library(deSolve)

sir_equations <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <-  b*N + E1*S*((a*P)/(1 + (a*h*P))) + E2*I*((a*P)/(1+(a*h*P))) + E3*S*((f*C)/(1+(f*h*C))) - v*S - beta*S*I - lambda*S*C + gamma*I
    dI <-  beta*I*S - gamma*I - (v+m)*I  + lambda*S*C - lambda*I*C
    dC <-  v*I + m*I - (S +I)((f*C)/(1 + f*h*C)) + lambda*I*C
    dP <-  v*S - (S + I)*((a*P)/(1 +(a*h*P)))
    dN <-  dS + dI
    return(list(c(dS, dI, dC, dP, dN)))
  })
}
parameters_values <- c(
  beta  = 0.005, 
  gamma = 0.2,    
  lambda = 0.1, 
  v = 0.09,  
  h = 0.3,  
  b = 0.01, 
  f = 0.4,  
  a = 0.6,  
  m = 0.9,   
  E1 = 1,  
  E2 = 0.1,  
  E3 = 0.2   
)

initial_values <- c(
  S = 99,  
  I =   1,  
  C =   0,
  N = 100,
  P = 99  
)

time_values <- seq(0, 50, 0.001) 

sir_values_1 <- ode(
  y = initial_values,
  times = time_values,
  func = sir_equations,
  parms = parameters_values 
)

sir_values_1 <- as.data.frame(sir_values_1)
sir_values_1

with(sir_values_1, {
  # plotting the time series:
  plot(time, S, type = "l", col = "blue",
       xlab = "time (days)", ylab = "number of people")
  # adding the time series of infectious:
  lines(time, I, col = "red")
  # adding the time series of recovered:
  lines(time, C, col = "purple")
  lines(time, P, col = "brown")
  lines(time, N, col = "black")
})

Answer 1

The problem is in the definition of sir_equations. In particular, the definition of dC is missing an asterisk, you should change it from

    dC <-  v*I + m*I - (S +I)((f*C)/(1 + f*h*C)) + lambda*I*C

to

    dC <-  v*I + m*I - (S +I)*((f*C)/(1 + f*h*C)) + lambda*I*C

note the asterisk after (S + I). After that, it should work:

sir_equations <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <-  b*N + E1*S*((a*P)/(1 + (a*h*P))) + E2*I*((a*P)/(1+(a*h*P))) + E3*S*((f*C)/(1+(f*h*C))) - v*S - beta*S*I - lambda*S*C + gamma*I
    dI <-  beta*I*S - gamma*I - (v+m)*I  + lambda*S*C - lambda*I*C
    dC <-  v*I + m*I - (S +I)*((f*C)/(1 + f*h*C)) + lambda*I*C
    dP <-  v*S - (S + I)*((a*P)/(1 +(a*h*P)))
    dN <-  dS + dI
    return(list(c(dS, dI, dC, dP, dN)))
  })
}
parameters_values <- c(
  beta  = 0.005, 
  gamma = 0.2,    
  lambda = 0.1, 
  v = 0.09,  
  h = 0.3,  
  b = 0.01, 
  f = 0.4,  
  a = 0.6,  
  m = 0.9,   
  E1 = 1,  
  E2 = 0.1,  
  E3 = 0.2   
)

initial_values <- c(
  S = 99,  
  I =   1,  
  C =   0,
  N = 100,
  P = 99  
)

time_values <- seq(0, 50, 0.001) 

sir_values_1 <- ode(
  y = initial_values,
  times = time_values,
  func = sir_equations,
  parms = parameters_values 
)

sir_values_1 <- as.data.frame(sir_values_1)
sir_values_1

with(sir_values_1, {
  # plotting the time series:
  plot(time, S, type = "l", col = "blue",
       xlab = "time (days)", ylab = "number of people")
  # adding the time series of infectious:
  lines(time, I, col = "red")
  # adding the time series of recovered:
  lines(time, C, col = "purple")
  lines(time, P, col = "brown")
  lines(time, N, col = "black")
})

You could make the plot a bit easier with ggplot2:

library(ggplot2)
library(tidyr)
sir_values_1 %>% 
  pivot_longer(-time, names_to="series", values_to="values") %>% 
  ggplot(aes(x=time, y=values, colour=series)) + 
  geom_line() + 
  theme_bw() + 
  theme(panel.grid=element_blank())