Class 19: Comments on Midterm Exam

Computing for Molecular Biology 2

Andrés Aravena, PhD

21 May 2021

Question 1

1.1 How many steps?

Write a function, called num_steps that takes one input called N, representing the number of discs.

The function must return the number of steps required to move N discs from peg 𝐴 to peg 𝐵 following the game’s rules.

A true but wrong answer

num_steps <- function(N) {
(2^N-1)
}

How do you know this is a valid answer?

Can you give a mathematical proof of your answer?

You can only give this answer at Question 4

A little better

If N=0 number of steps=0
If N=2 number of steps=3
If N=3 number of steps=7
If N=4 number of steps=15
If N=5 number of steps=31
If N=6 number of steps=63
Formula of number of steps that depend on N is “(2^N)-1”

num_steps <- function(N) {
  if(N==0){
  return(0)
  }
  else{
  return((2^N)-1)}
}

a good answer

num_steps <- function(N) {
  if (N == 1) {
    return(1)
  } else {
    return(num_steps(1) + 2*num_steps(N-1))
  }
}

a good answer

num_steps <- function(N) {
  if(N <= 0) { 
    ans <- 0
  } else if(N == 1) {  
    ans <- 1
  } else {
    ans <- num_steps(N-1)*2 + 1 
  }
 return(ans)
}

not understanding

num_steps <- function(N) {
    # could not understand the question
}
  • If you don’t understand, it’s probably my fault

  • If you don’t ask, then it’s certainly your fault

Formulas must be proved

Can you believe this?

Formula for Fibonacci numbers. If \(F_n=F_{n1}+F_{n-2}\) and \(F_1=F_2=1,\) then \[F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n \right)\]

How do you know if this is true? Can you test it for \(F_{10^{100}}\)

Try this one

Let’s say we have \[1, 4, 8, 48, 88, 488, …\]

What comes next? Why?

Question 2

1.2 Fit a linear model

Now you can apply the num_steps function to the values of N from 1 to 30, and store the results in a vector named S.

You can plot S depending on N and you will see this diagram.

N <- 1:30
S <- sapply(N, num_steps)
plot(S~N)

This plot will help you to decide which transformation to use for fitting a linear model. Please build a good linear model and print the resulting coefficients.

Good answer

plot(log(S)~N)

model <- lm(log(S)~N)
exp(coef(model))

Must be log(S)

x <- lm(S ~ N, data=num_steps)

and print the coefficients

Be coherent with log and exp

model <- lm(N~S)
exp(coef(model))

What is this?

coef(S~N)

Question 3

1.3 Find the solution

Write a function, called hanoi(), taking four inputs:

  • N: the number of discs to move
  • origin: the name of the peg where the discs are now
  • destination: the name of the peg where the discs must be at the end
  • temporary: the name of the peg that can be used as a temporary place for some discs

The function must return a vector of text. Each element of the vector has the name of two pegs, indicating the movement of a single disc. The easiest way to make such text is using the paste() function.

A good answer

hanoi <- function(N, origin, destination, temporary) {
  if (N == 1) {
    return(paste(origin, destination))
  } else {
    return(c(paste(hanoi(N-1, origin, temporary, destination)),
             paste(hanoi(1, origin, destination, temporary)),
             paste(hanoi(N-1, temporary, destination, origin))))
  }
}

We don’t need paste(x). Just x

(but paste(x, y) is useful)

A better answer

hanoi <- function(N, origin, destination, temporary) {
  if (N == 1) {
    return(paste(origin, destination))
  } else {
    return(c(hanoi(N-1, origin, temporary, destination),
             hanoi(1, origin, destination, temporary),
             hanoi(N-1, temporary, destination, origin)))
  }
}

A simpler answer

hanoi <- function(N, origin, destination, temporary) {
  if (N == 1) {
    return(paste(origin, destination))
  } else {
    return(c(hanoi(N-1, origin, temporary, destination),
             paste(origin, destination),
             hanoi(N-1, temporary, destination, origin)))
  }
}

Almost good

hanoi <- function( N, origin, destination, temporary) { 
    if( N == 1) {
        print(paste(origin, destination))
    }
  else{
    hanoi(N-1, origin, temporary, destination )
    print(paste(origin, destination ) ) 
    hanoi(N-1, temporary, destination, origin)} 
  }

Do not print. Don’t tell secrets. Don’t break the contract.

Is this honest?

hanoi <- function(N, origin, destination, temporary) {
    toh_in <-  function(n, s, d) {
      if (n == 1) {
        x[i] <<- paste(temp[1:3 == s], temp[1:3 == d])
        i <<- i + 1
      }
      else {
        toh_in(n - 1, s, 6 - s - d)
        x[i] <<- paste(temp[1:3 == s], temp[1:3 == d])
        i <<- i + 1
        toh_in(n - 1, 6 - s - d, d)
      }  
      return(x[(i - 2 ^ n + 1):(i - 1)])
    }
    temp <<- c("A", "B", "C")
    s = which(temp == origin)
    d = which(temp == destination)
    x <<-  i <<-  1
    toh_in(n = N, s, d)
}

What is this?

hanoi <- function(N, origin, destination, temporary) {
  origin <- rep(NA, N)
  destination <- rep(NA, N)
  temporary <- rep(NA, N)
  origin[1] <- temporary
  destination[1] <- 0
  for(i in 2:N) {
    origin[i] <- destination[i-1]*temporary
    destination[i] <- destination [i-1] + origin[i]
  }
  return(data.frame(origin, destination))
}
paste(origin, destination)

What is this?

hanoi <- function(N, origin, destination, temporary) {
    origin  <- d_origin <- rep(NA, N)
    destination    <- d_destination   <- rep(NA, N)
    temporary     <- d_temporary <- rep(NA, N)
    origin[1] <- N
    d_origin <- 0
    destination[1] <- d_destination <- 0
    temporary[1] <- d_temporary <- 0
    
    for(i in 2:N) {
      d_origin[i] <- -temporary[i-1] * destination[N]
      d_destination[i] <- + temporary[i-1] * origin[N]
      d_temporary[i] <- + origin[i-1] * - destination[i-1]
      origin <- origin[i-1] + d_origin[i]
      destination <- destination[i-1] + d_destination[i]
      temporary <- temporary[i-1] + d_temporary[i]
    }
    paste(origin, destination, temporary)
}

What is this?

hanoi <- function(N, origin, destination, temporary) {
    v1 <- c(("A"-> "B"), ("A" -> "C"), ("C" -> "A"), ("B" -> "A"), ("C" -> "B"))
    return(paste(v1))
}

Question 4

1.4 (bonus) Find the formula

What is the formula for S depending on N? To get an exact result, you may use first a linear model to find the relation between S+1 and N. How long would it take to move 64 discs if you move 1 every second?

# write your answer here in plain English

Good answer, but is this English?

When we add 1 to S, our coefficient numbers become 1 and 2.

model2 <- lm(log(S+1)~N)
exp(coef(model2))

In every stage, we multiply the previous stage result, and add 1 extra step. So the formula should be something like this: 2 * hanoi(N-1) + 1.

For a 64 discs game, we need to move:

num_steps(64)

1.844674e+19 times. If we move 1 disc in every second, it will take 1.844674e+19 second. A year, basically has 3e+06 seconds.

num_steps(64) / 3e+06

It would take 6e+12 year, or 60 billion centuries.

Without Markdown rules

Finding Formula: 1. Plotting S depending on N. - We can use plot() function. 2. Building a good linear model. - We should use log() function on S and use lm() function. We should store the results. 3. Finding coefficients of the good linear model. - We should use exp() function for the stored result. We will get 2 values. 4. We should write y= First value * second value^N. - If we write to all of codes correct, we will get this formula: y=0.8536174*2.0151370^N

Moving 64 Discs:

1-We will find out how many moves we have to move 64 disks. -For this, we can use the function we wrote in the first question. We will get the result 1.844674e + 19. This number gives us the total number of discs we will move. 2-We multiply this number by 1 second.

With Markdown rules

Finding Formula:

  1. Plotting S depending on N.
    • We can use plot() function.
  2. Building a good linear model.
    • We should use log() function on S and use lm() function. We should store the results.
  3. Finding coefficients of the good linear model.
    • We should use exp() function for the stored result. We will get 2 values.
  4. We should write y= First value * second value^N.
    • If we write to all of codes correct, we will get this formula: y=0.8536174*2.0151370^N

With Markdown rules (cont…)

Moving 64 Discs:

  1. We will find out how many moves we have to move 64 disks.
  • For this, we can use the function we wrote in the first question. We will get the result 1.844674e + 19. This number gives us the total number of discs we will move.
  1. We multiply this number by 1 second.

Answer

I compared the two linear models with the two following codes using the plot:

lines((exp(predict(lm(log(S)~N)))~N))

lines((exp(predict(lm(log(S+1)~N)))~N))

When I compared both outputs, I concluded that the second linear model yields much more accurate results.

Thus, I calculated the coefficients of the second linear model, which are 1 and 2, respectively. Therefore, the formula for S+1 depending on N is:

S+1=1*2^N
S=2^N-1

It would take 2^64-1 (1.844674e+19) seconds to move 64 discs.

Answer

For Every N, it can be calculated as “2 ^ N - 1”. So S is calculated. To move 64 disks, we need to move the disks 1,845 e+19 times. This process takes more than 500 million years.

Answer

Formula could be written  Tn = 2^n-1. 
n = number of discs.
So;
T64 = 2^64-1
    = 1.85x10^19

Is this English?

Formulation of Hanoi moves S=(N^2)-1. We can explain if we want to find the number of moving(N) for depends on S, we calculate square root of S+1 for find to N value or we can improve a code for count move discs timing.

find_time<-function(S) {
  second<-S%%60
  S<-S%/%60
  minute<-S%%60
  S<-S%/%60
  hour<-S%%24
  S<-S%/%24
  day<-S%%30
  S<-S%/%30
  month<-S%%12
  S<-S%/%12
  year<-S
  print(paste("if every move take a second it will take ",year,"years ",month,"month", day, "days", hour, "hours",minute, "minutes", second,"seconds"))
}
find_time(num_steps(64))
remove(find_time)

The Legend

According to the legend, at the beginning of time, the Hindu temple priests were given a tower of 64 disks of gold and they had the task of transferring the disks from one pole to the third pole on the other side of the temple

The priests worked day and night to complete the task, and when they did, legend has it that the temple crumbled into dust, and the world vanished before the priests could complete the task

The legend was spread by Edouard Lucas in hopes of helping to popularize the Tower of Hanoi game. He invented the “Tower of Hanoi” puzzle in 1883

Question 5

2.1 Get Eschericia coli genome

Read E.coli genome (accession NC_000913) and store it in genome.

To verify, print the genome length

length(genome)

Use the correct index

library(seqinr)
read.fasta("NC_000913.ffn", seqtype = "DNA", set.attributes = FALSE)
genome <- read.fasta("NC_000913.ffn", seqtype = "DNA", set.attributes = FALSE)

Function read.fasta gives a list

We want the first list element

Answer

library(seqinr)
genes <- read.fasta("NC_000913.ffn", seqtype = "DNA", set.attributes = FALSE)
genome <- genes[[1]]

Here we are reading genes ("FASTA features nucleotides, .ffn), not chromosomes (Fasta nucleotides, .fna)

A good answer

library(seqinr)
chromosomes <- read.fasta("NC_000913.fna", seqtype="DNA", set.attributes = FALSE)
genome <- chromosomes[[1]]

Do not include full paths

library(seqinr)
chromosomes<- read.fasta("C:/Users/Hp2020/Downloads/NC_000913.fna", seqtype= "DNA",  set.attributes = FALSE)
genome <- chromosomes[[1]]

This does not work in other people’s computers

Do not install packages

install.packages("seqinr")
library(seqinr)
setwd("C:/Users/Casper/Desktop")
chromosome <- read.fasta("sequence.fasta" , seqtype = "DNA") 
genome <- chromosome[[1]]

Never use install.packages or setwd in a Rmarkdown file

Question 6

2.2 A function to score a DNA fragment

The score of a DNA fragment is calculated as the sum of the value of the matrix on for each letter on each position. That is why this kind of matrices are called “Position Specific Scoring Matrix”, or PSSM.

The first letter takes the value from the first row, the second letter looks into the second row, and so on. Then all the values are summed together, that is the score. For example, the sequence ACTGACTGA is \[-186 + -127 + -186 + -286 + -86 + 184 + -127 + -186 + 167\]

Please write an R function, called calc_score to evaluate the score of a DNA fragment. The function takes two inputs: the DNA fragment, and a Position Specific Scoring Matrix, and returns the score of that fragment. That is, the sum of all values in the corresponding row and column.

Answer

calc_score <- function(dna, M) {
  V <- toupper(dna)
  scores <- rep(NA, length(V))
  for(i in 1:length(V)) {
    if(V[i] == "A") {
      scores[i] <- M$A[i]
    } else if(V[i] == "T") {
      scores[i] <- M$T[i]
    } else if(V[i] == "G") {
      scores[i] <- M$G[i]
    } else if(V[i] == "C") {
      scores[i] <- M$C[i]
    }
  }
  ans <- sum(scores)
  return(ans)
}

Answer

calc_score <- function(dna, M) {
  n<-toupper(dna)
  toplam<-0
  b<-0
  for (i in n){
    b<-b+1
    if (i=="A"){
      toplam =toplam + M[1,b]
    }else if(i=="c"){
      toplam=toplam + M[2,b]
    }else if(i=="G"){
      toplam=toplam + M[3,b]
    }else if (i=="T"){
      toplam = toplam + M[4,b]
    }
  }
  return(toplam)
}

Answer

calc_score <- function(dna, M) {
  a <- toupper(dna)
  fragment_score <- rep(1:length(a))
  for(i in 1:length(a)){
    fragment_score[i] <- M[i,a[i]]
  }
  return(sum(fragment_score))
}

Answer

calc_score <- function(dna, M) {
    dna=toupper(paste(dna,collapse = ""))
  score = 0
  for (i in 1:nchar(dna)) {
    score = score + M[i, substring(dna, i, i)]
  }
  score
}

a good solution

calc_score <- function(dna, M) {
  dna <- toupper(dna)
  score <- 0
  for (i in 1:nchar(dna)) {
    score <- score + M[i, dna[i]]
  }
  return(score)
}

Question 7

2.3 Another function to score any genomic position

Using the calc_score function from the previous question, we will now calculate the score for any position in the genome. Please write a function called score_of_position that takes three inputs: the position in the genome (a number), the PSSM matrix, and the genome.

Answer

score_of_position <- function(i, M, genome) {
Sum <-0
position <- i
genome <- toupper(genome)
for (j in 1:9) {
if (genome[position +j-1] == "A"){
   Sum <- Sum+M[j,1]}
 if (genome[position +j-1] == "C"){
    Sum <- Sum+M[j,2]}
if (genome[position +j-1] == "G"){
    Sum <- Sum+M[j,3]}
if (genome [position +j-1] == "T"){
    Sum <- Sum+M[j,4]}}
   return(Sum)}

Answer

score_of_position <- function(i, M, genome) {
  x  <- toupper(genome[seq(from=i,length.out=9)])
  score_vector <- rep(NA, length(x))
    for(i in 1:length(x)) {
      if (x[i] == "T") {
             score_vector[i] <- M$T[i]
      } else if(x[i] == "G") {
             score_vector[i] <- M$G[i]
      } else if(x[i] == "C") {
             score_vector[i] <- M$C[i]
      } else if(x[i] == "A") {
             score_vector[i] <- M$A[i]
      }
  }
  ans <- sum(score_vector)
  return(ans)
}

Almost good

score_of_position <- function(i, M, genome) {
  position <- seq(from=i, length.out=9)
  window <- genome[position]
  answer <- calc_score(window, M)
  return(answer)
}

Why 9?

Almost good

score_of_position <- function(i, M, genome) {
  calc_score(genome[i:(i + 8)], M)
}

Why 8?

My solution

score_of_position <- function(i, M, genome) {
  calc_score(genome[seq(from=i, length.out=nrow(M))])
}

Now it works with any matrix M

Question 8

2.4 Calculate all scores

In class 9 we used the GC skew to locate a genomic region where the replication origin may be located. The binding site for DnaA should be in the same range of positions.

position <- 1540000:1550000

Now you can apply the score_of_position function to evaluate the score for each element of position, and store the resulting values in a new vector called score. This answer is just code, not a function.

A good classical solution

score <- rep(NA, length(position))
for(i in 1:length(position)){
  score[i] <- score_of_position(position[i],
                                PSSM_dnaA, genome)
}

A good modern solution

score <- sapply(position, score_of_position,
                PSSM_dnaA, genome)

What is the function name?

score_of_position <- seq(from=1540000,
                          to=1550000,
                          by=position)
sapply(score_of_position, score, position, genome)

Question 9

2.5 Where is the largest score?

Finally, to find the binding site for DnaA, we need to find the place that has the largest score. We do not care about the score, only about its position in the genome.

This works

position[which.max(score)]

This also works

position[which(score==max(score))]

This is wrong

max(score)

This is very wrong

max(score_of_position)
max

Bonus:

To verify that you read this questions at the right time, please send me an email to the answer address with your student number in the subject and the text “I read it”.

Question 10

3.1 Simulate the system

Please complete the code of the food_mice_cats function. The basic parts of the code have already been filled for you. Do not change the code outside the answer region.

My answer

birth <- birth_rate * fat_mice[i-1]
growth <- growth_rate * food[i-1]
reprod <- reprod_rate * fat_cats[i-1]
catch <- catch_rate * fat_mice[i-1] * cats[i-1]
death <- death_rate * fat_cats[i-1]
eating <- eating_rate * food[i-1] * mice[i-1]
hunger <- hunger_rate * mice[i-1]
dying <- dying_rate * cats[i-1]

d_food[i] <- growth - eating
d_mice[i] <- birth - eating - hunger
d_cats[i] <- reprod - catch - dying
d_fat_cats[i] <- catch - death
d_fat_mice[i] <- eating - catch 

Question 11

3.2 Effect of eating rate

In the previous question we saw that mice eat the food until the cats catch most of the fat_mice. One strategy that mice can use to survive is to eat faster so they can grow faster and outcompete the cats. To explore this question you will write a function, called final_state, that takes one input called rate.

This function should simulate the food_mice_cats system for 180 days using the same initial conditions and rates as in the previous question, except that the eating_rate input of food_mice_cats should be made equal to rate.

3.2 Effect of eating rate

The output of food_mice_cats should be a vector with four elements. The element names should be growth, food, total_cats, and total_mice. The value of growth should be equal to rate. The value of food must be the quantity of food on the last day of simulation (day 180). The value of total_cats must be the sum of the number of cats and fat_cats on the last day. In the same way, total_mice must be the sum of the number of mice and fat_mice on the last day.

Answer

final_state <- function(rate) {
  result <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
                reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
    growth <- rate
    food <- result$food[180]
    total_cats <- result$cats[180] + result$fat_cats[180]
    total_mice <- result$mice[180] + result$fat_mice[180]
 return(c(growth, food, total_cats, total_mice))
}

Answer

final_state <- function(rate) {
    result <-
    food_mice_cats(
      N = 180, birth_rate = 0.1, growth_rate = 1e-3,
      reprod_rate = 0.04, catch_rate = 0.05, death_rate = 0.01,
      eating_rate = rate, hunger_rate = 1e-3, dying_rate = 0.05,
      food_ini = 100, mice_ini = 10, fat_mice_ini = 0,
      cats_ini = 10, fat_cats_ini = 0)
  result <- tail(result, 1)[c("food", "cats", "fat_cats",
                              "mice", "fat_mice")]
  result <- c(rate, result[1, 1], sum(result[2:3]), sum(result[4:5]))
  names(result) <- c("growth", "food" , "total_cats" , "total_mice")
  result
}

Answer

final_state <- function(rate) {
food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3, reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0,
                cats_ini=10, fat_cats_ini=0) -> a
  
  answer_final_state<- c(growth = rate, food = a$food[180],
              total_cats = sum(a$cats[180] + a$fat_cats[180]),
              total_mice = sum(a$mice[180] + a$fat_mice[180]))
   return(answer_final_state)}

Answer

final_state <- function(rate) {
  result <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
                reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
  return(c(growth = rate, 
           food = result$food[180], 
           total_cats = result$cats[180]+result$fat_cats[180], 
           total_mice = result$mice[180]+result$fat_mice[180]))
}

Answer

final_state <- function(rate) {
        keeper  <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
                reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
  growth<-rate
  food<-keeper["food"][[1]][180]
  total_cats<-keeper["fat_cats"][[1]][180]+keeper["cats"][[1]][180]
  total_mice<-keeper["fat_mice"][[1]][180]+keeper["mice"][[1]][180]
  return(tibble(growth,food,total_cats,total_mice))
}

Answer

final_state <- function(rate) {
  for(i in 1:length(rate)){
    data <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
                reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate[i],
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)}
  
  out <- c(rate[i], data[180,2,i], data[180,4,i]+data[180,5,i],
          data[180,3,i]+data[180,6,i])
  return(out)
}

Answer

final_state <- function(rate) {
  n<-180
  newResult <- food_mice_cats(N=n, birth_rate=0.1, growth_rate=1e-3,
                  reprod_rate=0.04,
                  catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                  hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                  mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
  namer<-list(rate,newResult$food[n],
        (newResult$cats[n] + newResult$fat_cats[n]),
        (newResult$mice[n] + newResult$fat_mice[n]))
  names(namer[[1]])<-"growth"
  names(namer[[2]])<-"food"
  names(namer[[3]])<-"total_cats"
  names(namer[[4]])<-"total_mice"
  return(c(namer[[1]],namer[[2]],namer[[3]],namer[[4]]))
}

Bad Answer

final_state <- function(rate) {
    food_mice_cats <- c("growth"== rate,
                        "food" == food(180),
                        "total_cats" == sum(cats(180), fat_cats(180)), "total_mice" == sum(mice(180), fat_mice(180)))
    food_mice_cats %>% mutate(eating_rate == rate)
    return(food_mice_cats)
}

Answer

final_state <- function(rate) {
   result2 <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
                reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0)
   
   answer <- c( growth = rate, food = result2$food[180], total_cats = sum(result2$cats[180], result2$fat_cats[180]), total_mice = sum(result2$mice[180], result2$fat_mice[180])
   )
   answer
}

Answer

final_state <- function(rate) {
  var1 <-tail(
    food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3, reprod_rate=0.04,
                catch_rate=0.05, death_rate=0.01, eating_rate=rate,
                hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
                mice_ini=10, fat_mice_ini=0, cats_ini=10, fat_cats_ini=0
    ), 1)
  food <- var1$food
  total_cats <- var1$cats + var1$fat_cats
  total_mice <- var1$mice + var1$fat_mice

  return (as.double(data.frame(rate, food, total_cats, total_mice)))
}

Answer

final_state <- function(rate) {
    temporary_result <- food_mice_cats(N=180, birth_rate=0.1, growth_rate=1e-3,
          reprod_rate=0.04, catch_rate=0.05, death_rate=0.01, eating_rate=rate,
          hunger_rate=1e-3, dying_rate=0.05, food_ini=100,
          mice_ini=10, fat_mice_ini=0, cats_ini=10,fat_cats_ini=0)
  growth<-rate
  food<-temporary_result$food[180]
  total_cats<-temporary_result$cats[180] + temporary_result$fat_cats[180]
  total_mice<-temporary_result$mice[180] + temporary_result$fat_mice[180]
  names(growth)<-"growth"
  names(food)<-"food"
  names(total_mice)<-"total_mice"
  names(total_cats)<-"total_cats"
  return(c(growth,food,total_cats,total_mice))
}