##Importing

Model1 <- read.csv("Model_1_r0.1.csv")
Model2 <- read.csv("Model_1_r0.5.csv")
Model3 <- read.csv("Model_1_r1.0.csv")

##Analysis

class(Model1)
## [1] "data.frame"
class(Model2)
## [1] "data.frame"
class(Model3)
## [1] "data.frame"
View(Model1)
View(Model2)
View(Model3)
head1 <- head(Model1, n = 3)
head1
##   ï..time parasites immune_cells
## 1     0.0  1.000000     1.000000
## 2     0.1  1.009900     1.000100
## 3     0.2  1.019898     1.000201
tail1 <- tail(Model1)
tail1
##     ï..time parasites immune_cells
## 496    49.5  128.8982     3.397924
## 497    49.6  130.1434     3.436722
## 498    49.7  131.4001     3.476298
## 499    49.8  132.6685     3.516672
## 500    49.9  133.9485     3.557862
## 501    50.0  135.2403     3.599890
sum1 <- summary(Model1)
sum1
##     ï..time       parasites        immune_cells  
##  Min.   : 0.0   Min.   :  1.000   Min.   :1.000  
##  1st Qu.:12.5   1st Qu.:  3.426   1st Qu.:1.025  
##  Median :25.0   Median : 11.728   Median :1.114  
##  Mean   :25.0   Mean   : 27.625   Mean   :1.380  
##  3rd Qu.:37.5   3rd Qu.: 40.057   3rd Qu.:1.473  
##  Max.   :50.0   Max.   :135.240   Max.   :3.600
max1 <- max(Model1$parasites)
max1
## [1] 135.2403
maxpara1 <- which(Model1$parasites == max(Model1$parasites))
maxpara1
## [1] 501

Question: Physiologically, why might the authors make this assumption?

Answer: The assumption that immune cell growth rate plateaus when parasite population grows sufficiently large, can be physiologically attributed to resource competition. After parasite invasion, resources in the host are gradually depleted until little is left for immune cells to access. Other physiological reasons may include: maximum possible immune cell turnover or immunity production.

Question: Why did we exponentiate the second term in the first equation?

Answer: Since the second term is negative, it represents exponential decay. The exponent accounts for immune cell population and its capable rate of destruction for each time step. The output must be subtracted from 1 to produce a rate of decrease.

##Plots!

library('ggplot2')

ggplot(Model1, aes(x= Model1$ï..time)) +
  geom_line(aes(y = Model1$parasites, color = "green"), size = 1.5) +
  geom_line(aes(y = Model1$immune_cells, color = "blue"), size =1.5) +
  theme_classic() + 
  labs(x = "Time", y = "Abundance", title = "Within-Host Population Dynamics for Model 1, r = 0.1") +
  scale_color_identity(name = "Legend",
                       breaks = c("green", "blue"),
                       labels = c("Parasites", "Immune Cells"),
                       guide = "legend")
## Warning: Use of `Model1$parasites` is discouraged. Use `parasites` instead.
## Warning: Use of `Model1$ï..time` is discouraged. Use `ï..time` instead.
## Warning: Use of `Model1$immune_cells` is discouraged. Use `immune_cells`
## instead.
## Warning: Use of `Model1$ï..time` is discouraged. Use `ï..time` instead.

ggplot(Model2, aes(x= Model2$ï..time)) +
  geom_line(aes(y = Model2$parasites, color = "green"), size = 1.5) +
  geom_line(aes(y = Model2$immune_cells, color = "blue"), size =1.5) +
  theme_classic() + 
  labs(x = "Time", y = "Abundance", title = "Within-Host Population Dynamics for Model 2, r = 0.5") +
  scale_color_identity(name = "Legend",
                       breaks = c("green", "blue"),
                       labels = c("Parasites", "Immune Cells"),
                       guide = "legend")
## Warning: Use of `Model2$parasites` is discouraged. Use `parasites` instead.
## Warning: Use of `Model2$ï..time` is discouraged. Use `ï..time` instead.
## Warning: Use of `Model2$immune_cells` is discouraged. Use `immune_cells`
## instead.
## Warning: Use of `Model2$ï..time` is discouraged. Use `ï..time` instead.
## Warning: Removed 4 row(s) containing missing values (geom_path).

## Warning: Removed 4 row(s) containing missing values (geom_path).

ggplot(Model3, aes(x= Model3$ï..time)) +
  geom_line(aes(y = Model3$parasites, color = "green"), size = 1.5) +
  geom_line(aes(y = Model3$immune_cells, color = "blue"), size =1.5) +
  theme_classic() + 
  labs(x = "Time", y = "Abundance", title = "Within-Host Population Dynamics for Model 3, r = 1.0") +
  scale_color_identity(name = "Legend",
                       breaks = c("green", "blue"),
                       labels = c("Parasites", "Immune Cells"),
                       guide = "legend")
## Warning: Use of `Model3$parasites` is discouraged. Use `parasites` instead.
## Warning: Use of `Model3$ï..time` is discouraged. Use `ï..time` instead.
## Warning: Use of `Model3$immune_cells` is discouraged. Use `immune_cells`
## instead.
## Warning: Use of `Model3$ï..time` is discouraged. Use `ï..time` instead.

Question: How might we interpret these data? Do these data seem to represent a chronic infection? An infection that can be cleared? An infection that ends with the death of the host? Explain.

Two of my simulations (r=0.5,1.0) resulted in the infection being cleared. My weakest simulation (r=0.1) likely required more time for immune cell growth. If intrinsic growth rate of parasite is higher, the infection may ultimately kill the host; however, my three simulations cannot disprove this idea. Assuming our discrete model replicates the original differential, the original paper by Antia et al. indicate host death at higher r values.

Another Question: How does the value of “r” change your time series? How might you interpret these dynamical changes physiologically?

Higher values of “r” evoke faster immune cell population growth. Physiologically, higher concentrations of toxins in the body produce a stronger immune response. This is simply modeled using parsite density.

Iterations

for (i in sequence){
  perform some operation 
}
vec <- c() # Makes an empty vector to store the outputs of our for loop in
for (i in 1:10){
  vec[i] <- i + 1 # The value at position i in our vector is i + 1
}
# Means that the value at position 1 in our vector is 1 + 1, the value at position 2 in our vectos is 2 + 1, etc.

print(vec) # To show our resulting vector
##  [1]  2  3  4  5  6  7  8  9 10 11
if (test_expression){
  statement
}

if (test_expression){
  statement
} else{
  statement2
}
x <- 995 # Assign x some value
if (x > 2){
  print("It's a bones day")
}
## [1] "It's a bones day"
if (x == 995){
  print("It's a bones day!")
} else {
  print("It's a no bones day!")
}
## [1] "It's a bones day!"
vec2 <- c(seq(1:10)) # Makes a vector with values 1-10
vec3 <- c() # Another empty vector to store the outputs of our for loop in 
for (i in vec2){
  if (vec2[i] <= 5){
    vec3[i] <- "It's a bones day!"
  } else {
    vec3[i] <- "It's a no bones day!"
  }
}
print(vec3)
##  [1] "It's a bones day!"    "It's a bones day!"    "It's a bones day!"   
##  [4] "It's a bones day!"    "It's a bones day!"    "It's a no bones day!"
##  [7] "It's a no bones day!" "It's a no bones day!" "It's a no bones day!"
## [10] "It's a no bones day!"
function_name <- function(argument){
  statement
}
Eqn_Func <- function(x){
  y <- 2*x + 4
  return(y)
}

# The function needs us to give a value of x. If we say x = 2, then the function calculates y according to the statment we wrote, y = 2x + 4
Eqn_Func(2)
## [1] 8
# This function takes initial parasite and immune cell population abundances, and, according to the model we've been using, determines the population abundances at the next time step
Sim_Func <- function(Para_Last, Immune_Last){
  # Parameter values
  ts = 0.1
  r = 0.3; k = 0.001; p = 1; o = 1000
  
  # A discrete model of within-host infectious disease dynamics
  Para_Next = (1 + r*ts - (1-exp(-k*Immune_Last*ts)))*Para_Last
  Immune_Next = (1 + p*(Para_Last/(Para_Last + o))*ts)*Immune_Last
  
  Para_Last <- Para_Next
  Immune_Last <- Immune_Next
  
  return(c(Para_Last, Immune_Last))
}

# To project this model into the future, we need to implement the function within a for loop that can feed the function the return values from the previous iteration of the loop. 
Abundances <- data.frame() # Empty data frame to hold our outputs
Para_Last <- 1; Immune_Last <- 1 # Initial population abundances
for (i in 1:500){ # Loop will run the function 500 times, but with new "initial" abundances every time to simulate the populations forwards in time
  Output = Sim_Func(Para_Last, Immune_Last)
  Para_Last = Output[1]
  Immune_Last = Output[2]
  Addition <- c(Para_Last, Immune_Last)
  Abundances <- rbind(Abundances, Addition)
  colnames(Abundances) <- c("Parasites", "Immune_Cells")
}

# Let's see if that worked:
head(Abundances)
##   Parasites Immune_Cells
## 1  1.029900     1.000100
## 2  1.060694     1.000203
## 3  1.092409     1.000309
## 4  1.125072     1.000418
## 5  1.158711     1.000530
## 6  1.193357     1.000646
# Let's plot our output!
Inits <- data.frame("1", "1")
colnames(Inits) <- c("Parasites", "Immune_Cells")
Abundances <- data.frame(rbind(Inits, Abundances))
Time <- seq(0,50,0.1)
df <- cbind(Time, Abundances)
plot(df$Time, df$Parasites, type = "l", col = "goldenrod1", lwd = 4, xlab="Time", ylab="Abundance")
lines(df$Time, df$Immune_Cells, col = "cornflowerblue", lwd = 4)
legend("topleft", legend=c("Parasites", "Immune Cells"),
       col=c("goldenrod1", "cornflowerblue"), lty = 1, lwd = 4)

Marking:

  1. Plot with new parameters??? And some aesthetic changes
  2. Answer all interpretation questions in EEB319IntroR.Rmd (highlighted in red)