Set up the global options

###Part 1 Tracking Three Populations (2 Points)

Set up a data frame to track three populations. We wil have a food chain in which we have mice, snakes, and hawks. Mice are our prey species. Snakes will feed on mice. Hawks will feed on both snakes and mice. Start by setting up a data frame that tracks the three populations over 100 time periods. All three species will follow exponential growth.

Starting Population(N0) for Mice is 100, their growth rate is 15% Starting Population(N0) for Snakes is 20, their growth rate is 5% Starting Population(N0) for Hawks is 5, their growth rate is 2%

The code below will need to be edited. Also be sure to add comments to explain what is being done in each step, especially for code that is provided.

mice<-c()
snakes<-c()
hawks<-c()

Mouse_N0<-100
Snake_N0<-20
Hawk_N0<-5
  
Mouse_r<-0.15
Snake_r<-0.05
Hawk_r<-0.02
  
#Use a for loop to track population. Below is an EXAMPLE and will need to be edited

  for (i in 2:100){
    mice[1]<-Mouse_N0
    snakes[1]<-Snake_N0
    hawks[1]<-Hawk_N0
    mice[i]<-Mouse_r*mice[i-1]
    snakes[i]<-Snake_r*snakes[i-1]
    hawks[i]<-Hawk_r*hawks[i-1]
  }


Tracking_Populations<-data.frame("Mouse_Population"=mice, "Snake_Population"=snakes, "Hawk_Population"=hawks, "Time" =1:100)

#Create a plot that plots the mouse, snake, and hawk populations on the same figure

#Hint start with a plot, then add points() and lines

plot(Tracking_Populations$Time, Tracking_Populations$Mouse_Population, type="l", col="blue", ylim=c(0, max(c(mice, snakes, hawks), na.rm=TRUE)), xlab="Time", ylab="Population Size", main="Population Dynamics")
lines(Tracking_Populations$Time, Tracking_Populations$Snake_Population, col="green")
lines(Tracking_Populations$Time, Tracking_Populations$Hawk_Population, col="red")
legend("topright", legend=c("Mice", "Snakes", "Hawks"), col=c("blue", "green", "red"), lty=1)

###Part 2 Adding Stochasticity (2 Points)

In our example, mice populations are highly variable but snakes and hawks have more stable populations. Paste your code from above into the chunk below. Edit your code to make the following changes

mice have growth rates that vary from -1% to +9% snakes have growth rates that vary from 3% to 6% hawks have growth rates that vary from 1.8% to 2.5%

#This example generates a variable growth rate between -5% and 5%. Hint, you will need to move your growth rate from outside of your for loop to inside of it.
mice<-c()
snakes<-c()
hawks<-c()

Mouse_N0<-100
Snake_N0<-20
Hawk_N0<-5
  
Mouse_r<-0.15
Snake_r<-0.05
Hawk_r<-0.02
  
#Use a for loop to track population. Below is an EXAMPLE and will need to be edited

  for (i in 2:100){
    mice[1]<-Mouse_N0
    snakes[1]<-Snake_N0
    hawks[1]<-Hawk_N0
    Mouse_r <- runif(1, -0.01, 0.09)  # Mice growth rate between -1% and +9%
    Snake_r <- runif(1, 0.03, 0.06)  # Snakes growth rate between 3% and 6%
    Hawk_r <- runif(1, 0.018, 0.025)  # Hawks growth rate between 1.8% and 2.5%
    mice[i]<-Mouse_r*mice[i-1]
    snakes[i]<-Snake_r*snakes[i-1]
    hawks[i]<-Hawk_r*hawks[i-1]
  }


Tracking_Populations<-data.frame("Mouse_Population"=mice, "Snake_Population"=snakes, "Hawk_Population"=hawks, "Time" =1:100)

# Create the plot
plot(Tracking_Populations$Time, Tracking_Populations$Mouse_Population, type="l", col="blue", 
     ylim=c(0, max(c(mice, snakes, hawks), na.rm=TRUE)), 
     xlab="Time", ylab="Population Size", main="Population Dynamics")

# Add lines for snakes and hawks
lines(Tracking_Populations$Time, Tracking_Populations$Snake_Population, col="green")
lines(Tracking_Populations$Time, Tracking_Populations$Hawk_Population, col="red")

# Add a legend
legend("topright", legend=c("Mice", "Snakes", "Hawks"), col=c("blue", "green", "red"), lty=1)

###Part 3 Adding a 5 year storm (2 points)

In this example, the area these populations live is storm prone. ON AVERAGE, a major storm occurs every five years. This does not mean a storm happens every fifth year

Mice and snakes are both especially vulenerable to storms since they live in burrows but hawks are less effected.

If a storm occurs, 25% of mice will die, 20% of snakes will die, and 5% of hawks will die.

Add an if statement that incorperates the storm dynamics. Hint, you will need one if statement that captures the storm dynamic that then influences the three populations

#There is a 33% chance a storm occurs in this example. This for loop will remove 50% of the example population if a storm occurs. This code will need to be modified and incorperated into your for loop which tracks populaiton


#####PLEASE NOTE, in the video, the code says Example_Population[i]<-0.5*Example_Population[i-1], it should sayExample_Population[i-1] as it does below

mice<-c()
snakes<-c()
hawks<-c()

Mouse_N0<-100
Snake_N0<-20
Hawk_N0<-5
  
Mouse_r<-0.15
Snake_r<-0.05
Hawk_r<-0.02
  
#Use a for loop to track population. Below is an EXAMPLE and will need to be edited

  for (i in 2:100){
  mice[1]<-Mouse_N0
  snakes[1]<-Snake_N0
  hawks[1]<-Hawk_N0
  
  # Determine if a storm occurs (on average every 5 years, but randomly)
  Storm_Chance <- runif(1, 0, 1)
  if (Storm_Chance < 0.20) {  # Adjusted probability to match storm frequency
    mice[i-1] <- mice[i-1] * 0.75  # 25% mortality for mice
    snakes[i-1] <- snakes[i-1] * 0.80  # 20% mortality for snakes
    hawks[i-1] <- hawks[i-1] * 0.95  # 5% mortality for hawks
  }
  
  # Randomized growth rates
  Mouse_r <- runif(1, -0.01, 0.09)
  Snake_r <- runif(1, 0.03, 0.06)
  Hawk_r <- runif(1, 0.018, 0.025)
  
  # Population growth equations
  mice[i] <- Mouse_r * mice[i-1]
  snakes[i] <- Snake_r * snakes[i-1]
  hawks[i] <- Hawk_r * hawks[i-1]
}

Tracking_Populations<-data.frame("Mouse_Population"=mice, "Snake_Population"=snakes, "Hawk_Population"=hawks, "Time" =1:100)

# Create a plot that plots the mouse, snake, and hawk populations on the same figure

plot(Tracking_Populations$Time, Tracking_Populations$Mouse_Population, type="l", col="blue", 
     ylim=c(0, max(c(mice, snakes, hawks), na.rm=TRUE)), 
     xlab="Time", ylab="Population Size", main="Population Dynamics")
lines(Tracking_Populations$Time, Tracking_Populations$Snake_Population, col="green")
lines(Tracking_Populations$Time, Tracking_Populations$Hawk_Population, col="red")
legend("topright", legend=c("Mice", "Snakes", "Hawks"), col=c("blue", "green", "red"), lty=1)

###Part 4 Linking the Food Chain (6 points) Until this point, all three of these populations have been unrelated (i.e. there is no term which connects there growth rates). In this example, snakes eat mice, and hawks eat both mice and snakes. The only way mice die is due to snakes and hawks. The only way snakes die is due to being eaten by hawks. Hawks have a natural death rate.

Using the predator-prey models from the previous section as an example, modifiy your code as follows.

The mice population is equal to:

\[ N_{mice_{t}} = r_{mice}*N_{mice{t-1}}- \beta_1*N_{mice{t-1}}* N_{snake_{t-1}}-\beta_2*N_{mice{t-1}}*N_{hawk_{t-1}} \]

Where beta1 is the encounter rate for mice and snakes and beta 2 is the encounter rate for mice and hawks

The snake population is equal to:

\[ N_{snake_{t}} = \delta_1*\beta_1*N_{mice{t-1}}*N_{snake_{t-1}}-\beta_3*N_{snakes{t-1}}*N_{hawk_{t-1}} \]

Where gamma1 is equalt to the how effectively snakes turn mice into food and beta3 is how often hawks encounter and eat snakes.

Finally, the hawk population is equal to:

\[ N_{hawk_{t}} = \delta_2*\beta_3*N_{snakes{t-1}}*N_{hawks_{t-1}}+\delta_3*\beta_2*N_{mice{t-1}}*N_{hawk_{t-1}}-\gamma*N_{hawk_{t-1}} \]

Where delta 2 and 3 are the conversation rate for hawks turning snakes and mice into energy / offspring respectively and gamma is the death rate of hawks

Use the following rates: Mice growth rate should stay the same as you have it set beta1<-0.03 beta2<-0.01 beta3<-0.02 delta1<-0.25 delta2<-0.3 delta3<-0.15 gamma<-0.3

#Change your population growth terms to include the interactions. You can copy and paste code from the example script from the previous competition code but it will need to be modified. Also note, you will only keep the random growth rate for mice, snakes and hawks growth rate will now just be due to how they catch and eat prey. 
mice <- c()
snakes <- c()
hawks <- c()

Mouse_N0 <- 100
Snake_N0 <- 20
Hawk_N0 <- 5

# Growth rate remains the same for mice
Mouse_r <- 0.15

# Predator-prey interaction parameters
beta1 <- 0.03  # Encounter rate between mice and snakes
beta2 <- 0.01  # Encounter rate between mice and hawks
beta3 <- 0.02  # Encounter rate between snakes and hawks
delta1 <- 0.25 # Conversion efficiency of snakes consuming mice
delta2 <- 0.3  # Conversion efficiency of hawks consuming snakes
delta3 <- 0.15 # Conversion efficiency of hawks consuming mice
gamma <- 0.3   # Natural death rate of hawks


for (i in 2:100) {
  mice[1] <- Mouse_N0
snakes[1] <- Snake_N0
hawks[1] <- Hawk_N0
  # Determine if a storm occurs (on average every 5 years, but randomly)
  Storm_Chance <- runif(1, 0, 1)
  if (Storm_Chance < 0.20) {  # Adjusted probability to match storm frequency
    mice[i-1] <- mice[i-1] * 0.75  # 25% mortality for mice
    snakes[i-1] <- snakes[i-1] * 0.80  # 20% mortality for snakes
    hawks[i-1] <- hawks[i-1] * 0.95  # 5% mortality for hawks
  }
  # Stochastic mouse growth rate
  Mouse_r <- runif(1, -0.01, 0.09)
  
  # Population equations incorporating predation
  mice[i] <- mice[i-1] + (Mouse_r * mice[i-1]) - (beta1 * mice[i-1] * snakes[i-1]) - (beta2 * mice[i-1] * hawks[i-1])
  snakes[i] <- (delta1 * beta1 * mice[i-1] * snakes[i-1]) - (beta3 * snakes[i-1] * hawks[i-1])
  hawks[i] <- (delta2 * beta3 * snakes[i-1] * hawks[i-1]) + (delta3 * beta2 * mice[i-1] * hawks[i-1]) - (gamma * hawks[i-1])
  
  # Ensure populations don't go negative
  mice[i] <- max(mice[i], 0)
  snakes[i] <- max(snakes[i], 0)
  hawks[i] <- max(hawks[i], 0)
}

Tracking_Populations <- data.frame("Mouse_Population"=mice, "Snake_Population"=snakes, "Hawk_Population"=hawks, "Time"=1:100)

# Plot population trends
plot(Tracking_Populations$Time, Tracking_Populations$Mouse_Population, type="l", col="blue", 
     ylim=c(0, max(c(mice, snakes, hawks), na.rm=TRUE)), 
     xlab="Time", ylab="Population Size", main="Predator-Prey Population Dynamics")
lines(Tracking_Populations$Time, Tracking_Populations$Snake_Population, col="green")
lines(Tracking_Populations$Time, Tracking_Populations$Hawk_Population, col="red")
legend("topright", legend=c("Mice", "Snakes", "Hawks"), col=c("blue", "green", "red"), lty=1)

###Part 5

If any population crashes to 0, that population can no longer grow. Add an if statement for the populations stops future population growth if a population goes to 0.

#You should have three if loops, one before each population equation
# Initialize populations
mice <- c()
snakes <- c()
hawks <- c()

# Starting populations
Mouse_N0 <- 100
Snake_N0 <- 20
Hawk_N0 <- 5

# Growth rates and parameters
Mouse_r <- 0.15
beta1 <- 0.03  # Encounter rate for mice and snakes
beta2 <- 0.01  # Encounter rate for mice and hawks
beta3 <- 0.02  # Encounter rate for snakes and hawks
delta1 <- 0.25 # Conversion rate of mice to snakes
delta2 <- 0.3  # Conversion rate of snakes to hawks
delta3 <- 0.15 # Conversion rate of mice to hawks
gamma <- 0.3   # Natural death rate for hawks

# Simulation loop
for (i in 2:100) {
  # Initial populations
  mice[1] <- Mouse_N0
  snakes[1] <- Snake_N0
  hawks[1] <- Hawk_N0
  # Determine if a storm occurs (on average every 5 years, but randomly)
  Storm_Chance <- runif(1, 0, 1)
  if (Storm_Chance < 0.20) {  # Adjusted probability to match storm frequency
    mice[i-1] <- mice[i-1] * 0.75  # 25% mortality for mice
    snakes[i-1] <- snakes[i-1] * 0.80  # 20% mortality for snakes
    hawks[i-1] <- hawks[i-1] * 0.95  # 5% mortality for hawks
  }
  # Stochastic mouse growth rate
  Mouse_r <- runif(1, -0.01, 0.09)
  # Mice population update
  if (mice[i-1] < 0) {
    mice[i] <- 0
  } else {
    mice[i] <- mice[i-1] + (Mouse_r * mice[i-1]) - (beta1 * mice[i-1] * snakes[i-1]) - (beta2 * mice[i-1] * hawks[i-1])
  }
  
  # Snakes population update
  if (snakes[i-1] < 0) {
    snakes[i] <- 0
  } else {
    snakes[i] <- (delta1 * beta1 * mice[i-1] * snakes[i-1]) - (beta3 * snakes[i-1] * hawks[i-1])
  }
  
  # Hawks population update
  if (hawks[i-1] < 0) {
    hawks[i] <- 0
  } else {
    hawks[i] <- (delta2 * beta3 * snakes[i-1] * hawks[i-1]) + (delta3 * beta2 * mice[i-1] * hawks[i-1]) - (gamma * hawks[i-1])
  }
}

# Create a data frame to track populations over time
Tracking_Populations <- data.frame(
  "Mouse_Population" = mice,
  "Snake_Population" = snakes,
  "Hawk_Population" = hawks,
  "Time" = 1:100
)

# Plot the populations over time
plot(Tracking_Populations$Time, Tracking_Populations$Mouse_Population, type="l", col="blue", 
     ylim=c(0, max(c(mice, snakes, hawks), na.rm=TRUE)), xlab="Time (Days)", ylab="Population", 
     main="Population Dynamics of Mice, Snakes, and Hawks")
lines(Tracking_Populations$Time, Tracking_Populations$Snake_Population, col="green")
lines(Tracking_Populations$Time, Tracking_Populations$Hawk_Population, col="red")
legend("topright", legend=c("Mice", "Snakes", "Hawks"), col=c("blue", "green", "red"), lty=1)

###Part 6 If you have not already, create a plot for your three populations based upon your code in part 5. Run your code several times to look at how things change with each run (between the random growth rate for mice and storms, your plots should look different). Do any of the populations crash in any of your runs? If so, which and how does this impact the other populations. What interventions might you take in order to maintain more stable populations and how would this be reflected in your code.

So far in every single run the snakes and hawks have crashed. Hawks crash immediately and snakes crash soon after. When both of these crash obviously the mouse population gains a big boost, which varies from run to run. Sometimes the mouse population will gain a big boost once they crash, go up a bit, then down, then maybe a big boost up again, and then down. Sometimes at the end it will crash, and other times it will keep going up until it hits the end of the plot.

For example, if you created secure burrows that mitigated the risk of storms for snakes, the death rate from storms might change from 20% to 5%. What would this change look like in your code?

I am going to give snakes more secure burrows that make their mortality rate from the storm 7%.It honestly didn’t change anything at all.

# Initialize populations
mice <- c()
snakes <- c()
hawks <- c()

# Starting populations
Mouse_N0 <- 100
Snake_N0 <- 20
Hawk_N0 <- 5

# Growth rates and parameters
Mouse_r <- 0.15
beta1 <- 0.03  # Encounter rate for mice and snakes
beta2 <- 0.01  # Encounter rate for mice and hawks
beta3 <- 0.02  # Encounter rate for snakes and hawks
delta1 <- 0.25 # Conversion rate of mice to snakes
delta2 <- 0.3  # Conversion rate of snakes to hawks
delta3 <- 0.15 # Conversion rate of mice to hawks
gamma <- 0.3   # Natural death rate for hawks

# Simulation loop
for (i in 2:100) {
  # Initial populations
  mice[1] <- Mouse_N0
  snakes[1] <- Snake_N0
  hawks[1] <- Hawk_N0
  # Determine if a storm occurs (on average every 5 years, but randomly)
  Storm_Chance <- runif(1, 0, 1)
  if (Storm_Chance < 0.20) {  # Adjusted probability to match storm frequency
    mice[i-1] <- mice[i-1] * 0.75  # 25% mortality for mice
    snakes[i-1] <- snakes[i-1] * 0.93  # Changed 20% mortality for snakes to 7%
    hawks[i-1] <- hawks[i-1] * 0.95  # 5% mortality for hawks
  }
  # Stochastic mouse growth rate
  Mouse_r <- runif(1, -0.01, 0.09)
  # Mice population update
  if (mice[i-1] < 0) {
    mice[i] <- 0
  } else {
    mice[i] <- mice[i-1] + (Mouse_r * mice[i-1]) - (beta1 * mice[i-1] * snakes[i-1]) - (beta2 * mice[i-1] * hawks[i-1])
  }
  
  # Snakes population update
  if (snakes[i-1] < 0) {
    snakes[i] <- 0
  } else {
    snakes[i] <- (delta1 * beta1 * mice[i-1] * snakes[i-1]) - (beta3 * snakes[i-1] * hawks[i-1])
  }
  
  # Hawks population update
  if (hawks[i-1] < 0) {
    hawks[i] <- 0
  } else {
    hawks[i] <- (delta2 * beta3 * snakes[i-1] * hawks[i-1]) + (delta3 * beta2 * mice[i-1] * hawks[i-1]) - (gamma * hawks[i-1])
  }
}

# Create a data frame to track populations over time
Tracking_Populations <- data.frame(
  "Mouse_Population" = mice,
  "Snake_Population" = snakes,
  "Hawk_Population" = hawks,
  "Time" = 1:100
)

# Plot the populations over time
plot(Tracking_Populations$Time, Tracking_Populations$Mouse_Population, type="l", col="blue", 
     ylim=c(0, max(c(mice, snakes, hawks), na.rm=TRUE)), xlab="Time (Days)", ylab="Population", 
     main="Population Dynamics of Mice, Snakes, and Hawks")
lines(Tracking_Populations$Time, Tracking_Populations$Snake_Population, col="green")
lines(Tracking_Populations$Time, Tracking_Populations$Hawk_Population, col="red")
legend("topright", legend=c("Mice", "Snakes", "Hawks"), col=c("blue", "green", "red"), lty=1)