[R-sig-dyn-mod] HELP for deSolve with SEIR model

[异昔时]

Dear colleagues, 

I¡¯m a new R-help user. I've read the advertisements about the good manners and I hope to propose a good question. 

I¡¯m using R to build an epidemiological SEIR model based on ODEs. The deSolve package is very useful to solve deterministic ODE systems but I'd like to perform a SEIR model with a 'time control variable 'on the beta parameter. 
           Because my study is in school, so every Saturday and Sunday is vacation. the value of beta (the rate of spread of the disease) is 0 on vacation. The value of beta is a piecewise function rather than a fixed value.
I don¡¯t know how to change the code of original R-CODE. Have you some suggestions about useful methods and/or function in R for reaching my aim. 

 

Thanks in advance. Jack 
 
 

 

library (deSolve)

seir_model = function (current_timepoint, state_values, parameters)

{

  # create state variables (local variables)

  S = state_values [1]        # susceptibles

  E = state_values [2]        # exposed

  I = state_values [3]        # infectious

  R = state_values [4]        # recovered

  

  with ( 

    as.list (parameters),     # variable names within parameters can be used 

         {

           # compute derivatives

           dS = (-beta * S * I)

           dE = (beta * S * I) - (delta * E)

           dI = (delta * E) - (gamma * I)

           dR = (gamma * I)

           

           # combine results

           results = c (dS, dE, dI, dR)

           list (results)

         }

    )

}

 

 

beta_value<-function(timepoints)              # beta=0.7 in original

{

result=ifelse(timepoints %in% c(6,7,12,14)==T,0,0.7)

return(result)

}

beta_value(timepoints)

 

gamma_value = 0.2

delta_value = 1 / 21

 

parameter_list = c (beta = beta_value(timepoints), gamma = gamma_value, delta = delta_value)

V = 5000        # exposed hosts

X = 9999        # susceptible hosts

Y = 1           # infectious hosts

Z = 0           # recovered hosts

 

N = V + X + Y + Z

initial_values = c (S = X/N, E = V/N, I = Y/N, R = Z/N)

timepoints = seq (0, 50, by=1)

output = lsoda (initial_values, timepoints, seir_model, parameter_list)

 

head(output,50)
	[[alternative HTML version deleted]]