Hi David, Thank you for your suggestion- I will look into it some more.
Thanks again, Hannah On Wed, Apr 18, 2018 at 5:09 PM, David Winsemius <dwinsem...@comcast.net> wrote: > > > On Apr 18, 2018, at 1:04 AM, Hannah Meredith <hrmeredit...@gmail.com> > wrote: > > > > Hello, > > > > I am solving a set of ODEs using deSolve and have run into a problem I > > would appreciate some advice on. One of the parameters (m) in the ODEs > > changes between two states when one of the variables (D) crosses a > > threshold (D_T) for the first time in either direction. Additionally, > when > > the variable crosses the threshold (either by increasing or decreasing), > > there is a time delay (delay) before the parameter shifts values. > > I think you migh not have fully reviewed the available functions in that > package. Doing an sos::findFn-search on ODE's with time delay brings up > help pages for the `dede` function in deSolve. > > For the conditional value of the within regime-fixed parameter `m`, > couldn't you just use this: > > ... + m*(6.8+ 3.2*(D < D_T)) > > Caveat: very lightly tested in your toy example, which did generate a > result where D cross that threshold without blowing up. > > -- > > David. > > > > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > Toy model: > > # define parameters > > P <- c(m = 10, > > R_1 = 0.5, > > R_2 = 0.1, > > D_0 = 50, > > D_T = 5, > > delay = 3 > > ) > > # set initial values > > y_0 <- c(D = 0, Y = 1) > > > > # set timesteps > > times <- seq(0, 100, length = 101) > > > > # define function > > bit <- function(t,y,parms){ > > with(as.list(c(parms,y)),{ > > > > > > # How I would implement the change if there was no time delay > > # m <- ifelse(D>=D_T, m*.68, m) > > > > # Option 1: nested if statements with time delay > > # if (D >= D_T) { > > # t_start1 <- t + delay > > # if (t >= t_start1){ > > # m <- 6.8 > > # } > > # } > > # if (D < D_T) { > > # t_start2<- t + delay > > # if (t >= t_start2){ > > # m <- 10 > > # } > > # } > > > > > > dD <- D_0 * (R_1 / (R_1 - R_2)) * (-R_2 * exp(-R_2 * t) + R_1 * > > exp(-R_1 * t)) > > dY <- Y + m; > > res <- c(dD, dY) > > list(res, m = m) > > }) > > } > > > > > > # Solve model > > library(deSolve) > > out <- ode(y = y_0, times = times, func = bit, parms = P) > > out.df <- as.data.frame(out) > > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > I have come up with a number of ideas, but I either haven't been able to > > get them to work or don't know if R has a command for it: > > > > 1. Nested "if" statements with something to track time (as shown above). > > One issue here is* t-start *is redefined with each time step in the ode > > solver so *t* never reaches *t_start +delay*. > > > > 2. Event triggered changes. From the examples I have found (link > > <https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf> > p. > > 24-25), the time at which the event happens is either pre-set or occurs > at > > a root of the ODE. Also- would this strategy require multiple events? > One > > event in which the variable crossing the threshold triggers a timer, > > another for the timer triggering the parameter change? > > > > 3. Lagvalue. This might take care of the time lag component, but I am > still > > struggling to understand (1) how lagvalue works and (2) how to integrate > > both the threshold and timer components. > > > > if ( D < D_T ) > > m<- ifelse(lagvalue(t-delay)<=0, 6.8, 10) > > end > > > > 4. Sigmoidal function. Someone suggested this to me earlier, but I am not > > sure how this would work. Is it possible to have a sigmoidal function > that > > switches repeatedly between two states (as opposed to the characteristic > S > > curve that transitions once from one state to another)? > > > > 5. Mathmematica has a WhenEvent function- is there something similar in > R? > > > > > > Thank you for your time and any tips you might have! > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide http://www.R-project.org/ > posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > David Winsemius > Alameda, CA, USA > > 'Any technology distinguishable from magic is insufficiently advanced.' > -Gehm's Corollary to Clarke's Third Law > > > > > > -- Hannah Meredith | PhD Whitaker Scholar London School of Hygiene and Tropical Medicine [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.