Homework 1 Course 71256
Yuval Neumann 318375631
2024-12-02
Question 1
1) Updated Flu Model
2)
We have three populations:
- \(S(t)\): People without the flu (susceptible population)
- \(N(t)\): People with the flu (infected population)
- \(I(t)\): Recovered people who are temporarily immune but lose immunity over time.
For \(S(t)\) (Susceptibles):
Susceptibles decrease due to infection at rate \(a \cdot c \cdot S(t) \cdot N(t)\) and increase due to loss of immunity at rate \(\nu \cdot I(t)\): \[ \frac{dS}{dt} = \nu I(t) - a \cdot c \cdot S(t) \cdot N(t) \]
For \(N(t)\) (Infected):
Infected increase due to infection at rate \(a \cdot c \cdot S(t) \cdot N(t)\) and decrease due to recovery at rate \(r \cdot N(t)\): \[ \frac{dN}{dt} = a \cdot c \cdot S(t) \cdot N(t) - r \cdot N(t) \]
For \(I(t)\) (Immune):
Immune increase due to recovery at rate \(r \cdot N(t)\) and decrease due to loss of immunity at rate \(\nu \cdot I(t)\): \[ \frac{dI}{dt} = r \cdot N(t) - \nu I(t) \]
3)
I’m solving the equations with the ‘deSolve’ package in R. All my code is written in R - I’m more comfortable with it.
# Load required packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(deSolve)
# Parameters
a <- 0.1 # Infection rate
c <- 0.5 # Contact rate
r <- 0.3 # Recovery rate
nu <- 0.1 # Immunity loss rate
# Initial conditions - all three populations' fractions must add up to 1
S0 <- 0.7 # Initial susceptible population
N0 <- 0.3 # Initial infected population - I like 0.3, it makes the figure more interesting
I0 <- 0.0 # Initial immune population - everyone's still sick at the beginning
y0 <- c(S = S0, N = N0, I = I0)
# Making sure they sum to 1
(S0 + N0 + I0) == 1## [1] TRUE# Time points
t <- 1:100
# Differential equations
flu_model <- function(t, y, params) {
S <- y["S"]
N <- y["N"]
I <- y["I"]
with(as.list(params), {
dS <- nu * I - a * c * S * N
dN <- a * c * S * N - r * N
dI <- r * N - nu * I
list(c(dS, dN, dI))
})
}
# Parameters for the model
params <- list(a = a, c = c, r = r, nu = nu)
# Solve the ODEs
solution <- ode(y = y0, times = t, func = flu_model, parms = params)
# Convert the solution to a data frame
solution_df <- as.data.frame(solution)
# Pivot the DF longer so I can plot it conveniently
solution_df_long <- solution_df %>%
pivot_longer(cols = 2:4,
names_to = "Population",
values_to = "frac")
# Factoring the population column so it's easier to plot
solution_df_long$Population <- as.factor(solution_df_long$Population)
# Plotting
ggplot(solution_df_long, aes(x = time, y = frac, color = Population)) + # Basic command
theme_bw() + # Visually pleasing theme
geom_line(lwd = 1.2, alpha = 0.7) + # Adding the line - it will be subsetted as defined in ggplot
scale_color_manual(values = c("I" = "navy",
"N" = "forestgreen",
"S" = "firebrick3"),
labels = c("I(t)", "N(t)", "S(t)")) + # Custom colors and legend labels
labs(x = "Time",
y = "Population Fraction",
title = "Flu Model Dynamics") + # Setting the title and axis titles
scale_y_continuous(breaks = seq(0, 1, 0.2)) + # Making the y-axis more sensible
theme(aspect.ratio = 0.75,# Ratio of the figure - I find continuous figures to benefit from a 0.75 ratio
plot.title = element_text(size = 28, face = "bold"),
legend.title = element_text(size = 26, face = "bold"),
legend.text = element_text(size = 24),
axis.title = element_text(size = 26),
axis.text = element_text(size = 24)) # A bunch of text adjustments to make the figure readable
Question 2
1)
- \(X(t)\): Reproductive population, since it’s the only one out of the three that has an independent growth ( \(b X(t)\) ).
- \(Y(t)\): Workers, since the soldiers are a specialized group of workers, that means the workers must come directly from the reproductive population, hence \(k X(t)\).
- \(Z(t)\): The soldiers, which are derived directly from the workers (\(Y(t)\)) at a rate of \(a Z(t) \frac{Y(t)}{X(t)+Y(t)+Z(t)}\).
2)
Reproductive Population (\(X(t)\)): \[ \frac{dX(t)}{dt} = b X(t) + \theta - k X(t) \]
Workers Population (\(Y(t)\)): \[ \frac{dY(t)}{dt} = k X(t) + g Z(t) - a Z(t)\left[ \frac{Y(t)}{X(t)+Y(t)+Z(t)} \right] \] Let’s define \(\frac{Y(t)}{X(t)+Y(t)+Z(t)}\) as \(\gamma\): \[ \frac{dY(t)}{dt} = k X(t) + g Z(t) - a Z(t)\left[ \gamma \right] = \] \[ \frac{dY(t)}{dt} = k X(t) + Z(t)\left[g - a\gamma \right] \]
Soldiers Population (\(Z(t)\)): \[ \frac{dZ(t)}{dt} = a Z(t)\left[ \frac{Y(t)}{X(t)+Y(t)+Z(t)} \right] - d Z(t) - g Z(t) \] Again I define \(\frac{Y(t)}{X(t)+Y(t)+Z(t)}\) as \(\gamma\): \[ \frac{dZ(t)}{dt} = a Z(t)\left[ \gamma \right] - d Z(t) - g Z(t) = \] \[ \frac{dZ(t)}{dt} = Z(t)\left[ a \gamma - d - g \right] \]
Question 3
1)
# Opening the .csv file
rachel_data <- read.csv("Rachel_experiment.csv")
# Looking at the data
head(rachel_data)## plant_1 plant_2 plant_3 plant_4 plant_5 plant_6 plant_7 plant_8
## 1 4.346797 4.219714 4.319874 4.211997 4.391348 4.494053 4.351277 4.166667
## 2 4.349052 4.222344 4.323806 4.211521 4.395961 4.493523 4.350097 4.167167
## 3 4.349159 4.223810 4.322862 4.211082 4.393910 4.494502 4.354086 4.165323
## 4 4.351118 4.222565 4.321319 4.210129 4.393380 4.493495 4.355966 4.169929
## 5 4.352799 4.222877 4.319346 4.212100 4.391878 4.492605 4.355026 4.171899
## 6 4.354834 4.224095 4.323825 4.212798 4.390542 4.495211 4.354023 4.169921
## plant_9 plant_10
## 1 4.424110 4.301932
## 2 4.424865 4.304454
## 3 4.425007 4.304397
## 4 4.423972 4.303569
## 5 4.422331 4.301668
## 6 4.424594 4.304116# Now I'll add a time column
# First, the number of hours in total
time_of_experiment <- nrow(rachel_data)*12 - 12
# Now adding the column
rachel_data$time <- seq(0, time_of_experiment, by = 12)
# Another column, in units of weeks
hours_in_day <- 24
days_in_week <- 7
hours_in_week <- hours_in_day * days_in_week
rachel_data$time_weeks <- rachel_data$time/hours_in_week
# Showing my final column
rachel_data$time_weeks## [1] 0.00000000 0.07142857 0.14285714 0.21428571 0.28571429 0.35714286
## [7] 0.42857143 0.50000000 0.57142857 0.64285714 0.71428571 0.78571429
## [13] 0.85714286 0.92857143 1.00000000 1.07142857 1.14285714 1.21428571
## [19] 1.28571429 1.35714286 1.42857143 1.50000000 1.57142857 1.64285714
## [25] 1.71428571 1.78571429 1.85714286 1.92857143 2.00000000 2.07142857
## [31] 2.14285714 2.21428571 2.28571429 2.35714286 2.42857143 2.50000000
## [37] 2.57142857 2.64285714 2.71428571 2.78571429 2.85714286 2.92857143
## [43] 3.00000000 3.07142857 3.14285714 3.21428571 3.28571429 3.35714286
## [49] 3.42857143 3.50000000 3.57142857 3.64285714 3.71428571 3.78571429
## [55] 3.85714286 3.92857143 4.00000000 4.07142857 4.14285714 4.21428571
## [61] 4.28571429 4.35714286 4.42857143 4.50000000 4.57142857 4.64285714
## [67] 4.71428571 4.78571429 4.85714286 4.92857143 5.00000000 5.07142857
## [73] 5.14285714 5.21428571 5.28571429 5.35714286 5.42857143 5.50000000
## [79] 5.57142857 5.64285714 5.71428571 5.78571429 5.85714286 5.92857143
## [85] 6.00000000 6.07142857 6.14285714 6.21428571 6.28571429 6.35714286
## [91] 6.42857143 6.50000000 6.57142857 6.64285714 6.71428571 6.78571429
## [97] 6.85714286 6.92857143 7.00000000 7.07142857 7.14285714 7.21428571
## [103] 7.28571429 7.35714286 7.42857143 7.50000000 7.57142857 7.64285714
## [109] 7.71428571 7.78571429 7.85714286 7.92857143 8.00000000 8.07142857
## [115] 8.14285714 8.21428571 8.28571429 8.35714286 8.42857143 8.50000000
## [121] 8.57142857 8.64285714 8.71428571 8.78571429 8.85714286 8.92857143
## [127] 9.00000000 9.07142857 9.14285714 9.21428571 9.28571429 9.35714286
## [133] 9.42857143 9.50000000 9.57142857 9.64285714 9.71428571 9.78571429
## [139] 9.85714286 9.92857143 10.00000000 10.07142857 10.14285714 10.21428571
## [145] 10.28571429 10.35714286 10.42857143 10.50000000 10.57142857 10.64285714
## [151] 10.71428571 10.78571429 10.85714286 10.92857143 11.00000000 11.07142857
## [157] 11.14285714 11.21428571 11.28571429 11.35714286 11.42857143 11.50000000
## [163] 11.57142857 11.64285714 11.71428571 11.78571429 11.85714286 11.92857143
## [169] 12.00000000 12.07142857 12.14285714 12.21428571 12.28571429 12.35714286
## [175] 12.42857143 12.50000000 12.57142857 12.64285714 12.71428571 12.78571429
## [181] 12.85714286 12.92857143 13.00000000 13.07142857 13.14285714 13.21428571
## [187] 13.28571429 13.35714286 13.42857143 13.50000000 13.57142857 13.64285714
## [193] 13.71428571 13.78571429 13.85714286 13.92857143 14.00000000 14.07142857
## [199] 14.14285714 14.214285712)
# I need to subtract the initial weight of each column from all the column
# I'll make a new DF and then subtract manually, will take very little code and less messy than a loop (since it's only 10 columns)
normalized_df <- rachel_data
normalized_df$plant_1 <- normalized_df$plant_1 - normalized_df$plant_1[1] # In R, indexing starts at 1
normalized_df$plant_2 <- normalized_df$plant_2 - normalized_df$plant_2[1]
normalized_df$plant_3 <- normalized_df$plant_3 - normalized_df$plant_3[1]
normalized_df$plant_4 <- normalized_df$plant_4 - normalized_df$plant_4[1]
normalized_df$plant_5 <- normalized_df$plant_5 - normalized_df$plant_5[1]
normalized_df$plant_6 <- normalized_df$plant_6 - normalized_df$plant_6[1]
normalized_df$plant_7 <- normalized_df$plant_7 - normalized_df$plant_7[1]
normalized_df$plant_8 <- normalized_df$plant_8 - normalized_df$plant_8[1]
normalized_df$plant_9 <- normalized_df$plant_9 - normalized_df$plant_9[1]
normalized_df$plant_10 <- normalized_df$plant_10 - normalized_df$plant_10[1]
# Showing the first few rows:
head(normalized_df)## plant_1 plant_2 plant_3 plant_4 plant_5
## 1 0.000000000 0.000000000 0.0000000000 0.0000000000 0.000000000
## 2 0.002255426 0.002629493 0.0039316009 -0.0004762568 0.004612573
## 3 0.002362859 0.004095678 0.0029873360 -0.0009152076 0.002562236
## 4 0.004321315 0.002850630 0.0014446686 -0.0018683082 0.002031779
## 5 0.006002816 0.003163091 -0.0005283232 0.0001032461 0.000529461
## 6 0.008037703 0.004380937 0.0039504147 0.0008006985 -0.000805877
## plant_6 plant_7 plant_8 plant_9 plant_10 time
## 1 0.0000000000 0.000000000 0.0000000000 0.0000000000 0.0000000000 0
## 2 -0.0005294736 -0.001179955 0.0005002467 0.0007549990 0.0025221866 12
## 3 0.0004490139 0.002809636 -0.0013434276 0.0008963889 0.0024649058 24
## 4 -0.0005579254 0.004689435 0.0032622556 -0.0001378918 0.0016367117 36
## 5 -0.0014472897 0.003749088 0.0052325101 -0.0017796126 -0.0002637369 48
## 6 0.0011578310 0.002746338 0.0032539274 0.0004841392 0.0021843563 60
## time_weeks
## 1 0.00000000
## 2 0.07142857
## 3 0.14285714
## 4 0.21428571
## 5 0.28571429
## 6 0.357142863)
# I need to multiply all the values in all plant columns
# I'll multiply manually, again, not worth it to write a loop for 10 columns
oz_to_grams <- 28.35
normalized_df$plant_1 <- normalized_df$plant_1 * oz_to_grams
normalized_df$plant_2 <- normalized_df$plant_2 * oz_to_grams
normalized_df$plant_3 <- normalized_df$plant_3 * oz_to_grams
normalized_df$plant_4 <- normalized_df$plant_4 * oz_to_grams
normalized_df$plant_5 <- normalized_df$plant_5 * oz_to_grams
normalized_df$plant_6 <- normalized_df$plant_6 * oz_to_grams
normalized_df$plant_7 <- normalized_df$plant_7 * oz_to_grams
normalized_df$plant_8 <- normalized_df$plant_8 * oz_to_grams
normalized_df$plant_9 <- normalized_df$plant_9 * oz_to_grams
normalized_df$plant_10 <- normalized_df$plant_10 * oz_to_grams
# Showing the first few rows:
head(normalized_df)## plant_1 plant_2 plant_3 plant_4 plant_5 plant_6
## 1 0.00000000 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
## 2 0.06394132 0.07454612 0.11146089 -0.013501880 0.13076644 -0.01501058
## 3 0.06698706 0.11611249 0.08469097 -0.025946136 0.07263940 0.01272954
## 4 0.12250927 0.08081537 0.04095636 -0.052966538 0.05760094 -0.01581719
## 5 0.17017984 0.08967363 -0.01497796 0.002927026 0.01501022 -0.04103066
## 6 0.22786889 0.12419958 0.11199426 0.022699801 -0.02284661 0.03282451
## plant_7 plant_8 plant_9 plant_10 time time_weeks
## 1 0.00000000 0.00000000 0.000000000 0.00000000 0 0.00000000
## 2 -0.03345173 0.01418199 0.021404223 0.07150399 12 0.07142857
## 3 0.07965318 -0.03808617 0.025412626 0.06988008 24 0.14285714
## 4 0.13294547 0.09248495 -0.003909231 0.04640078 36 0.21428571
## 5 0.10628665 0.14834166 -0.050452017 -0.00747694 48 0.28571429
## 6 0.07785868 0.09224884 0.013725348 0.06192650 60 0.357142864)
# I'll pivot my DF longer, to have an easier time with plotting
long_normalized <- normalized_df %>%
pivot_longer(cols = 1:10,
names_to = "plant",
values_to = "weight")
long_normalized$plant <- factor(long_normalized$plant,
levels = c("plant_1",
"plant_2",
"plant_3",
"plant_4",
"plant_5",
"plant_6",
"plant_7",
"plant_8",
"plant_9",
"plant_10"))
ggplot(long_normalized, aes(x = time_weeks, y = weight)) +
theme_bw() +
geom_line(color = "firebrick3", lwd = 1.2) +
scale_x_continuous(breaks = seq(0, 14, by = 2)) +
facet_wrap(~plant, ncol = 5, nrow = 2) +
labs(x = "Time (Weeks)",
y = "Weight (grams)") +
theme(strip.text = element_text(size = 26, face = "bold"),
axis.title = element_text(size = 26),
axis.text = element_text(size = 24))
5)
I need to keep plants 5, 6, 9 and 10. Let’s use the ‘filter’ function from dplyr:
filtered_data <- long_normalized %>%
filter(plant %in% c("plant_5", "plant_6", "plant_9", "plant_10"))6)
Let’s plot:
ggplot(filtered_data, aes(x = time_weeks, y = weight, color = plant)) +
theme_bw() +
geom_line(lwd = 1.2) +
scale_x_continuous(breaks = seq(0, 14, by = 1)) +
labs(x = "Time (Weeks)",
y = "Weight (grams)",
color = "Plant Number") +
scale_color_manual(values = c("plant_5" = "#F87A53",
"plant_6" = "#E6C767",
"plant_9" = "#898121",
"plant_10" = "#4C4B16"),
labels = c("5", "6", "9", "10")) +
theme(axis.title = element_text(size = 26),
axis.text = element_text(size = 24),
legend.title = element_text(size = 26, face = "bold"),
legend.text = element_text(size = 24)) 
I’ve decided to calculate \(EC_{50}\), as it’s the concentration (here - number of weeks) where half the growth happens. I’ll use the ‘drc’ package:
library(drc)
library(broom)
model_ec50 <- drm(data = filtered_data,
formula = weight ~ time_weeks,
curveid = plant,
fct = LL.4(names = c("Hill slope", "Min", "Max", "EC50")))
summary_df <- as.data.frame(tidy(model_ec50))
summary_df[summary_df$term == "EC50" , ]## term curve estimate std.error statistic p.value
## 13 EC50 plant_5 7.663238 0.01550401 494.2746 0
## 14 EC50 plant_6 7.094828 0.01639223 432.8166 0
## 15 EC50 plant_9 9.241837 0.01635726 564.9992 0
## 16 EC50 plant_10 9.036597 0.01452617 622.0910 0After calculating \(EC_{50}\), I look at the data of each plant, and take the weight value where my \(EC_{50}\) line should end as the weight value of the nearest time value.
ec50_df <- summary_df %>% filter(term == "EC50")
plant_5_df <- normalized_df[ , c("plant_5", "time_weeks")]
plant_6_df <- normalized_df[ , c("plant_6", "time_weeks")]
plant_9_df <- normalized_df[ , c("plant_9", "time_weeks")]
plant_10_df <- normalized_df[ , c("plant_10", "time_weeks")]
ggplot(filtered_data, aes(x = time_weeks, y = weight, color = plant)) +
theme_bw() +
geom_line(lwd = 1.2) +
geom_segment(aes(x = 7.663238, y = 0, xend = 7.663238, yend = 7.84564088), color = "#F87A53", lwd = 1.2) +
geom_segment(aes(x = 7.094828, y = 0, xend = 7.094828, yend = 7.22206923), color = "#E6C767", lwd = 1.2) +
geom_segment(aes(x = 9.241837, y = 0, xend = 9.241837, yend = 7.8118574235), color = "#898121", lwd = 1.2) +
geom_segment(aes(x = 9.036597, y = 0, xend = 9.036597, yend = 8.3022362482), color = "#4C4B16", lwd = 1.2) +
labs(x = "Time (Weeks)",
y = "Weight (grams)",
color = "Plant Number") +
scale_x_continuous(breaks = seq(0, 14, by = 1)) +
scale_color_manual(values = c("plant_5" = "#F87A53",
"plant_6" = "#E6C767",
"plant_9" = "#898121",
"plant_10" = "#4C4B16"),
labels = c("5", "6", "9", "10")) +
theme(axis.title = element_text(size = 26),
axis.text = element_text(size = 24),
legend.title = element_text(size = 26, face = "bold"),
legend.text = element_text(size = 24))
Question 4
1) Solve for \(n(1)\)
The recursion equation for the logistic model is:
\[ n(t+1) = n(t) + r \cdot n(t) \left( 1 - \frac{n(t)}{K} \right), \]
where \(r = 1\) and \(K = 100\). Using the initial values:
For \(n(0) = 50\): \[ n(1) = 50 + 1 \cdot 50 \left( 1 - \frac{50}{100} \right) = 50 + 50 \cdot \frac{1}{2} = 50 + 25 = 75. \]
For \(n(0) = 100\): \[ n(1) = 100 + 1 \cdot 100 \left( 1 - \frac{100}{100} \right) = 100 + 100 \cdot 0 = 100. \]
For \(n(0) = 200\): \[ n(1) = 200 + 1 \cdot 200 \left( 1 - \frac{200}{100} \right) = 200 + 200 \cdot (-1) = 200 - 200 = 0. \]
For \(n(0) = 500\): \[ n(1) = 500 + 1 \cdot 500 \left( 1 - \frac{500}{100} \right) = 500 + 500 \cdot (-4) = 500 - 2000 = -1500. \]
2) Determine \(n^*\) Above Which \(n(t+1) < 0\)
The recursion equation becomes negative when:
\[ n(t+1) = n(t) + r \cdot n(t) \left( 1 - \frac{n(t)}{K} \right) < 0. \]
Rearranging for \(n(t)\): \[ r \cdot n(t) \left( 1 - \frac{n(t)}{K} \right) < -n(t). \]
Divide through by \(n(t)\) (assuming \(n(t) > 0\)): \[ r \left( 1 - \frac{n(t)}{K} \right) < -1. \]
Simplify: \[ 1 - \frac{n(t)}{K} < -\frac{1}{r}. \]
\[ -\frac{n(t)}{K} < -\frac{1}{r} - 1. \]
Multiply through by \(-1\) (flipping the inequality): \[ \frac{n(t)}{K} > \frac{1}{r} + 1. \]
Solve for \(n(t)\): \[ n(t) > K \cdot \left( \frac{1}{r} + 1 \right). \]
Thus: \[ n^* = K \cdot \left( \frac{1}{r} + 1 \right). \]
For \(r = 1\) and \(K = 100\): \[ n^* = 100 \cdot \left( \frac{1}{1} + 1 \right) = 100 \cdot 2 = 200. \]
3) Verifying
From Part 1, when \(n(0) > 200\), the population \(n(t+1)\) becomes negative: - For \(n(0) = 500\), we found \(n(1) = -1500\).
This matches the result from Part 2, where \(n^* = 200\) is the threshold above which \(n(t+1)\) becomes negative.