

Hello Thomas,

I have now a working example for integrating Bauch et. al. model using DeSolve without it returning
worrisome warning messages.

I indeed use lsodar -  I do not use the events option but do use the root finding to stop integration at a 
particular event. What I do is integrate the model step wise (changing the vaccination rate nu) for each step.

You can find my code here below and - should it be viewed as a nice example - please by all means
do include it as a worked example.

Thank you.

Raffaello.



  
remove(list=ls())
setwd("C:/Documents and Settings/rvardava/My Documents/Projects_2009/Infectious_disease_Deductive_Games")
library(deSolve)


#############################################################################
##          SEIDV Epidemic Model
#############################################################################

## Bauch et. al. PNAS 2003 smallpox vaccination model
SEIDV_model <- function(t, Omega, parameters) {
 with(as.list(c(Omega,parameters)), {
   dV =  nu
   dS = -beta*S*I-dV
   dE =  beta*S*I-sigma*E
   dI =  sigma*E-gamma*I
   dD =  gamma*I
   return(list(c(dS, dE, dI, dD, dV)))
   })
   }
   
# event triggered if state variable S<=0 
# in this case we stop integrating and carry on with a different 
# vaccination rate nu.
rootfun <- function (t,  Omega, parameters) { 
with(as.list(Omega),{
  ifelse(S<=0, r.event <- 0 , r.event  <- 1 )
return(r.event )})}

####################################################################

EvolveSEIDVmodel <- function(p, final.time, t.step, parameters){
with(as.list(c(ode.parameters, pop.parameters)),{


  ## First phase: integration during the period the initial phase of a smallpox outbreak. 
  ## The outbreak is detected but vaccines are not mobilized yet and thus unavailable.
  times <- seq(from=0,to=tres,by=0.1) 
  Omega1 <- c(S = (1-p)*(N-alpha), E = (1-p)*alpha, I = 0, D = 0,V = 0)/N
  ode.parameters <- c(beta = beta, gamma = gamma, sigma = sigma, nu = 0, tres = tres)  
  Omega1 <- as.data.frame(lsodar(Omega1, times, func = SEIDV_model, parms = ode.parameters, hmax = t.step))
  
  ## Second phase: integration during the period of vaccination.
  l<-nrow(Omega1) 
  Omega2<-as.vector(as.numeric(Omega1[l,c(2:6)]))
  names(Omega2)<-names(Omega1[c(2:6)])
  ode.parameters <- c(beta = beta, gamma = gamma, sigma = sigma, nu = nu, tres = tres)
  
  times <-  seq(from=tres,to=final.time,by=0.1)  
  Omega2 <- as.data.frame(lsodar(Omega2 , times, func = SEIDV_model, parms = ode.parameters, 
              rootfun = rootfun, hmax = t.step)) 
  
  ## Third phase: integration during end phase. All the initially susceptible population has either been 
  ## infected or has been vaccinated.            
  l<-nrow(Omega2) 
  Omega3<-as.vector(as.numeric(Omega2[l,c(2:6)]))
  names(Omega3)<-names(Omega2[c(2:6)])
  if(Omega3["S"]<0) Omega3["S"] <- 0 
  ode.parameters <- c(beta = beta, gamma = gamma, sigma = sigma, nu = 0, tres = tres)
  
  times <-  seq(from=Omega2[l,"time"],to=final.time,by=0.1)  
  Omega3 <- as.data.frame(lsodar(Omega3 , times, func = SEIDV_model, parms = ode.parameters, hmax = t.step)) 
  
  Omega<- rbind(Omega1,Omega2,Omega3)
  
  ## remove duplicated points that occurred by rbinding the final state of one phase and the initial state
  ## of the successive phase.
  Omega <- Omega[!duplicated(Omega$time,fromLast = T),]             
  
  return(Omega)

})
}

####################################################################

## Parameters are those used by Bauch et. al. PNAS 2003
pop.parameters <- c(N = 2.85*10^8, alpha = 5*10^3, d_v = 10^(-6), d_s = 0.3 , r = 0.01)
ode.parameters <- c(beta = 5/3, gamma = 1/3, sigma = 1/11, nu = 1/10, tres = 14)

parameters <- c(pop.parameters, ode.parameters)

  ## p represents the pre-emptice vaccination coverage
  p <- 0.315
  final.time <- 25
  t.step <- 0.01
  
  Omega <- EvolveSEIDVmodel(p, final.time, t.step, parameters)
  
  Omega





> From: r-sig-dynamic-models-request@r-project.org
> Subject: R-sig-dynamic-models Digest, Vol 4, Issue 2
> To: r-sig-dynamic-models@r-project.org
> Date: Fri, 8 Jan 2010 12:00:02 +0100
> 
> Send R-sig-dynamic-models mailing list submissions to
> 	r-sig-dynamic-models@r-project.org
> 
> To subscribe or unsubscribe via the World Wide Web, visit
> 	https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> or, via email, send a message with subject or body 'help' to
> 	r-sig-dynamic-models-request@r-project.org
> 
> You can reach the person managing the list at
> 	r-sig-dynamic-models-owner@r-project.org
> 
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of R-sig-dynamic-models digest..."
> 
> 
> Today's Topics:
> 
>    1. Re: System with an event leading to a discontinuity in the
>       ODE equations: fear of non convergence and (Thomas Petzoldt)
> 
> 
> ----------------------------------------------------------------------
> 
> Message: 1
> Date: Fri, 08 Jan 2010 09:48:25 +0100
> From: Thomas Petzoldt <Thomas.Petzoldt@tu-dresden.de>
> Subject: Re: [R-sig-dyn-mod] System with an event leading to a
> 	discontinuity in the ODE equations: fear of non convergence and
> To: Special Interest Group for Dynamic Simulation Models in R
> 	<r-sig-dynamic-models@r-project.org>
> Message-ID: <4B46F159.6000600@tu-dresden.de>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
> I wondered a little bit why you use:
> 
> dV = nu
> 
> and not something like:
> 
> dV = nu * S
> 
> or another more smooth function to stop vaccination when S goes to zero.
> However, after googling, I found the paper of Bauch et al (2003), PNAS 
> 100(18), 10564-10567 that used a constant rate nu and a hard stopping 
> function f(S,t) because it's controlled by external intervention. Seems
> that this is indeed an example for combining root finding and events. 
> Note also that there is an initial response time t_res in that paper.
> 
> Raffaelo, if you have a working solution (e.g. using Karlines approach),
> it would be nice if we can include a full working example in our model 
> collection.
> 
> Thomas P.
> 
> 
> 
> ------------------------------
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> 
> 
> End of R-sig-dynamic-models Digest, Vol 4, Issue 2
> **************************************************
 		 	   		  
_________________________________________________________________

ail you.

cial-network-basics.aspx?ocid=PID23461::T:WLMTAGL:ON:WL:en-xm:SI_SB_3:092
	[[alternative HTML version deleted]]

