[R-sig-dyn-mod] from exponetial dwelling times to gamma dwellingtimes

[Soetaert, Karline]

Raffaello,
 
What you are trying to implement is a partial differential equation. 
In the environmental sciences we use such models when we describe processes occurring vertically in sediments, or how constituents change along a river, and so on. 
 
I had a look at this paper, and their equation 3 is really an advective equation in disguise:
 
It says for the dynamics of I_2:
dI_2/dt = n*ga*I_1 - n*ga*I_2 + mu.
 
Rearranging:  
dI_2/dt = -n*ga*(I_2 - I_1) + mu.
 
Defining dx = 1/n :
dI_2/dt = -ga*(I_2 - I_1)/dx + mu. 
 
or, for all I:
dI/dt =  -ga * dI/dx + mu
 
This is an advection-reaction equation.
 
 
Here "dx" is the "grid" size, and the grid goes from 0 to 1. "ga" is the "velocity" at which the individuals propagate through the infected states.
 
 
It is relatively simple and very efficient to implement numerical differences in R using the "diff" operator where 
d ga*I/dx is represented for the internal cells as: diff(ga*I)/dx.  Here ga*I is the "flux". You also need to specify what happens in the "upstream" compartment, where the flux is given by beta*I*S.
.
So your model would look more or less like:
 
myfunc <- function (t, Y, p, n) {
 dx <- 1/n
 S <- Y[1}] 
 I <- Y[2:(n+1)]  
 Influx <- beta * S* sum(I)*dx
 dI <- diff(c(Influx,ga*I))/dx + mu
 dS <- ....
 
 list(c(dS, dI))
}
 
 
This is also how the "advection" equation has been discretised in the paper you cite.
One note: this kind of discretisation has a rather large numerical error, called numerical diffusion , and that depends a.o. on the grid size, so when playing with "n" you will see smaller and larger values of these numerical errors.  
 
 
Also, there is an R-package called ReacTran which deals with "reactive transport" equations. The examples in the vignette are from the environmental sciences, but maybe you can get some inspiration there.
 
 
Hope this helps,
 
Karline

________________________________

From: r-sig-dynamic-models-bounces at r-project.org on behalf of Raffaello Vardavas
Sent: Tue 8/24/2010 10:52 PM
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] from exponetial dwelling times to gamma dwellingtimes




Dear All,

I have a general question. Lets say I'm constructing an SIR type infectious disease model. Although deterministically set, the differential equations assume that  the dwelling time in each compartment is exponential.
I.e., it doesn't matter if you were first or last to have been infected - the chances of progressing to the recovered stage is the same - thus there is the implied markov assumption and thus exponential distributed waiting times.

However biologically - other waiting/progression time distributions  (eg. gamma or weibull) make more biological sense.  Indeed, when modeling the SIR dynamics via an individual-based model - then  especially at disease invasion stage (i.e., far from equilibrium) the dynamics can be very different if one that uses a gamma distribution instead of an exponential distribution for progression rates.

One way to fix this in a SIR ODE model would be to use many concatenated similar compartments: So for the infected compartment we could split it up into n sequential compartments, I_1, I_2,..I_n - and use n sets of similar differential equations for I where all the rates in each single equation get multiplied by n.  This would lead to an implied gamma distribution for the dwelling or progression times in the infected stage - which as mentioned is biologically more sensible.

Question: is there anyway to do this in deSolve? Are there any example?

What I'm looking for is some function that could do this automatically:  splits each inputted ODE into many similar sequential ones with the scaled up rates - then integrates the model  - lastly sums back the dynamics for the sub-compartments (i.e., I_1, I_2,..I_n ) to obtain the dynamics of the overall infected compartment I. My guess is that this I would need to write my own function for this which takes as input an ODE and n.

Has anyone worked on something similar - and could given me an indication on how to start?
I basically want to avoid writing n sets of similar equations - also because I may want to experiment with n - as this sets the form of the beta distribution of the waiting times.

Thanks.

Raff

PS. I may not have been very clear in explaining myself. I therefore provide a reference which explains the problem:
Proc Biol Sci. 2001 May 7;268(1470):985-93. Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods.
Lloyd AL. 
                                         
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 8802 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20100826/e2aba46f/attachment.bin>