Dynamical Systems in Biology

Exercise 1:

The ODEs of the FitzHugh-Nagumo model, \(\dot{\vec x} = f(\vec x)\), are given by: \[ \begin{align} \dot{v}(t) & = v(t)\cdot \left( a-v(t)\right)\cdot\left(v(t)-1\right)-w(t)+I_{\mathrm{app}}, \\ \dot{w}(t) & = \epsilon\cdot\left(v(t)-\gamma\cdot w(t)\right) \end{align} \] where \(\mathbb{R} \ni t \mapsto v(t),w(t) \in \mathbb{R}\) and the parameters \((a, \epsilon, \gamma) \in \mathbb{R}_{+}^3\) are positive.

Write a function FN <- function(time, ini, pars) that can be passed as func argument to ode() of the package deSolve. Add two functions out2config <- function(out, parameter) and out2phase <- function(out, parameter) that arrange the output of ode() in a config or phase space data.frame.

# load packages needed for the exercises
library(deSolve)
library(ggplot2)

theme_set(theme_bw())

## FitzHugh-Nagumo ODE
FN <- function(time, ini, pars) {
  with(as.list(c(ini, pars)), {
  
    dv <- v * (a-v) * (v-1) - w + Iapp
    dw <- e * (v - g*w)
    
    return(list(cbind(dv, dw)))
  })  
}

## Write result of ode() in config space data.frame
out2config <- function(out, parameter="none") {
  ncols <- dim(out)[2]
  nrows <- dim(out)[1]
  data.frame(time = out[,1], 
             value = as.vector(out[,2:ncols]), 
             species = rep(colnames(out)[-1], each=nrows), 
             parameter = parameter)
}

## Write result of ode() in phase space data.frame 
out2phase <- function(out, parameter="none") {
  cbind(as.data.frame(out[,-1]), parameter=parameter)
}

Use ode() to integrate the system for \(a = 0.25\), \(\epsilon=0.002\), \(\gamma = 1.1\) and \(I_{\mathrm{app}}=0\) for different initial values \((v(0), w(0))\) varying \(v(0)\) around \(a = 0.25\) and keeping \(w(0) = 0\).
- Plot solutions in configuration space.
- Plot solutions in phase space.

# Set parameter values
Iapp <- 0

pars <- c(
  a = .25,
  e = .002,
  g = 1.1,
  Iapp = Iapp)

# Set vector of timepoints
times <- seq(0, sqrt(1000), len=2000)^2

# List of initial v values
v0 <- c(0.1, 0.245, 0.255, 0.3, 0.8)

# Loop over initial values and integrate FN ODEs
out <- lapply(v0, function(myv) {
  yini <- c(
    v = myv,
    w = 0
  )
  
  result <- ode(yini, times, FN, pars)
  return(result)
})

# apply out2config to each element of out and rbind them
outconfig <- do.call(rbind, lapply(1:length(out), function(i) out2config(out[[i]], paste0("v0=", v0[i]))))

# apply out2phase to each element of out and rbind them
outphase <- do.call(rbind, lapply(1:length(out), function(i) out2phase(out[[i]], paste0("v0=", v0[i]))))

## Plot results ---------------------------------

P1 <- ggplot(outconfig, aes(x=time, y=value, color=parameter)) + 
  facet_wrap(~species, scales="free", ncol = 1) + 
  geom_line()

print(P1)

P2 <- ggplot(outphase, aes(x=v, y=w)) + facet_wrap(~parameter) + 
  geom_path() 

print(P2)

Write a function FNnc <- function(v, pars) returning the nullclines of the system \(w(v)\) for \(\dot{v}(t) = 0\) and \(w(v)\) for \(\dot{w}(t) = 0\). Add them to your phase-space plot.

## FitzHugh-Nagumo nullclines
FNnc <- function(v, pars) {
  with(as.list(pars), {
    
    w_dveq0 <- v * (a-v) * (v-1) + Iapp
    w_dweq0 <- v/g
    
    result <- cbind(v, w_dveq0, w_dweq0)
    colnames(result) <- c("v", "dv = 0", "dw = 0")  
    
    return(result)
  })  
  
}

## calculate nullclines 
v <- seq(-0.2, 1, len=200)
nc <- FNnc(v, pars)

 P3 <- P2 + annotate("line", x=nc[,1], y=nc[,2], lty=2, col="red") + 
  annotate("line", x=nc[,1], y=nc[,3], lty=2, col="red") +
   coord_cartesian(ylim = c(-.05, .25))
 
 print(P3)

Now test \(I_{app}>0\). Plot solutions in config space for Iapp = c(0, 0.02, 0.1, 0.5, 0.7) and use the nullclines to understand the behaviour in phase space. The initial value \((v(0), w(0))\) can be fixed to \((0, 0)\).

# Set vector of timepoints
times <- seq(0, sqrt(5000), len=2000)^2

# List of initial v values
yini <- c(v = 0, w = 0)
Iapp <- c(0, .02, .1, .5, .7)

# Loop over initial values and integrate FN ODEs
out <- lapply(Iapp, function(I) {
  
  pars["Iapp"] <- I
  
  result <- ode(yini, times, FN, pars)
  return(result)
})

# apply out2config to each element of out and rbind them
outconfig <- do.call(rbind, lapply(1:length(out), function(i) out2config(out[[i]], paste0("Iapp = ", Iapp[i]))))

# apply out2phase to each element of out and rbind them
outphase <- do.call(rbind, lapply(1:length(out), function(i) out2phase(out[[i]], paste0("Iapp = ", Iapp[i]))))

# compute null-clines (different for each Iapp. Therefore,
# an output with the structure of the FN-solution is generated
# for later plotting)

outnc <- do.call(rbind, lapply(Iapp, function(I) {
  
  pars["Iapp"] <- I
  
  v <- seq(-0.5, 2, len=200)
  nc <- FNnc(v, pars) 
  
  result <- data.frame(v = v, w = as.vector(nc[, -1]), 
                       nullcline = rep(c("v", "w"), each = length(v)),
                       parameter = paste("Iapp =", I))
  
  return(result)
}))



## Plot results ---------------------------------

P1 <- ggplot(outconfig, aes(x=time, y=value)) + 
  facet_grid(parameter~species, scales="free") + 
  geom_line()

print(P1)

P2 <- ggplot(outphase, aes(x=v, y=w)) + facet_wrap(~parameter, scales = "free") + 
  geom_path() + 
  geom_line(data = subset(outnc, nullcline == "v"), color = "red", lty = 2) +
  geom_line(data = subset(outnc, nullcline == "w"), color = "red", lty = 2) +
  ylim(c(-.5, NA))
   

print(P2)

## Warning: Removed 34 rows containing missing values (geom_path).

Dynamical Systems in Biology

Tutorial 4

Daniel Kaschek, Mirjam Fehling-Kaschek

Summer term 2017

Exercise 1:

Cathedral exercise