Simulation5D.R

sim.png

Simulation5D.R

library(deSolve)
library(dplyr)
library(ggplot2)

# Model
f <- function(x) 258 / (1 + x ^ 7)
model <- function(t, state, parameters) {
  with(as.list(c(state)), {
    dx1 <- f(x5) - x1
    dx2 <- f(x1) - x2
    dx3 <- f(x2) - x3
    dx4 <- f(x3) - x4
    dx5 <- f(x4) - x5
    list(c(dx1, dx2, dx3, dx4, dx5))
  })
}

# Initial data
p1 <- c(x1 = 0.04013048, x2 = 257.972,  x3 = 0.01652379, x4 = 253.59240000, x5 = 2.60166400)
p2 <- c(x1 = 8.579741,   x2 = 6.763996, x3 = 3.692900,   x4 = 2.019579,     x5 = 1.336559)
times1 <- seq(0, 100, by = 0.1)
times2 <- seq(0, 13, by = 0.01)

# Simulation
out1 <- data.frame(lsoda(p1, times1, model, NULL))
out2 <- data.frame(lsoda(p2, times2, model, NULL))

# Projection
eigen2 <- c(-0.2628656, -0.4253254, 0.0000000, 0.4253254,  0.2628656)
eigen3 <- c(-0.3618034,  0.1381966, 0.4472136, 0.1381966, -0.3618034)
proj <- function(out, e) out$x1 * e[1] + out$x2 * e[2] + out$x3 * e[3] + out$x4 * e[4] + out$x5 * e[5]
t1 <- data.frame(e2 = proj(out1, eigen2), e3 = proj(out1, eigen3), name = "t1")
t2 <- data.frame(e2 = proj(out2, eigen2), e3 = proj(out2, eigen3), name = "t2")
df <- rbind(t1, t2)

# Plotting
plot2d <-
  df %>%
  ggplot(aes(x = e2, y = e3, group = name, colour = name)) +
  geom_path() +
  theme_bw() +
  ggtitle("Cycles in 5-dimensional systems") +
  scale_colour_discrete(name = "Trajectories",
                        labels = c("t1 (stable)\nValency=1", "t2 (unstable)\nValency=3"))
plot2d
ggsave("sim.png", plot2d)