Dynamical Systems in Biology
Tutorial 5
Daniel Kaschek, Mirjam Fehling-Kaschek
Summer term 2017
Reaction-diffusion systems
Consider a reaction-diffusion system with two components \(u(x, t)\) and \(v(x, t)\): \[ \begin{align} \frac{\partial u(x,t)}{\partial t} & = \gamma \cdot f(u,v) + \frac{\partial^2 u(x,t)}{\partial x^2}, \\ \frac{\partial v(x,t)}{\partial t} & = \gamma \cdot g(u,v) + d \cdot \frac{\partial^2 v(x,t)}{\partial x^2} \end{align} \] At the boundary, \(u\) and \(v\) satisfy the zero-flux condition. We will simulate the system numerically in discrete space. Based on \(N\) equidistant supporting points, the distance between neighboring points is \(\Delta x = \frac{1}{N-1}\). The resulting approximation of the diffusion term reads
\[ \begin{align} \frac{\partial^2 u(x,t)}{\partial x^2} \quad \longrightarrow \quad \frac{u(x - \Delta x) - 2 u(x) + u(x + \Delta x)}{(\Delta x)^2}, \end{align} \] and analogously for \(v\). Thereby, the partial differential equation turns into a system of ordinary differential equations.
## Plotting functions
plotOut <- function(out, mode = c("time", "space")) {
t <- out[,1]
x <- 1:N
u <- out[,1:N+1]
v <- out[,(N+1):(2*N)+1]
mydata <- rbind(
data.frame(time = t, x = rep(x, each = length(t)), name = "u", value = as.vector(u)),
data.frame(time = t, x = rep(x, each = length(t)), name = "v", value = as.vector(v)))
mode <- match.arg(mode)
if (mode == "time") {
P <- ggplot(mydata, aes(x = time, y = value, group = x, color = x)) +
facet_wrap(~name, scales="free") +
geom_line()
}
if (mode == "space") {
P <- ggplot(mydata, aes(x = x, y = value, group = time, color = time)) +
facet_wrap(~name, scales="free") +
geom_line()
}
return(P)
}Exercise 1: Turing pattern formation
Implement the linear Turing system with the reaction functions
\[ \begin{align} f(u,v) & = a \cdot u + b \cdot v, \\ g(u,v) & = c \cdot u + e \cdot v, \end{align} \] with parameters \(a = -2\), \(b = 2.5\), \(c = -1.25\), \(\gamma = 1000\), \(d = 0.5\), \(e = 1.5\), and \(N = 250\). Let the initial values \(u_{1...N} \propto N(0,1)\) and \(v_{1...N} \propto N(0,1)\) be normally distributed random numbers.
# load packages needed for the exercises
library(deSolve)
library(ggplot2)
theme_set(theme_bw())
## Turing system
Turing <- function(Time, State, Pars, N) {
with(as.list(Pars), {
# Recover u and v from State vector
u <- State[1:N]
v <- State[(N+1):(2*N)]
dia <- c(-1, rep(-2, N-2), -1)
# Boundary conditions for u and v
uL <- c(u[-1], 0)
uR <- c(0, u[-length(u)])
vL <- c(v[-1], 0)
vR <- c(0, v[-length(v)])
# Reaction part of the system
f1 <- a*u + b*v
g1 <- c*u + e*v
# Reaction-diffusion system
deltax <- 1/(N-1)
du <- g * f1 + (uL + dia*u + uR)/deltax^2
dv <- g * g1 + d * (vL + dia*v + vR)/deltax^2
return(list(c(du, dv)))
})
}
# Set times and parameters and compute result
N <- 250
pars <- c(
a = -2,
b = 2.5,
c = -1.25,
d = 0.5,
e = 1.5,
g = 1000
)
yini <- rnorm(2*N, 0, 1)
times <- c(0, exp(seq(log(1e-3), log(1e1), length.out=200)))
out <- ode(method="lsodes", func = Turing, y = yini, parms = pars, times = times, N = N)
# Plot the result
plotOut(out, "time") + scale_x_log10()## Warning: Transformation introduced infinite values in continuous x-axis
plotOut(out, "space")
Exercise 2: Activator-inhibitor model of Gierer and Meinhardt
Implement the (dimensionless) activator-inhibitor model with the reaction functions
\[ \begin{align} f(u,v) & = a - b\cdot u + \frac{u^2}{v}, \\ g(u,v) & = u^2 - v, \end{align} \]
and parameters \(a = 0.1\), \(b = 1\), \(\gamma = 100\), \(d = 10\), and \(N = 250\). Use the homogeneous equilibrium states, \[ \begin{align} u_{1...N} = \frac{a + 1}{b} \quad \text{and} \quad v_{1...N} = \left( \frac{a + 1}{b} \right)^2, \end{align} \] and add noise to obtain initial values.
# Gierer-Meinhardt model
Gierer <- function(Time, State, Pars, N) {
with(as.list(Pars), {
# Recover u and v from State vector
u <- State[1:N]
v <- State[(N+1):(2*N)]
dia <- c(-1, rep(-2, N-2), -1)
# Boundary conditions for u and v
uL <- c(u[-1], 0)
uR <- c(0, u[-length(u)])
vL <- c(v[-1], 0)
vR <- c(0, v[-length(v)])
# Reaction part of the system
f2 <- a-b*u+(u^2)/v
g2 <- u^2 - v
# Reaction-diffusion system
deltax <- 1/(N-1)
du <- g * f2 + (uL + dia*u + uR)/deltax^2
dv <- g * g2 + d * (vL + dia*v + vR)/deltax^2
return(list(c(du, dv)))
})
}
# Set times and parameters and compute result
N <- 250
pars <- c(
a = .1,
b = 1,
g = 100,
d = 10
)
times <- c(0, exp(seq(log(1e-3), log(1e0), length.out=200)))
steady <- with(as.list(pars), c(u = (a+1)/b, v = (a+1)^2/b^2))
solveGierer <- function(n, gamma = 100, noise = 0) {
pars["g"] <- gamma
a <- pars["a"]
b <- pars["b"]
yini <- rep(c((a+1)/b, (a+1)^2/b^2), each = N)
if (n >= 0)
yini <- yini + c(0.1*cos(2*n*pi*c(1:N)/N), rep(0, each = N))
if (noise > 0)
yini <- yini + rnorm(length(yini), 0, noise)
out <- ode(method = "lsodes", func = Gierer, y = yini, parms = pars, times = times, N = N)
}- In how far does the result depend on the initial conditions?
set.seed(0)
out <- solveGierer(n = -1, noise = .1)
plotOut(out, "space")
out <- solveGierer(n = -1, noise = .1)
plotOut(out, "space")
out <- solveGierer(n = -1, noise = .1)
plotOut(out, "space")
out <- solveGierer(n = -1, noise = .1)
plotOut(out, "space")
- The size of the domain is explicitly represented by \(\gamma\). How does \(\gamma\) scale with the size of the domain, i.e., with the number of maximums of the activator.
set.seed(0)
out <- solveGierer(n = -1, noise = .1, gamma = 100)
plotOut(out, "space")
set.seed(0)
out <- solveGierer(n = -1, noise = .1, gamma = 400)
plotOut(out, "space")
set.seed(0)
out <- solveGierer(n = -1, noise = .1, gamma = 900)
plotOut(out, "space")
- Perturb the homogeneous equilibrium state by small modes of the form \(\cos(\frac{2n \pi x}{N})\), where \(x = \{1,...,N\}\) and \(n \in \mathbb{N}_{0}\). Check for which wave number \(n\) the equilibrium is instable and a pattern forms. To this end compute the dispersion relation \(\Re\{\lambda(k)\}\) where \(\lambda\) is an eigenvalue of \(J|_{(u_0, v_0)} - k^2 D\) and \(D\) is the diffusion matrix \(D = {\rm diag}(1, d)\).
library(numDeriv)
# Get the dispersion relation
dispersion <- function(n, steady, pars) {
# Reaction function
reaction <- function(state, pars) {
with(as.list(c(state, pars)), {
f <- a - b*u + u^2/v
g <- u^2 -v
return(c(f, g))
})
}
# Compute eigenvalues of J - k² D
with(as.list(pars), {
unlist(lapply(n, function(myn) {
M <- g*jacobian(reaction, x = steady, pars = pars) - (2*pi*myn)^2 * diag(c(1, d))
Re(eigen(M, only.values = TRUE)$values[2])
}))
})
}
n <- seq(0, 2, length.out = 200)
lambda <- dispersion(n, steady, pars)
roots <- n[diff(sign(lambda)) != 0]
qplot(x = n, y = lambda) + geom_hline(yintercept = 0) + geom_vline(xintercept = roots)
out <- solveGierer(n = 0, noise = 0, gamma = 100)
plotOut(out, "time") + scale_x_sqrt()
out <- solveGierer(n = 1, noise = 0, gamma = 100)
plotOut(out, "time") + scale_x_sqrt()
plotOut(out, "space")
out <- solveGierer(n = 2, noise = 0, gamma = 100)
plotOut(out, "time") + scale_x_sqrt()
plotOut(out, "space")
out <- solveGierer(n = 3, noise = 0, gamma = 100)
plotOut(out, "time") + scale_x_sqrt()
out <- solveGierer(n = 4, noise = 0, gamma = 100)
plotOut(out, "time") + scale_x_sqrt()
Cathedral exercise
What distinguished our cathedral from all other German gothic cathedrals?