library(package = "deSolve")
library(package = "viridis")
## Loading required package: viridisLite

1. Comparing and entering the numerical and analytic solutions to differential equations

Now that you have started to use R’s numerical solver, I want you to be able to see its solutions compared to the analytic one that is found by integration. To this end, and to give you some practice using ode() and plotting functions in R, please do the following exercises.

1.1. Let’s first start with density-independent growth. Please plot 4 time series in the same plot: a numerical solution for \(r>0\), an analytically solution for \(r>0\), a numerical solution for \(r<0\), and an analytically solution for \(r<0\). First, plot the time series for density-independent growth using ode(). Set \(N(0) = 10\) and \(r = 0.2\). Then, use the analytic solution to the density-independent differential equation to generate a time series and add it to the plot using lines() (hint: see problem 2 from the last problem set to plot a vector of an independent variable against a dependent one.) In lines(), add , col = "red", lty = 2, lwd = 5 to make the line dotted and thick (lwd is the argument for line width). Next, repeat what you just did in the previous 2 steps, but with \(r = -0.2\). Add them to the same graph and make sure that your vertical axis shows all of the curves completely.

1.2. Do the same thing as above, but with the logistic equation. Fist, plot the logistic equation’s numerical solution (black, solid line) and analytic solution (red, dotted line). The logistic Do this with \(N(0) = 10\), \(r = 0.2\), and \(K = 75\).

1.1.

# Numerical solution, r > 0
densInd <- function(time, state, parms) {
  with(as.list(c(state, parms)), {
    dN <- r*N
    return(list(c(dN)))
  })
}

N0 <- c(N = 10)
time_vec <- seq(from = 0, to = 10, length.out = 1e3)
parm_vec1 <- c(r = 0.2)

densIndNS1 <- ode(y = N0, func = densInd, parms = parm_vec1, times = time_vec)

plot(x = densIndNS1[,1], y = densIndNS1[,2], ylim = c(0, 75), las = 1, xlab = "Time", ylab = "Density", type = "l")

# Analytic solution, r > 0
Nt <- N0*exp(parm_vec1["r"]*time_vec)
lines(x = time_vec, y = Nt, col = "red", lty = 3, lwd = 5)

# Numerical solution, r < 0
parm_vec2 <- c(r = -0.2)

densIndNS2 <- ode(y = N0, func = densInd, parms = parm_vec2, times = time_vec)

lines(x = densIndNS2[,1], y = densIndNS2[,2])

# Analytic solution, r < 0
Nt <- N0*exp(parm_vec2["r"]*time_vec)
lines(x = time_vec, y = Nt, col = "blue", lty = 3, lwd = 5)

1.2.

# Numerical solution, r > 0
densDep <- function(time, state, parms) {
  with(as.list(c(state, parms)), {
    dN <- r*N*(K - N)/K
    return(list(c(dN)))
  })
}

N0 <- c(N = 10)
time_vec <- seq(from = 0, to = 35, length.out = 1e3)
parm_vec1 <- c(r = 0.2, K = 75)

densDepNS1 <- ode(y = N0, func = densDep, parms = parm_vec1, times = time_vec)

plot(x = densDepNS1[,1], y = densDepNS1[,2], ylim = c(0, 80), las = 1, xlab = "Time", ylab = "Density", type = "l")

# Analytic solution, r > 0
Nt <- parm_vec1["K"]*N0/((parm_vec1["K"] - N0)*exp(-parm_vec1["r"]*time_vec) + N0)
lines(x = time_vec, y = Nt, col = "red", lty = 3, lwd = 5)

# Numerical solution, r < 0
parm_vec2 <- parm_vec1
parm_vec2["r"] <- -0.2

densDepNS2 <- ode(y = N0, func = densDep, parms = parm_vec2, times = time_vec)

lines(x = densDepNS2[,1], y = densDepNS2[,2])

# Analytic solution, r < 0
Nt <- parm_vec2["K"]*N0/((parm_vec2["K"] - N0)*exp(-parm_vec2["r"]*time_vec) + N0)
lines(x = time_vec, y = Nt, col = "blue", lty = 3, lwd = 5)

2. Parameter analysis

In this course we are studying general population models. You have seen various models used to represent different types of biological scenarios, like populations that change in a density-independent manner, populations that change in a density-dependent manner, populations that change in steps, populations that change continuously, etc. But once we have a model, what then? Well, then we study the model.

In Botsford et al. (the chapter that we read in week 2), they argued that understanding and prediction are the two most important things that we we get out of models. Prediction is important, but beyond the scope of this course, unless this course was years long and we would build models, evaluate their predictive power, modify the model, evaluate, etc. But to understand population ecology through models, this entails analyzing them, which can really mean anything.

So far in this course you have done some flavors of analysis. First, you’ve done graphical analysis of our models. You have taken the dynamic models and their solutions and graphed them. Time series are a form of analysis, which many of you presumably saw since many of your potential remodeling projects had time series figures. Second, you’ve taken your dynamic models and plotted them as functions of density. That is a form of analysis, as important population biology is revealed (e.g., \(r\) is the maximum per-capita growth rate, population growth rate and per-capita growth rate both intercept the horizontal axis as the same point [\(K\), \(r/\alpha\), etc.], etc.). Third, you’ve done equilibrium analysis by finding \(N\) for which the rate of change is 0; i.e., where the population is not changing. Fourth and last, you’ve taken dynamic models and their solutions, and you have varied parameter values to see how that changes the model behavior/output. In the last problem set, for instance, you took the logistic model from Stevens text and found that by modifying \(r_d\) that the population may approach a single value, that it may oscillate around a value before settling on it, that it may oscillate around the value forever, or it may change chaotically/aperiodically. Note you will be learning another type of analysis this week where we evaluate a model’s behavior around equilibria, which is what we call stability analysis. So, pat yourselves on the back for analyzing models without being explicit about it.

For this problem, I will ask you to a more extensive and more formal analysis of model parameters. This is, essentially, the most significant type of analysis that you can do to learn from most population models because mathematical analysis is often limited and not useful. An example of a mathematical analysis that is not useful might be when you find an equilibrium value for a model, but it has a dozen+ terms that are difficult to make sense of mathematically or biologically.

2.1. Let’s first pick a model to analyse: we’ll use the \(r\)-\(\alpha\) logistic equation: \[ \frac{\mathrm{d}N}{\mathrm{d}t} = N(r - \alpha N).\] Imagine that this is your first time seeing this model and you don’t know anything about it. If you were trying to understand the model, an important and primary question might be: how do the parameters, \(r\) and \(\alpha\) affect the model behavior? Well, let’s see. First, take any 3 values of \(r\) and plot the three numerical solutions on the same time series. Use whatever \(N_0\), \(\alpha\), and length of time you desire/is useful for observing the model’s behavior. Second, make a plot with the 3 values of \(r\) on the horizontal axis and the final density from the 3 numerical solutions on the vertical axis. As suggestion for a better practice, try not just copying and pasting the final value, but index the last value by taking the second column of the last row. So for instance, if you saved your numerical solution as ralpha_1, I would index the value like this: ralpha_1[nrow(ralpha_1), 2]. Last, comment in one or more sentences on the effect of \(r\) on density.

2.2. Now, do the same thing with \(\alpha\). First, take any 3 values of \(\alpha\) and plot the three numerical solutions on the same time series. Use whatever \(N_0\), \(r\), and length of time you desire/is useful for observing the model’s behavior. Second, make a plot with the 3 values of \(\alpha\) on the horizontal axis and the final density from the 3 numerical solutions on the vertical axis. Last, comment in one or more sentences on the effect of \(\alpha\) on density.

2.3. You’re going to be annoyed with me for this, for which I apologize in advance 😬. For the 3 values of \(r\) and 3 values of \(\alpha\) above, I want you to first plot all 9 combinations. Second, I want you to save all nine final density values in a matrix. To do this, first, make a matrix. I’ll do this for you: ralpha_mat <- matrix(data = NA, nrow = 3, ncol = 3). Now, just for the sake of helping you know where to save the final densities, let’s label the rows and columns like this: rownames(ralpha_mat) <- paste0("r", 1:3) and colnames(ralpha_mat) <- paste0("alpha", 1:3). Now, when you look at ralpha_mat, you see:

##    alpha1 alpha2 alpha3
## r1     NA     NA     NA
## r2     NA     NA     NA
## r3     NA     NA     NA

You know how to index rows and columns, but as a reminder, if I had a matrix called le_matrix, I could save a value in its third row and fourth column like le_matrix[3, 4]. Third and last, once you have the matrix filled, I want you to plot the final densities of \(r\) against \(\alpha\) using a heat map. Do this by using the function image() that comes in R’s base graphics. If you just throw in your matrix, you’ll notice the axes are flipped: that’s annoying but I have no control over this. So, a necessary step you’ll need to take is to first rotate the matrix using (called a transposition in matrix-speak): le_matrix1 <- t(le_matrix1). Then, inside of image(), give the following arguments:

  1. y =: put a vector of your \(\alpha\) values in; e.g., image(y = c(0.1, 1, 10))
  2. x =: put a vector of your \(r\) values in
  3. z =: this is where you enter your matrix
  4. col = heat.colors(n = 100): this sets the color palette such that large values are white and small values are red; I can show you how to make your own color palettes easily, later, if you’d like
  5. las =: because we always do this for the sake of people’s necks

After you make your beautiful heat map, take a moment to revel and reflect. Write a sentence or few on it.

Code used across question 2

# r-alpha function
  ralpha <- function(time, state, pars) {
    with(as.list(c(state, pars)), {
      dN <- N*(r - alpha*N)
      return(list(c(dN)))
    })
  }

# Time sequence and N0 for ode()
  time_vec <- seq(from = 0, to = 100, length.out = 1000)
  N0 <- c(N = 0.05)

2.1.

# Time series

  r_vec <- c(0.1, 0.2, 0.4)
  parm_vec1 <- c(alpha = 0.01, r = r_vec[1])
  parm_vec2 <- c(alpha = 0.01, r = r_vec[2])
  parm_vec3 <- c(alpha = 0.01, r = r_vec[3])
  
  ode_out1 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec1)
  ode_out2 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec2)
  ode_out3 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec3)
  
  plot(x = ode_out1[,1], ode_out1[,2], xlab = "Time", ylab = "Density", las = 1, type = "l", lwd = 2, ylim = c(0, 40))
    lines(x = ode_out2[,1], ode_out2[,2], col = "tomato", lwd = 2)
    lines(x = ode_out3[,1], ode_out3[,2], col = "skyblue", lwd = 2)

# r v. final densities
  # Find the number of rows from your output b/c we want the last one (all 3 will be the same)
  last_time <- nrow(ode_out1)
  
  final_dens <- c(ode_out1[last_time, 2], ode_out2[last_time, 2], ode_out3[last_time, 2])

  plot(x = r_vec, y = final_dens, pch = 19, cex = 1.25, las = 1, xlab = "Intrinsic rate of growth (r)", ylab = "Final density (N)", col = "seagreen3")
    lines(x = r_vec, y = final_dens, lwd = 2, col = "seagreen3")

2.2.

# Time series

  alpha_vec <- c(0.1, 0.3, 0.9)
  parm_vec1 <- c(alpha = alpha_vec[1], r = 1)
  parm_vec2 <- c(alpha = alpha_vec[2], r = 1)
  parm_vec3 <- c(alpha = alpha_vec[3], r = 1)
  
  ode_out1 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec1)
  ode_out2 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec2)
  ode_out3 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec3)
  
  plot(x = ode_out1[,1], ode_out1[,2], xlab = "Time", ylab = "Density", las = 1, type = "l", lwd = 2, ylim = c(0, 10))
    lines(x = ode_out2[,1], ode_out2[,2], col = "tomato", lwd = 2)
    lines(x = ode_out3[,1], ode_out3[,2], col = "skyblue", lwd = 2)

# alpha v. final densities
  # Find the number of rows from your output b/c we want the last one (all 3 will be the same)
  last_time <- nrow(ode_out1)
  
  final_dens <- c(ode_out1[last_time, 2], ode_out2[last_time, 2], ode_out3[last_time, 2])

  plot(x = alpha_vec, y = final_dens, pch = 19, cex = 1.25, las = 1, xlab = "Crowding coefficient (alpha)", ylab = "Final density (N)", col = "steelblue3")
    lines(x = alpha_vec, y = final_dens, lwd = 2, col = "steelblue3")

2.3.

# Time series

  alpha_vec <- c(0.1, 0.3, 0.9)
  parm_vec1 <- c(alpha = alpha_vec[1], r = r_vec[1])
  parm_vec2 <- c(alpha = alpha_vec[2], r = r_vec[1])
  parm_vec3 <- c(alpha = alpha_vec[3], r = r_vec[1])
  parm_vec4 <- c(alpha = alpha_vec[1], r = r_vec[2])
  parm_vec5 <- c(alpha = alpha_vec[2], r = r_vec[2])
  parm_vec6 <- c(alpha = alpha_vec[3], r = r_vec[2])
  parm_vec7 <- c(alpha = alpha_vec[1], r = r_vec[3])
  parm_vec8 <- c(alpha = alpha_vec[2], r = r_vec[3])
  parm_vec9 <- c(alpha = alpha_vec[3], r = r_vec[3])

  ode_out1 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec1)
  ode_out2 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec2)
  ode_out3 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec3)
  ode_out4 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec4)
  ode_out5 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec5)
  ode_out6 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec6)
  ode_out7 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec7)
  ode_out8 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec8)
  ode_out9 <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec9)
  
  col_vec <- rainbow(n = 9)
  plot(x = ode_out1[,1], ode_out1[,2], xlab = "Time", ylab = "Density", las = 1, type = "l", lwd = 2, ylim = c(0, 5), col = col_vec[1])
    lines(x = ode_out2[,1], ode_out2[,2], col = col_vec[1], lwd = 2)
    lines(x = ode_out3[,1], ode_out3[,2], col = col_vec[2], lwd = 2)
    lines(x = ode_out4[,1], ode_out4[,2], col = col_vec[3], lwd = 2)
    lines(x = ode_out5[,1], ode_out5[,2], col = col_vec[4], lwd = 2)
    lines(x = ode_out6[,1], ode_out6[,2], col = col_vec[5], lwd = 2)
    lines(x = ode_out7[,1], ode_out7[,2], col = col_vec[6], lwd = 2)
    lines(x = ode_out8[,1], ode_out8[,2], col = col_vec[7], lwd = 2)
    lines(x = ode_out9[,1], ode_out9[,2], col = col_vec[8], lwd = 2)

# alpha v. final densities
  # Find the number of rows from your output b/c we want the last one (all 3 will be the same)
  last_time <- nrow(ode_out1)
  
  # Create matrix
  ralpha_mat <- matrix(data = NA, nrow = 3, ncol = 3)
  rownames(ralpha_mat) <- paste0("r", 1:3)
  colnames(ralpha_mat) <- paste0("alpha", 1:3)
  
  # Save values in matrix
  ralpha_mat[1, 1] <- ode_out1[last_time, 2]
  ralpha_mat[1, 2] <- ode_out2[last_time, 2]
  ralpha_mat[1, 3] <- ode_out3[last_time, 2]
  ralpha_mat[2, 1] <- ode_out4[last_time, 2]
  ralpha_mat[2, 2] <- ode_out5[last_time, 2]
  ralpha_mat[2, 3] <- ode_out6[last_time, 2]
  ralpha_mat[3, 1] <- ode_out7[last_time, 2]
  ralpha_mat[3, 2] <- ode_out8[last_time, 2]
  ralpha_mat[3, 3] <- ode_out9[last_time, 2]

  # Heat map
  ralpha_mat_t <- t(ralpha_mat)
  image(y = r_vec, x = alpha_vec, z = ralpha_mat, xlab = "Intrinsic rate of growth (r)", ylab = "Crowding coefficiet (alpha)", col = heat.colors(n = 100), las = 1)

3. Harnessing the power of computational iteration

So, hopefully you’re seeing the point of what we’re trying to do with parameter analysis. But, you might be thinking: Chris, this is amazing, but WTH am I going to do with more than a handful a parameters? It seems like it’d be a huge pain in the neck! And with all of your preaching about las = 1, you are against pains in the neck! It is true, I am against pains in the neck. So, let’s scale this 💩 up! (Yes, I’ve been drinking a lot of ☕: how can you tell?)

Our computers can be our friends; especially when we harness their power to do work for us! We’re going to do this below using a type of iterative iterative process called a loop. Loops are your friends: you give them a condition, like Hey loop, go fetch me the final densities of a million values of \(r\) using my numerical solver., and they just do it. They’re really good friends.

Loops come in a variety of forms, but we are essentially just going to be using what we call for loops to repeat a calculations over a set of values, like in a vector. You basically tell R to perform some set of actions and over how many times. A simple loop to, say, calculate the area of a rectangle if we have a length of 2 and 5 lengths would be:

# One length and 5 width values
length_rec <- 2
width_rec <- 2:6

# I have the number of iterations I want to do, which is the length of the width vector
n_widths <- length(width_rec)

# Create a data structure to save the area values
area_vec <- rep(x = NA, times = n_widths)

# The for loop
for (i in 1:n_widths) {
  width_i <- width_rec[i]
  area_i <- length_rec*width_i
  area_vec[i] <- area_i
}

# Let's plot and look at what we just did
# (Note I make the axes 1:1 to better observe how a 1 unit change in width corresponds to change in area)
plot(x = width_rec, y = area_vec, las = 1, xlab = "Width", ylab = "Area", pch = 19, xlim = c(2, 12), ylim = c(2, 12))
  lines(x = width_rec, y = area_vec)

Cool, huh! I’ll first say that for this very simple problem, the simplest way to actually do this is just length_rec*width_rec: that will give you the same result as the for loop. But, this example was supposed to be illustrative or what you could potentially do. Second, study the code. Several things should pop out:

3.1. Hopefully when reading the example loop above you were thinking about how you might do this with the parameter values from question 2. Well, that’s what I want you to do. For this problem, do the same thing that you did in 2.1, but instead of choosing 3 arbitrary \(r\) values, create a vector of 100 values. Start and end at some arbitrary value of \(r\) and create a sequence of 100. If I arbitrarily started at -1 and ended at 10, I would create my sequence as: r_seq <- seq(from = -1, to = 10, length.out = 100).

3.2. Similar to the previous problem, recreate 2.2, but with \(\alpha\). Choose arbitrary values to start and stop, and analyze 100 values. Please use \(\alpha >0\), because, unlike \(r\), switching the sign of \(\alpha\) to be \(>0\) doesn’t make biological sense for the logistic model.

3.3. (Optional and recommended if you can make the time) I imagine that it took a bit of effort to complete 3.1. and 3.2. Pat yourself on the back for a job well done: you deserve it. For this final problem, as you guessed, I want you to do the same thing as 2.3. To do this, we need to do a just few things differently. First, you need to write a loop within a loop—this is what we call a nested loop. It looks like this:

for (i in 1:10) {
  # Here is where i is being indexed
  for (j in 1:10) {
    # Here's where j is being indexed
    # Here's where you put your calculations
  }
}

You can see that the first difference compared with above is that there is now a loop nested in another. Another thing that we need to do differently is to give the nested loop variable a different value; here, we use j. What happens here is that i takes the value 1, then j is looped over to take 1, 2, … 10. Then, i takes the value 2, then j is looped over to take 1, 2, … 10. This happens until i takes on the value 10, then j is looped over to take 1, 2, … 10. The last thing that we need to do differently is create a data structure that can take all the combinations of i and j: you probably guessed it, a matrix.

Here is an example again using nested loops to find the area of a rectangle, where I loop over length and width:

# I create a named value for the number of length and width values I want to use
n_vals <- 100

# Create vectors for the lengths and widths
length_vec <- seq(from = 0, to = 20, length.out = n_vals)
width_vec <- seq(from = 0, to = 20, length.out = n_vals)

# Create a matrix to save the values in
area_mat <- matrix(data = NA, nrow = n_vals, ncol = n_vals)
  rownames(area_mat) <- paste0("l", 1:n_vals)
  colnames(area_mat) <- paste0("w", 1:n_vals)
  # Uncomment to check out the first 5 rows and columns and see that we want lengths in rows, widths in columns
  # area_mat[1:5, 1:5]

# The loop
for (i in 1:n_vals) {
  # Here is where length i is being indexed
  l_i <- length_vec[i]
  for (j in 1:n_vals) {
    # Here's where width j is being indexed
    w_j <- length_vec[j]
    # Here's where you put your calculations
    area <- l_i*w_j
    # Here is where we save the ith row for length and jth column for width
    area_mat[i, j] <- area
  }
}

# Transpose the matricx and plot
area_mat1 <- t(area_mat)
image(x = width_vec, y = length_vec, z = area_mat1, las = 1, col = heat.colors(n = 100), xlab = "Width", ylab = "Length")

# Note there are different functions that allow you to plot 3D data like these, including contour() and filled.contour(). Uncomment to see what they look like:
# contour(x = width_vec, y = length_vec, z = area_mat1, las = 1, xlab = "Width", ylab = "Length")
#filled.contour(x = width_vec, y = length_vec, z = area_mat1, las = 1, xlab = "Width", ylab = "Length")

Isn’t that just sooooo coooooool 😎. So, do this with the \(r\)-\(\alpha\) logistic equation. I’m eager to see your heat maps!

Code used across question 3

# r-alpha function
  ralpha <- function(time, state, pars) {
    with(as.list(c(state, pars)), {
      dN <- N*(r - alpha*N)
      return(list(c(dN)))
    })
  }

# Time sequence and N0 for ode()
  time_vec <- seq(from = 0, to = 100, length.out = 1000)
  N0 <- c(N = 0.05)

# Create a vector for parameters
  parm_vec <- c(alpha = 0.01, r = 0.01)

3.1.

# Time series
 length_rs <- 100
 r_vec <- seq(from = 0.1, to = 5.1, length.out = length_rs)
 
 length_times <- length(time_vec)
 
 rdens_vec <- rep(x = NA, times = length_rs)

 # For loop
 for (i in 1:length_rs) {
   # Save the ith r into the parameter vector
   r_i <- r_vec[i]
   parm_vec["r"] <- r_i
 
   # Run ode()
   ode_outi <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec)
   
   # Save ode results
   rdens_vec[i] <- ode_outi[length_times, 2]
 }

# r v. final densities
 # Find the number of rows from your output b/c we want the last one (all 3 will be the same)

 plot(x = r_vec, y = rdens_vec, pch = 19, cex = 1, las = 1, xlab = "Intrinsic rate of growth (r)", ylab = "Final density (N)", col = "maroon2")
   lines(x = r_vec, y = rdens_vec, lwd = 2, col = "maroon2")

3.2

# Time series
  length_alphas <- 100
  alpha_vec <- seq(from = 0.1, to = 0.6, length.out = length_alphas)
  
  length_times <- length(time_vec)
  
  alphadens_vec <- rep(x = NA, times = length_alphas)

  # For loop
  for (i in 1:length_alphas) {
    # Save the ith r into the parameter vector
    alpha_i <- alpha_vec[i]
    parm_vec["alpha"] <- alpha_i
  
    # Run ode()
    ode_outi <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec)
    
    # Save ode results
    alphadens_vec[i] <- ode_outi[length_times, 2]
  }

# r v. final densities
  # Find the number of rows from your output b/c we want the last one (all 3 will be the same)

  plot(x = r_vec, y = alphadens_vec, pch = 19, cex = 1, las = 1, xlab = "Crowding coefficient (alpha)", ylab = "Final density (N)", col = "slateblue3")
    lines(x = r_vec, y = alphadens_vec, lwd = 2, col = "slateblue3")

3.3.

  ralpha_mat <- matrix(data = NA, nrow = length_rs, ncol = length_alphas)
    rownames(ralpha_mat) <- paste0("r", 1:length_rs)
    colnames(ralpha_mat) <- paste0("c", 1:length_alphas)
  
  for (i in 1:length_rs) {
    r_i <- r_vec[i]
    parm_vec["r"] <- r_i
    for (j in 1:length_alphas) {
      alpha_j <- alpha_vec[j]
      parm_vec["alpha"] <- alpha_j
    
      ode_out <- ode(y = N0, times = time_vec, func = ralpha, parms = parm_vec)
      ralpha_mat[i, j] <- ode_out[length_times, 2]
    }
  }

  image(x = r_vec, y = alpha_vec, z = ralpha_mat, las = 1, xlab = "Intrinsic rate of grwoth (r)", ylab = "Crowding coefficient (alpha)", col = viridis(n = 100))