The Lotka-Volterra model

LotkaVolterra <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    # ------------------------------------------------
    # Define the equations here
    # ------------------------------------------------
    killingRate       <- Prey * Predator
    preyGrowthRate    <- rGrowth * Prey
    predatorDeathRate <- rDeath * Predator
    dPreydt     <- preyGrowthRate - killingRate
    dPredatordt <- killingRate - predatorDeathRate
    return(list(c(dPreydt, dPredatordt)))
    # ------------------------------------------------
  })
}
pars  <- c(rGrowth  = 0.5,    # per day, growth rate of prey
           rDeath   = 0.2 )   # per day, death rate of predator
y0  <- c(Prey = 10, Predator = 2)
timespan <- seq(0, 200, by = 1)
LotkaVolterra.sol <- ode(y0, timespan, LotkaVolterra, pars)
summary(LotkaVolterra.sol)
## Plot the solution of the Lotka-Volterra model
plot(LotkaVolterra.sol)

The Birth/Death model

birthDeath <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    # ------------------------------------------------
    # Define the equations here
    # ------------------------------------------------
    deathRate         <- rDeath * x
    productionRate    <- rProduction
    proliferationRate <- rProlif * x
    dxdt     <- productionRate + proliferationRate - deathRate
    
    return(list(c(dxdt)))
    # ------------------------------------------------
  })
}
pars  <- c(rDeath      = 0.5,     # per day, death rate 
           rProduction = 0.2,     # individuals per day, production from outside source
           rProlif     = 0.1 )    # per day proliferation (or reproduction) rate
y0  <- c(x = 1.0)
timespan <- seq(0, 20, by = 0.1)
birthDeath.sol <- ode(y0, timespan, birthDeath, pars)
summary(birthDeath.sol)
plot(birthDeath.sol)

Negative feedback loop – Goodwin model

Goodwin <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    # ------------------------------------------------
    # Define the equations here
    # ------------------------------------------------
    mRNAproductionRate <- k0*K^h/(K^h + Z^h)
    mRNAdegradationRate    <- alpha * X
    dXdt <- mRNAproductionRate - mRNAdegradationRate
    dYdt <- beta * ( X - Y )
    dZdt <- beta * ( Y - Z )
    return(list(c(dXdt, dYdt, dZdt)))
    # ------------------------------------------------
  })
}
pars  <- c(k0      = 2,    # transcripts per hour, max mRNA synthesis rate
           alpha   = 1.0,  # per hour, mRNA degradation rate
           beta    = 1.0,  # per hour, kinetic rate
           K       = 1,    # mmol, half-repression concentration
           h       = 20 )  # no unit, Hill coefficient  
y0  <- c(X = 1, Y = 0, Z = 0)
timespan <- seq(0, 20, by = 0.1)
Goodwin.sol <- ode(y0, timespan, Goodwin, pars)
summary(Goodwin.sol)
plot(Goodwin.sol)

Positive feedback loop – Bistability

Bistability <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    # ------------------------------------------------
    # Define the equations here
    # ------------------------------------------------
    mRNAproductionRate <- k0*X^h/(K^h + X^h)
    mRNAdegradationRate    <- a * X
    dXdt <- mRNAproductionRate - mRNAdegradationRate
    return(list(c(dXdt)))
    # ------------------------------------------------
  })
}
pars  <- c(k0      = 2,    # transcripts per hour, max mRNA synthesis rate
           a       = 1.0,  # per hour, mRNA degradation rate
           K       = 1,    # mmol, half-repression concentration
           h       = 20 )  # no unit, Hill coefficient  
y0  <- c(X = 1.1)
timespan <- seq(0, 20, by = 0.1)
Bistability.sol1 <- ode( c(X = 0.9), timespan, Bistability, pars)
Bistability.sol2 <- ode( c(X = 1.1), timespan, Bistability, pars)
Bistability.sol3 <- ode( c(X = 1.0), timespan, Bistability, pars)
plot(Bistability.sol1, Bistability.sol2, Bistability.sol3)

Logistic map

LogisticMap <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    # ------------------------------------------------
    # Define the equations here
    # ------------------------------------------------
    xiter     <- r * x * ( 1 - x )
      
    return(list(c(xiter)))
    # ------------------------------------------------
  })
}
pars  <- c(r  = 3.76 )     # basal growth parameter 
y0  <- c(x = 0.2)
timespan <- seq(0, 50, by = 1)
LogisticMap.sol <- ode(y0, timespan, LogisticMap, pars, method = "iteration")
summary(LogisticMap.sol)
plot(LogisticMap.sol)

LS0tCnRpdGxlOiAiRHluYW1pY2FsIG1vZGVscyBpbiBTeXN0ZW1zIEJpb2xvZ3kgLS0tIFIgZXhhbXBsZXMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMgVGhlIExvdGthLVZvbHRlcnJhIG1vZGVsCgpgYGB7cn0KTG90a2FWb2x0ZXJyYSA8LSBmdW5jdGlvbihUaW1lLCBTdGF0ZSwgUGFycykgewogIHdpdGgoYXMubGlzdChjKFN0YXRlLCBQYXJzKSksIHsKICAgIAogICAgIyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KICAgICMgRGVmaW5lIHRoZSBlcXVhdGlvbnMgaGVyZQogICAgIyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KICAgIGtpbGxpbmdSYXRlICAgICAgIDwtIFByZXkgKiBQcmVkYXRvcgogICAgcHJleUdyb3d0aFJhdGUgICAgPC0gckdyb3d0aCAqIFByZXkKICAgIHByZWRhdG9yRGVhdGhSYXRlIDwtIHJEZWF0aCAqIFByZWRhdG9yCgogICAgZFByZXlkdCAgICAgPC0gcHJleUdyb3d0aFJhdGUgLSBraWxsaW5nUmF0ZQogICAgZFByZWRhdG9yZHQgPC0ga2lsbGluZ1JhdGUgLSBwcmVkYXRvckRlYXRoUmF0ZQoKICAgIHJldHVybihsaXN0KGMoZFByZXlkdCwgZFByZWRhdG9yZHQpKSkKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgfSkKfQoKcGFycyAgPC0gYyhyR3Jvd3RoICA9IDAuNSwgICAgIyBwZXIgZGF5LCBncm93dGggcmF0ZSBvZiBwcmV5CiAgICAgICAgICAgckRlYXRoICAgPSAwLjIgKSAgICMgcGVyIGRheSwgZGVhdGggcmF0ZSBvZiBwcmVkYXRvcgoKeTAgIDwtIGMoUHJleSA9IDEwLCBQcmVkYXRvciA9IDIpCnRpbWVzcGFuIDwtIHNlcSgwLCAyMDAsIGJ5ID0gMSkKTG90a2FWb2x0ZXJyYS5zb2wgPC0gb2RlKHkwLCB0aW1lc3BhbiwgTG90a2FWb2x0ZXJyYSwgcGFycykKc3VtbWFyeShMb3RrYVZvbHRlcnJhLnNvbCkKYGBgCgpgYGB7cn0KIyMgUGxvdCB0aGUgc29sdXRpb24gb2YgdGhlIExvdGthLVZvbHRlcnJhIG1vZGVsCnBsb3QoTG90a2FWb2x0ZXJyYS5zb2wpCmBgYAoKIyBUaGUgQmlydGgvRGVhdGggbW9kZWwKCmBgYHtyfQpiaXJ0aERlYXRoIDwtIGZ1bmN0aW9uKFRpbWUsIFN0YXRlLCBQYXJzKSB7CiAgd2l0aChhcy5saXN0KGMoU3RhdGUsIFBhcnMpKSwgewogICAgCiAgICAjIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQogICAgIyBEZWZpbmUgdGhlIGVxdWF0aW9ucyBoZXJlCiAgICAjIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQogICAgZGVhdGhSYXRlICAgICAgICAgPC0gckRlYXRoICogeAogICAgcHJvZHVjdGlvblJhdGUgICAgPC0gclByb2R1Y3Rpb24KICAgIHByb2xpZmVyYXRpb25SYXRlIDwtIHJQcm9saWYgKiB4CgogICAgZHhkdCAgICAgPC0gcHJvZHVjdGlvblJhdGUgKyBwcm9saWZlcmF0aW9uUmF0ZSAtIGRlYXRoUmF0ZQogICAgCiAgICByZXR1cm4obGlzdChjKGR4ZHQpKSkKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgfSkKfQoKcGFycyAgPC0gYyhyRGVhdGggICAgICA9IDAuNSwgICAgICMgcGVyIGRheSwgZGVhdGggcmF0ZSAKICAgICAgICAgICByUHJvZHVjdGlvbiA9IDAuMiwgICAgICMgaW5kaXZpZHVhbHMgcGVyIGRheSwgcHJvZHVjdGlvbiBmcm9tIG91dHNpZGUgc291cmNlCiAgICAgICAgICAgclByb2xpZiAgICAgPSAwLjEgKSAgICAjIHBlciBkYXkgcHJvbGlmZXJhdGlvbiAob3IgcmVwcm9kdWN0aW9uKSByYXRlCgp5MCAgPC0gYyh4ID0gMS4wKQp0aW1lc3BhbiA8LSBzZXEoMCwgMjAsIGJ5ID0gMC4xKQpiaXJ0aERlYXRoLnNvbCA8LSBvZGUoeTAsIHRpbWVzcGFuLCBiaXJ0aERlYXRoLCBwYXJzKQpzdW1tYXJ5KGJpcnRoRGVhdGguc29sKQoKYGBgCgpgYGB7cn0KcGxvdChiaXJ0aERlYXRoLnNvbCkKYGBgCgojIE5lZ2F0aXZlIGZlZWRiYWNrIGxvb3AgLS0gR29vZHdpbiBtb2RlbAoKYGBge3J9Ckdvb2R3aW4gPC0gZnVuY3Rpb24oVGltZSwgU3RhdGUsIFBhcnMpIHsKICB3aXRoKGFzLmxpc3QoYyhTdGF0ZSwgUGFycykpLCB7CiAgICAKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgICAjIERlZmluZSB0aGUgZXF1YXRpb25zIGhlcmUKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgICBtUk5BcHJvZHVjdGlvblJhdGUgPC0gazAqS15oLyhLXmggKyBaXmgpCiAgICBtUk5BZGVncmFkYXRpb25SYXRlICAgIDwtIGFscGhhICogWAoKICAgIGRYZHQgPC0gbVJOQXByb2R1Y3Rpb25SYXRlIC0gbVJOQWRlZ3JhZGF0aW9uUmF0ZQogICAgZFlkdCA8LSBiZXRhICogKCBYIC0gWSApCiAgICBkWmR0IDwtIGJldGEgKiAoIFkgLSBaICkKCiAgICByZXR1cm4obGlzdChjKGRYZHQsIGRZZHQsIGRaZHQpKSkKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgfSkKfQoKcGFycyAgPC0gYyhrMCAgICAgID0gMiwgICAgIyB0cmFuc2NyaXB0cyBwZXIgaG91ciwgbWF4IG1STkEgc3ludGhlc2lzIHJhdGUKICAgICAgICAgICBhbHBoYSAgID0gMS4wLCAgIyBwZXIgaG91ciwgbVJOQSBkZWdyYWRhdGlvbiByYXRlCiAgICAgICAgICAgYmV0YSAgICA9IDEuMCwgICMgcGVyIGhvdXIsIGtpbmV0aWMgcmF0ZQogICAgICAgICAgIEsgICAgICAgPSAxLCAgICAjIG1tb2wsIGhhbGYtcmVwcmVzc2lvbiBjb25jZW50cmF0aW9uCiAgICAgICAgICAgaCAgICAgICA9IDIwICkgICMgbm8gdW5pdCwgSGlsbCBjb2VmZmljaWVudCAgCgp5MCAgPC0gYyhYID0gMSwgWSA9IDAsIFogPSAwKQp0aW1lc3BhbiA8LSBzZXEoMCwgMjAsIGJ5ID0gMC4xKQpHb29kd2luLnNvbCA8LSBvZGUoeTAsIHRpbWVzcGFuLCBHb29kd2luLCBwYXJzKQpzdW1tYXJ5KEdvb2R3aW4uc29sKQpgYGAKCmBgYHtyfQpwbG90KEdvb2R3aW4uc29sKQpgYGAKCiMgUG9zaXRpdmUgZmVlZGJhY2sgbG9vcCAtLSBCaXN0YWJpbGl0eQoKYGBge3J9CkJpc3RhYmlsaXR5IDwtIGZ1bmN0aW9uKFRpbWUsIFN0YXRlLCBQYXJzKSB7CiAgd2l0aChhcy5saXN0KGMoU3RhdGUsIFBhcnMpKSwgewogICAgCiAgICAjIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQogICAgIyBEZWZpbmUgdGhlIGVxdWF0aW9ucyBoZXJlCiAgICAjIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQogICAgbVJOQXByb2R1Y3Rpb25SYXRlIDwtIGswKlheaC8oS15oICsgWF5oKQogICAgbVJOQWRlZ3JhZGF0aW9uUmF0ZSAgICA8LSBhICogWAoKICAgIGRYZHQgPC0gbVJOQXByb2R1Y3Rpb25SYXRlIC0gbVJOQWRlZ3JhZGF0aW9uUmF0ZQoKICAgIHJldHVybihsaXN0KGMoZFhkdCkpKQogICAgIyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KICB9KQp9CgpwYXJzICA8LSBjKGswICAgICAgPSAyLCAgICAjIHRyYW5zY3JpcHRzIHBlciBob3VyLCBtYXggbVJOQSBzeW50aGVzaXMgcmF0ZQogICAgICAgICAgIGEgICAgICAgPSAxLjAsICAjIHBlciBob3VyLCBtUk5BIGRlZ3JhZGF0aW9uIHJhdGUKICAgICAgICAgICBLICAgICAgID0gMSwgICAgIyBtbW9sLCBoYWxmLXJlcHJlc3Npb24gY29uY2VudHJhdGlvbgogICAgICAgICAgIGggICAgICAgPSAyMCApICAjIG5vIHVuaXQsIEhpbGwgY29lZmZpY2llbnQgIAoKeTAgIDwtIGMoWCA9IDEuMSkKdGltZXNwYW4gPC0gc2VxKDAsIDIwLCBieSA9IDAuMSkKQmlzdGFiaWxpdHkuc29sMSA8LSBvZGUoIGMoWCA9IDAuOSksIHRpbWVzcGFuLCBCaXN0YWJpbGl0eSwgcGFycykKQmlzdGFiaWxpdHkuc29sMiA8LSBvZGUoIGMoWCA9IDEuMSksIHRpbWVzcGFuLCBCaXN0YWJpbGl0eSwgcGFycykKQmlzdGFiaWxpdHkuc29sMyA8LSBvZGUoIGMoWCA9IDEuMCksIHRpbWVzcGFuLCBCaXN0YWJpbGl0eSwgcGFycykKYGBgCgpgYGB7cn0KcGxvdChCaXN0YWJpbGl0eS5zb2wxLCBCaXN0YWJpbGl0eS5zb2wyLCBCaXN0YWJpbGl0eS5zb2wzKQpgYGAKCiMgTG9naXN0aWMgbWFwCgpgYGB7cn0KTG9naXN0aWNNYXAgPC0gZnVuY3Rpb24oVGltZSwgU3RhdGUsIFBhcnMpIHsKICB3aXRoKGFzLmxpc3QoYyhTdGF0ZSwgUGFycykpLCB7CiAgICAKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgICAjIERlZmluZSB0aGUgZXF1YXRpb25zIGhlcmUKICAgICMgLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCiAgICB4aXRlciAgICAgPC0gciAqIHggKiAoIDEgLSB4ICkKICAgICAgCiAgICByZXR1cm4obGlzdChjKHhpdGVyKSkpCiAgICAjIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQogIH0pCn0KCnBhcnMgIDwtIGMociAgPSAzLjc2ICkgICAgICMgYmFzYWwgZ3Jvd3RoIHBhcmFtZXRlciAKCnkwICA8LSBjKHggPSAwLjIpCnRpbWVzcGFuIDwtIHNlcSgwLCA1MCwgYnkgPSAxKQpMb2dpc3RpY01hcC5zb2wgPC0gb2RlKHkwLCB0aW1lc3BhbiwgTG9naXN0aWNNYXAsIHBhcnMsIG1ldGhvZCA9ICJpdGVyYXRpb24iKQpzdW1tYXJ5KExvZ2lzdGljTWFwLnNvbCkKYGBgCgpgYGB7cn0KcGxvdChMb2dpc3RpY01hcC5zb2wpCmBgYAoK