Dynamical Systems in Biology

Exercise 1: matrix manipulation

There are a couple of pre-defined vectors and numbers in R, e.g., pi, letters, LETTERS, etc..

Create a 26 x 26 matrix C where each row contains all letters from A to Z. Do not use a loop!

C <- t(matrix(LETTERS, length(LETTERS), length(LETTERS)))
print(C)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [2,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [3,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [4,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [5,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [6,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [7,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [8,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [9,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [10,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [11,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [12,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [13,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [14,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [15,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [16,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [17,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [18,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [19,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [20,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [21,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [22,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [23,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [24,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [25,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [26,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
##  [1,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [2,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [3,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [4,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [5,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [6,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [7,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [8,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##  [9,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [10,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [11,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [12,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [13,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [14,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [15,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [16,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [17,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [18,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [19,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [20,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [21,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [22,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [23,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [24,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [25,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
## [26,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "V"   "W"   "X"  
##       [,25] [,26]
##  [1,] "Y"   "Z"  
##  [2,] "Y"   "Z"  
##  [3,] "Y"   "Z"  
##  [4,] "Y"   "Z"  
##  [5,] "Y"   "Z"  
##  [6,] "Y"   "Z"  
##  [7,] "Y"   "Z"  
##  [8,] "Y"   "Z"  
##  [9,] "Y"   "Z"  
## [10,] "Y"   "Z"  
## [11,] "Y"   "Z"  
## [12,] "Y"   "Z"  
## [13,] "Y"   "Z"  
## [14,] "Y"   "Z"  
## [15,] "Y"   "Z"  
## [16,] "Y"   "Z"  
## [17,] "Y"   "Z"  
## [18,] "Y"   "Z"  
## [19,] "Y"   "Z"  
## [20,] "Y"   "Z"  
## [21,] "Y"   "Z"  
## [22,] "Y"   "Z"  
## [23,] "Y"   "Z"  
## [24,] "Y"   "Z"  
## [25,] "Y"   "Z"  
## [26,] "Y"   "Z"

Change all “V”-entries to the square numbers 1 to 676. Determine by the == operator which row contains the number 100.

select <- which(C == "V")
C[select] <- (1:length(LETTERS))^2
print(C)

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [2,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [3,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [4,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [5,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [6,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [7,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [8,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##  [9,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [10,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [11,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [12,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [13,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [14,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [15,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [16,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [17,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [18,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [19,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [20,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [21,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [22,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [23,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [24,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [25,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
## [26,] "A"  "B"  "C"  "D"  "E"  "F"  "G"  "H"  "I"  "J"   "K"   "L"   "M"  
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
##  [1,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "1"   "W"   "X"  
##  [2,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "4"   "W"   "X"  
##  [3,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "9"   "W"   "X"  
##  [4,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "16"  "W"   "X"  
##  [5,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "25"  "W"   "X"  
##  [6,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "36"  "W"   "X"  
##  [7,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "49"  "W"   "X"  
##  [8,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "64"  "W"   "X"  
##  [9,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "81"  "W"   "X"  
## [10,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "100" "W"   "X"  
## [11,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "121" "W"   "X"  
## [12,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "144" "W"   "X"  
## [13,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "169" "W"   "X"  
## [14,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "196" "W"   "X"  
## [15,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "225" "W"   "X"  
## [16,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "256" "W"   "X"  
## [17,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "289" "W"   "X"  
## [18,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "324" "W"   "X"  
## [19,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "361" "W"   "X"  
## [20,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "400" "W"   "X"  
## [21,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "441" "W"   "X"  
## [22,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "484" "W"   "X"  
## [23,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "529" "W"   "X"  
## [24,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "576" "W"   "X"  
## [25,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "625" "W"   "X"  
## [26,] "N"   "O"   "P"   "Q"   "R"   "S"   "T"   "U"   "676" "W"   "X"  
##       [,25] [,26]
##  [1,] "Y"   "Z"  
##  [2,] "Y"   "Z"  
##  [3,] "Y"   "Z"  
##  [4,] "Y"   "Z"  
##  [5,] "Y"   "Z"  
##  [6,] "Y"   "Z"  
##  [7,] "Y"   "Z"  
##  [8,] "Y"   "Z"  
##  [9,] "Y"   "Z"  
## [10,] "Y"   "Z"  
## [11,] "Y"   "Z"  
## [12,] "Y"   "Z"  
## [13,] "Y"   "Z"  
## [14,] "Y"   "Z"  
## [15,] "Y"   "Z"  
## [16,] "Y"   "Z"  
## [17,] "Y"   "Z"  
## [18,] "Y"   "Z"  
## [19,] "Y"   "Z"  
## [20,] "Y"   "Z"  
## [21,] "Y"   "Z"  
## [22,] "Y"   "Z"  
## [23,] "Y"   "Z"  
## [24,] "Y"   "Z"  
## [25,] "Y"   "Z"  
## [26,] "Y"   "Z"

which(C[select] == "100") # 100 is a character!

## [1] 10

Restrict the matrix C to a 9 x 9 square matrix containing only the letters E, …, I and W, …, Z.

myind <- c("E", "F", "G", "H", "I", "W", "X", "Y", "Z")
rownames(C) <- colnames(C) <- LETTERS
C <- C[myind, myind]
print(C)

##   E   F   G   H   I   W   X   Y   Z  
## E "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## F "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## G "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## H "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## I "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## W "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## X "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## Y "E" "F" "G" "H" "I" "W" "X" "Y" "Z"
## Z "E" "F" "G" "H" "I" "W" "X" "Y" "Z"

Generate a new matrix E with entries TRUE and FALSE. The entries should be TRUE if the corresponding letter of the matrix C appears before K in the alphabet.

# The "intuitive" approach
Cnum <- match(C, LETTERS)
Knum <- match("K", LETTERS)
E <- matrix(Cnum < Knum, 9, 9)
print(E)

##       [,1] [,2] [,3] [,4] [,5]  [,6]  [,7]  [,8]  [,9]
##  [1,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [2,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [3,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [4,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [5,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [6,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [7,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [8,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
##  [9,] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE

# The elegant approach
E <- C<"K"
print(E)

##      E    F    G    H    I     W     X     Y     Z
## E TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## F TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## G TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## H TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## I TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## W TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## X TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## Y TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE
## Z TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE

Exercise 2: functions and loops

Write a function f <- function(n) which returns a vector for the integers 1,…,n with entries TRUE or FALSE for even or odd numbers, respectively.

Implement the function by looping over i=1,…,n. Have a look at is.integer() or the %% operator for testing if a number is odd or even (?is.integer, help("%%")).

myfunc <- function(n) {
  x <- 1:n
  iseven <- logical(length(x))
  for (i in 1:n)
    iseven[i] <- ceiling(x[i]/2) != ceiling((x[i]+1)/2)
  return(iseven)
}

myfunc(5)

## [1] FALSE  TRUE FALSE  TRUE FALSE

A key feature of R is vectorization, enabling the calculation of operations on vectors and matrices and making loops mostly unnecessary: implement the function without using a loop.

# Here x can be any vector of numbers
myfunc <- function(x) {
  iseven <- ceiling(x/2) != ceiling((x+1)/2)
  return(iseven)
}

myfunc(c(1, 4, 2, 6, 7))

## [1] FALSE  TRUE  TRUE  TRUE FALSE

Exercise 3: Newton method

Given the ODE for the Verhulst model \[\dot{x} = f(x,\vec{p}) = p_1 x \left(1-\frac{x}{p_2}\right)\] find the points \(\bar{x}\) where \(\dot{\bar{x}} = 0\) numerically.

Write a function fg <- function(x,p) returning the vector \(\big(f(x,\vec p), f'(x, \vec p)\big)\).

fg <- function(x, p) {
  fout1 <- p[1]*x*(1-x/p[2])
  fout2 <- p[1]-2*x*p[1]/p[2]
  return(c(fout1, fout2))
}

Implement the newton method \(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\) in a function newton <- function(x0,prec) with initial value x0 and precision prec. Stop the iterations when \(|x_{n+1} - x_n| <\) prec is reached.

newton <- function(x0, p = c(2, 3), precision = 1e-8, maxiter = 1000) {
  xn <- x0
  feval <- fg(xn, p)
  xnp1 <- xn-feval[1]/feval[2]
  for(i in 1:maxiter) {
    xn <- xnp1
    feval <- fg(xn, p)
    xnp1 <- xn-feval[1]/feval[2]
    if(abs(xn-xnp1) <= precision) break
  }

  return(xnp1)
}

newton(x0 = 5)

## [1] 3

For an illustration of Newton’s method see: https://en.wikipedia.org/wiki/Newton's_method#/media/File:NewtonIteration_Ani.gif

Dynamical Systems in Biology

Tutorial 0

Daniel Kaschek, Mirjam Fehling-Kaschek

Summer term 2017

Exercise 1: matrix manipulation

Exercise 2: functions and loops

Exercise 3: Newton method