Programming Models in R

Question

I asked my original question incorrectly, so here is the better version.

I would like to generate a model using R. The gist of the model --> A polymer can grow or shrink at different rates. Every now and then the rate of shrinking increases by 20 fold while the growth rate stays constant, this is considered the "catastrophic state". The state switches to and from the "catastrophic state" with certain rates. The question then becomes how does the length of the polymer change with respect to time??? Here is what I am thinking:

Initialization:

5 polymers of length 0 (indicated by column index)

rnumbers <- data.frame(replicate(5,runif(20000, 0, 1)))

rstate <- rnumbers  # copy the structure
rstate[] <- NA      # preserve structure with NA's
# Init:
rstate[1, ] <- c(0)

> head(rstate)
  X1 X2 X3 X4 X5
1  0  0  0  0  0
2 NA NA NA NA NA
3 NA NA NA NA NA
4 NA NA NA NA NA
5 NA NA NA NA NA
6 NA NA NA NA NA

I want to run the simulation for 200s

Setting the rates :

dt <- c(.01) # how time is being divided
A <- dt*1 # probability growth will occur under normal conditions 
B <- dt*.25 # probability shrinking will occur under normal conditions
C <- dt*.03 # probability "catastrophic state" will occur
D <- dt*.3 # probability normal state will occur once in "catastrophic state"
E <- dt*5 # probability shrinking will occur "catastrophic state"

You notice that under normal conditions the probability of growth exceeds the probability of shrinkage, however when in "catastrophic state" shrinkage dominates. Also dt = .01 for 20000 rows in the data frame adds up to 200s

Not taking into account the switching to catastrophic state, this is what the code looks like :

library("Rcell")
rnumbers <- data.frame(replicate(5,runif(20000, 0, 1)))
dt <- c(.01)
A <- dt*1
B <- dt*.25
C <- dt*.03
D <- dt*.3
E <- dt*5

rstate <- rnumbers  # copy the structure
rstate[] <- NA      # preserve structure with NA's
# Init:
rstate[1, ] <- c(0)

step_generator <- function(col, rnum){
    for (i in 2:length(col) ){
            if( rnum[i] < B) { 
                col[i] <- -1
                }
                       else { 
                        if (rnum[i] < A) {
                            col[i] <- 1 
                            }
                              else {
                                col[i] <- 0 
                                }
                                }
                        }
    return(col)
    }
#  Run for each column index:
for(cl in 1:5){ rstate[ , cl] <- 
                        step_generator(rstate[,cl], rnumbers[,cl]) }

rstate1 <- transform(rstate, time = rep(dt))
rstate2 <- transform(rstate1, cumtime = cumsum(time))

cum_sum <- apply(rstate2, 2, cumsum)
cum_sum <- as.data.frame(cum_sum)

dev.new(width=5, height=5)  
cplot(cum_sum, X2 ~ time, geom = "line")

If you run this code, a jagged line with a positive slope is plotted over a 200 units of time. The package "Rcell" is required to use the plot I used.

My difficulty arises when I try to incorporate the catastrophic state. How can I make this code incorporate the catastrophic state? I imagine something like this, but I am unsure how to translate the syntax:

step_generator <- function(col, rnum)

for (i in 2:length(col)

start in normal growth conditions (growth prob = A ; shrinkage prob = B)

if rnum[i] < C switch to catastrophic state (
              growth prob = A ; 
              shrinkage prob = E), 
else stay in normal growth conditions (i.e. if rnum[i] >= C)

stay in catastrophic state until rnum[i] < D, when this happens switch back to normal growth conditions. 

repeat through the entire 20000 rows

thanks for any and all help!

Answer 1

Here is a solution to my question done on one column of random numbers, it can be easily iterated to multiple columns.

rnumbers <- data.frame(rnum = runif(20000, 0, 1)) #generates random #df 
dt <- c(.01) #time interval
A <- dt*1 #probability for growth
B <- dt*.25 # probability for shrinkage during growth phase
C <- dt*.03 # probability to enter catastrophe phase
D <- dt*.3 # probability to exit catastrophe phase and back into growth phase
E <- dt*5 # probability of shrinkage while in growth phase
rnumbers <- transform(rnumbers, cat = ifelse(rnum < C, .5, 0)) # generate column of times to enter catastrophe
rnumbers <- transform(rnumbers, rescue = ifelse(rnum < D, 1, 0)) # generate column of times to exit catastrophe
rnumbers <- transform(rnumbers, state = rowSums(rnumbers[,2:3])) # adds together two previous columns now (1.5 = catastrophe and 1 = growth and 0 = do nothing)
rnumbers[1,4] <- 1 #initialize in state 1 (growth phase)
rnumbers$state1 <- rnumbers$state # copies state into state 1
rnumbers$state1[rnumbers$state1 == 0] <- NA #makes any 0 in state1 NA
rnumbers$state1 <- na.locf(rnumbers$state1) # replaces NA with above row value
rnumbers <- transform(rnumbers, dynamics = ifelse(state1 < 1.5 , ifelse(rnum < B, -1, ifelse(rnum<A,1,0)), ifelse(rnum < A, 1, ifelse(rnum < E, -1,0))), interval = rep(dt)) # determines if polymer will grow or not
sum(rnumbers[,6])

rnumbers <- transform(rnumbers, time = cumsum(interval), step = cumsum(dynamics)) #cumulative sums of interval and growth or shrinkage
dev.new(width=5, height=5) #makes new plot
cplot(rnumbers,step ~ time, geom = "line") #plots length of polymer vs time

If anyone besides me is interested - this is a simulation of microtubule growth taking into account dynamic instability. The polymer starts at length initial (0 in this case), then can either grow or shrink at a certain probability. At any point during the simulation catastrophe can happen with a certain rate, this increases the probability of shrinkage by 20-fold and keeps growth probability the same. In order to get out of catastrophe rescue needs to occur which happens with a certain probability.

this is what polymer growth looks like with no catastrophe included :

library("Rcell")
require("zoo")
rnumbers <- data.frame(rnum = runif(20000, 0, 1))
dt <- c(.01)
A <- dt*1
B <- dt*.25
rnumbers <- transform(rnumbers, poly = ifelse(rnum < A, 1, ifelse(rnum < B, -1, 0)), dt = rep(dt))
rnumbers <- transform(rnumbers, time = cumsum(dt), step = cumsum(poly))
dev.new(width=5, height=5)
cplot(rnumbers,step ~ time, geom = "line")