Mid-project update.Rmd

---
title: "Mid-project update"
author: "Rachel Rodriguez, Riho Ishida, and Alicia Magsano"
date: "2024-11-19"
output: pdf_document
---

```{r}
library(tidyverse)
```

# Previous Work/Information

*Batrachochytrium dendrobatidis* (*Bd*) is a fungal pathogen that causes chytridiomycosis which is a disease that is lethal to a large number of amphibian species (Scheele et al., 2019). Chytridiomycosis has had detrimental effects on amphibian species as it has caused global declines and extinctions in many susceptible species (Scheele et al., 2019; Skerratt et al., 2007). As a result, the declines caused by *Bd* have been called, "the greatest recorded loss of biodiversity attributable to a disease" (Scheele et al., 2019). 

Given how severe the effects of *Bd* are it is important to understand what factors can help to limit the effects of this deadly pathogen. Temperature may be one such factor as previous research has shown that *Bd*, in its stationary phase, thrives in colder temperatures of 13-15°C, with a rapid decline of density at around 26°C (Stevenson et al., 2013). The same study demonstrated that the parasite grows best at around 21°C (Stevenson et al., 2013). Furthermore, there is evidence that *Bd* zoospore activity decreases as temperature increases, which correlates with a reduced prevalence of *Bd* in hosts (Sapsford et al., 2013; Stevenson et al., 2013). In vivo experiments support this, as reduced mortality is seen in subjects receiving higher temperatures for treatment (Berger et al., 2004). Therefore, given that *Bd* favours cooler temperatures, it is important to understand how changes in temperature due to climate change may affect *Bd* because this can help us better understand how to facilitate future conservation efforts for amphibians. 

Other studies have used modelling to determine the distribution of *Bd* in a certain location or to predict how *Bd*'s range may shift as a result of climate change (Sun et al., 2023; Tytar et al., 2023). However, there are no studies on how the change in temperature due to climate change may change the dynamics of *Bd* and how it influences susceptible amphibian populations.

# Our Question and Prediction

Question: How will changes in temperature due to climate change affect the probability that a population of frogs suffering from *Bd* infections will go extinct? 

Prediction: Since *Bd* favours cooler temperatures, we predict that at higher temperatures there will be a lower probability that a frog population infected with *Bd* will go extinct. 

# Analysis Plan

To answer this question we created a thermal performance curve (TPC) from data containing zoospore count and temperature. Then we made a stage-structured SIR model for a frog population which we will solve for at different temperatures by using the TPC we created to indicate the initial number of zoospores present at a given temperature. 

## Data Used

The data we used to make the TPC was created by combining the data for three different studies including two field studies and one lab study.

1. First Field Study (Kriger & Hero, 2007)
We used the data collected by Kriger and Hero (2007) from their study on the seasonal variation of chytridiomycosis prevalence and severity in a single population of *Litoria wilcoxii* (eastern stony creek frogs). In this study, they captured and sampled *L. wilcoxii* over 21 months in 6-week intervals within a 1 km area of the Nerang River in Queensland, Australia. During each of the 13 sampling sessions, they swabbed frogs for chytridiomycosis, recorded the air and water temperature at the site, and obtained rainfall data from a nearby weather station. The swabs were tested for *Bd* presence using quantitative polymerase chain reaction (qPCR) techniques which were used to determine the prevalence of *Bd* (divided the number of *Bd* positive frogs by the total number of frogs sampled that session) and the geometric mean number of zoospores (parasite load) found in infected frogs. Additionally, they determined each session's air and water temperature by calculating the mean air temperature and total amount of rain in the 30 days leading up to the sampling day. Finally, since the water temperature was manually sampled at the beginning and end of each sampling session, the mean of these two measurements was used as the water temperature for a given sampling session. Therefore, we used the zoospores and 30-day air temperature (°C) data from this study. 

2. Second Field Study (Kriger et al., 2007)
We also used the data collected by Kriger et al. (2007) whose study focused on how the prevalence and intensity of *Bd* is affected by latitudinal variation. In this study, they captured and swabbed *Litoria lesueuri* (stony creek frogs) for *Bd* at 31 sites along the east coast of Australia which encompassed a gradient of latitudes and the sampling from each site was done within 42 days (20 September 2005 – 1 November 2005). The swabs collected at each site were tested for *Bd* presence and quantity using qPCR. Furthermore, in addition to recording the latitude, longitude, and elevation at each sampling site, Kriger et al. (2007) also recorded different short and long-term climatic conditions. From the climatic conditions recorded, we were interested in the mean temperature of the warmest quarter (MTWQ). Therefore, we used the zoospores per frog (mean number of zoospores of all frogs sampled at a site) and MTWQ (°C) data from this study. 

3. Lab Study (Sheets et al., 2021)
Additionally, we used the data collected by Sheets et al. (2021) who conducted a study to determine how temperature affects different *Bd* isolates. In this study, Sheets et al. (2021) specifically tested isolates from the *Bd* lineage Global Panzootic Lineage (*Bd*GPL) which is hypervirulent (Farrer et al., 2011). In their study, Sheets et al. (2021) collected 5 different *Bd* isolates from amphibians around the United States which they genotyped using an amplicon sequencing approach. They then tested how temperature affected different traits of each isolate including viability, zoospore density, growth rate, and carrying capacity by incubating the isolates at various temperatures (4, 12, 17, 21, 25, 26, and 27 °C). Using this information, they generated thermal performance curves for each genotype and isolate. Therefore, we used the temperature and zoospore count (calculated as the average number of zoospores counted from two sides of a hemocytometer) data from the Temp and Counts columns. However, since the zoospore count was recorded over multiple days for each isolate and temperature combination we only included the highest number of zoospores for each isolate and temperature combination in our combined dataset. We decided to do this because before the highest number of zoospores was reached the number of zoospores was still able to grow and after the highest number was reached the number of zoospores likely began to decline due to not being able to persist in the lab setting anymore. 

## TPC

```{r}
data <- read.csv("Combined Chytrid Data.csv") 

# Max Temperature in the combined dataset
max(data$Temperature)
data$Zoospores <- as.integer(data$Zoospores)

# Plot of the combined zoospore and temperature data
plot(data$Temperature,data$Zoospores, 
     main = "Change in Zoospore Count Across Different Temperatures", 
     xlab = "Temperature (°C)",
     ylab = "Number of Zoospores")

# Distribution of zoospore count in the combined dataset 
hist(data$Zoospores, main = "Distribution of Zoospore Count", 
     xlab = "Number of Zoospores")
```

Since our data had so much variance we used a negative binomial distribution to construct the TPC. 
```{r}
# Function for evaluating the log-likelihood of all observations 
LL_eval <- function(data, A, Topt, w, size){
  
  mean <- round(A*exp(-(data$Temperature-Topt)^2/(2*w^2)))
  
 return(sum(dnbinom(x = data$Zoospores, size = size, mu = mean, log = T)))
  
}

# Testing if the function works 
LL_eval(data, A = 10, Topt = 15, w = 1, size = 1)
```

The parameters of our TPC that we want to find the maximum likelihood estimate for include:
- A: the average max number of zoospores that can be produced at the Topt.
- Topt: the optimum temperature for *Bd* zoospores.
- w: the width around the thermal optimum.
- size: the dispersion coefficient
```{r}
# Parameters to test 
combinations <- expand.grid(A = seq(1000,5000,100), 
                            Topt = seq(10,20,0.1), 
                            w = seq(4, 8, 0.1),
                            size = seq(0.01, 0.75, 0.05))

# Identifying the parameters that are most likely to describe our data
LL <- NULL

for (i in 1:nrow(combinations)){
  LL[[i]] <- LL_eval(data, A = combinations$A[i], 
                     Topt = combinations$Topt[i], 
                     w = combinations$w[i],
                     size = combinations$size[i])
}

LL <- as.data.frame(do.call(rbind, LL))
colnames(LL) <- c("LL")

# The maximum likelihood estimate of each parameter
MLE <- data.frame(combinations, LL) %>% subset(LL == max(LL)) 
MLE
```
Therefore, given the TPC the optimum temperature for *Bd* zoospores is 15°C. This thermal optimum is slightly lower than the thermal optimum generally associated with *Bd* which is 17-25°C (Bradley et al., 2019; Piotrowski et al., 2004). This difference might be a result of the specific data we used to generate the TPC. First, the frog species tested in the field studies are reservoir hosts for chytridiomycosis meaning they can survive with a higher prevalence of the disease compared to more susceptible species (Retallick et al., 2004). Furthermore, the lab data we included comes from different isolates of *Bd* and there is evidence to suggest that different isolates of *Bd* have different thermal optimum (Voyles et al., 2017). Therefore, including the data obtained from different isolates could have also influenced the optimum temperature we obtained. 

```{r}
# TPC over the data
temperature <- seq(0, 30, length = 100)
plot(temperature, MLE$A*exp(-(temperature-MLE$Topt)^2/(2*MLE$w^2)), type = "l", 
     ylim = c(0,max(data$Zoospores)), main = "TPC for Bd Zoospores", 
     xlab = "Temperature (°C)", ylab = "Number of Zoospores", col = "red")
points(data$Temperature,data$Zoospores)
legend(0, 12000, legend = "TPC", fill = "red")
```

Confidence interval for each parameter: 
```{r}
# For Topt
data.frame(combinations, LL) %>% 
  subset(A == MLE$A & w == MLE$w &  size == MLE$size) %>% 
  subset(abs(LL - MLE$LL) < qchisq(0.95,df=1)/2) %>% 
  summarise(upperCI_Topt = max(Topt), lowerCI_Topt = min(Topt))

# For A
data.frame(combinations, LL) %>% 
  subset(Topt == MLE$Topt & w == MLE$w &  size == MLE$size) %>% 
  subset(abs(LL - MLE$LL) < qchisq(0.95,df=1)/2) %>% 
  summarise(upperCI_A = max(A), lowerCI_A = min(A))

# For w
data.frame(combinations, LL) %>% 
  subset(Topt == MLE$Topt & A == MLE$A &  size == MLE$size) %>% 
  subset(abs(LL - MLE$LL) < qchisq(0.95,df=1)/2) %>% 
  summarise(upperCI_w = max(w), lowerCI_w = min(w))

# For size
data.frame(combinations, LL) %>% 
  subset(Topt == MLE$Topt & A == MLE$A &  w == MLE$w) %>% 
  subset(abs(LL - MLE$LL) < qchisq(0.95,df=1)/2) %>% 
  summarise(upperCI_size = max(size), lowerCI_size = min(size))
```

Profile of log-likelihoods for each parameter
```{r}
# For Topt
data.frame(combinations, LL) %>% 
  subset(A == MLE$A & w == MLE$w &  size == MLE$size) %>% 
  subset(is.finite(LL)) %>% ggplot(aes(x = Topt, y = LL)) + geom_point() +
  labs(title = "Topt", y = "log-likelihood") +
  geom_vline(xintercept = MLE$Topt, color = "red") +
  geom_vline(xintercept = 13.2, color = "red", linetype = "dashed") +
  geom_vline(xintercept = 17.2, color = "red", linetype = "dashed")

# For  A
data.frame(combinations, LL) %>% 
  subset(Topt == MLE$Topt & w == MLE$w & size == MLE$size) %>%
  subset(is.finite(LL)) %>% ggplot(aes(x = A, y = LL)) + geom_point() +
  labs(title = "Average Max Number of Zoospores Produced at the Topt", 
        y = "log-likelihood") +
  geom_vline(xintercept = MLE$A, color = "red") +
  geom_vline(xintercept = 1100, color = "red", linetype = "dashed") +
  geom_vline(xintercept = 2700, color = "red", linetype = "dashed")

# For w
data.frame(combinations, LL) %>% 
  subset(A == MLE$A & Topt == MLE$Topt &  size == MLE$size) %>%
  subset(is.finite(LL)) %>% ggplot(aes(x = w, y = LL)) + geom_point() +
  labs(title = "Width Around the Thermal Optimum", y = "log-likelihood") +
  geom_vline(xintercept = MLE$w, color = "red") +
  geom_vline(xintercept = 4.9, color = "red", linetype = "dashed") +
  geom_vline(xintercept = 7.8, color = "red", linetype = "dashed")

# For size
data.frame(combinations, LL) %>%
  subset(A == MLE$A & Topt == MLE$Topt &  w == MLE$w) %>% 
  subset(is.finite(LL)) %>% ggplot(aes(x = size, y = LL)) + geom_point() +
  labs(title = "Dispersion Coefficient", y = "log-likelihood") +
  geom_vline(xintercept = MLE$size, color = "red") +
  geom_vline(xintercept = 0.21, color = "red", linetype = "dashed") +
  geom_vline(xintercept = 0.36, color = "red", linetype = "dashed")
```

## Mathematical Model

We created a stage-structured SIR model including differential equations for zoospores, susceptible and infected tadpoles, and susceptible and infected frogs.

**Zoospores**
$$\frac{d Z}{d t} = -d_z Z + p_f I_f Z + p_t I_t Z$$

**Susceptible/Infected Tadpoles**
$$\frac{d S_t}{d t} = b S_t - \beta_t S_t I_t - \beta_f S_t I_f - \beta_z S_t Z - m_S S_t - d_{St} S_t + r_t I_t$$

$$\frac{d I_t}{d t} = \beta_t S_t I_t + \beta_f S_t I_f + \beta_z S_t Z - d_{It} I_t - v_t I_t - m_I I_t - r_t I_t$$

**Susceptible/Infected Frogs**
$$\frac{d S_f}{d t} = m_S S_t - \beta_t S_f I_t - \beta_f S_f I_f - \beta_z S_f Z - d_{Sf} S_f + r_f I_f$$

$$\frac{d I_f}{d t} = m_I I_t + \beta_t S_f I_t + \beta_f S_f I_f + \beta_z S_f Z - d_{If} I_f - v_f I_f - r_f I_f$$

**List of Variables**
- Z = number of zoospores
- $S_t$ = number of susceptible tadpoles
- $I_t$ = number of infected tadpoles 
- $S_f$ = number of susceptible frogs
- $I_f$ = number of infected frogs
- $p_f$ = production rate of zoospores from infected frogs
- $p_t$ = production rate of zoospores from infected tadpoles
- $m_S$ = maturation rate of susceptible tadpoles
- $m_I$ = maturation rate of infected tadpoles
- b = birth rate of tadpoles
- $r_t$ = recovery rate of tadpoles
- $r_f$ = recovery rate of frogs
- $v_t$ = virulence in tadpoles (death rate caused by the disease)
- $v_f$ = virulence in frogs (death rate caused by the disease)
- $\beta_f$ = rate of infection from frogs
- $\beta_t$ = rate of infection from tadpoles
- $\beta_z$ = rate of infection from environmental zoospores 
- $d_z$ = death rate of zoospores 
- $d_{St}$ = death rate of susceptible tadpoles
- $d_{It}$ = death rate of infected tadpoles 
- $d_{Sf}$ = death rate of susceptible frogs
- $d_{If}$ = death rate of infected frogs 

## Assumptions of the Model

We assume:
- There is a homogeneous mixing of the population, in other words, all infected and susceptible specimens have an equal chance of interaction (Tolles & Luong, 2020). This may not account for ecological barriers or behavioural changes induced by infection. 

- The chance of zoospore infection is equal for all zoospores, which is unrealistic in the real world. The possibility of infection will depend on random chance, the isolate, as well as the geographic location of the zoospores.

- There is no vertical inheritance of the disease. This assumption is in line with observations from field studies (Infection with *Batrachochytrium Dendrobatidis (Bd)*, 2020).

- There’s no emigration or immigration in the population.

- Individuals who recover from the disease do not have immunity and thus can be reinfected. 

- The parameters used in our model are estimates. For example, we assume the maturation rate for infected and susceptible tadpoles (moving from $S_t$ to $S_f$, or $I_t$ to $I_f$) is different because the trade-off between growth and defence will lead infected individuals to allocate more resources to fighting the infection and thus grow slower (Rauw, 2012).

## Limitations of the Model 

- Our model uses an initial number of zoospores from our TPC. The TPC is created using data based on *L. wilcoxii*/*L. lesueuri*, which are reservoir species, and isolates of *Bd* found in the US. This was an attempt to make a generalized estimate of the number of zoospores for any *Bd* isolate from the *Bd*GPL lineage in *L. wilcoxii*/*L. lesueuri*. Therefore, the initial number of zoospores used in our model is based solely on these species and thus cannot account for any other frog species. Additionally, each frog species plays a unique role in the disease dynamics. For example, there are differences in zoospore shedding rate and resilience between such species (Lips, 2016).

- Studies have pointed out that moisture influences Bd survival (Lips, 2016; Raffel et al., 2015). We only built a TPC, which does not reflect other factors that affect Bd survival in the real world.

- Our model does not account for the epidemic spread of Bd that expends across different communities and even countries. This can occur naturally or with influence from human activities (i.e., trade) (Infection with *Batrachochytrium Dendrobatidis (Bd)*, 2020; Lips, 2016). 

## Remaining Plan 

The next step is to determine the values of the rate parameters in our model and to translate our deterministic model into a stochastic one. Then we can use our TPC to determine a series of initial zoospore counts (Z0) over a sequence of temperature values which we will use to solve our mathematical model. Based on our predictions, we expect to see that with rising temperatures due to climate change, there is a lower prevalence of *Bd* within a population of frogs thus reducing the probability of extinction caused by *Bd* for that population. 

# References 

Berger, L., Speare, R., Hines, H., Martelli, G., Hyatt, A., McDonald, K., Skerratt, L., Olsen, 
  V., Clarke, J., Gillespie, G., Mahony, M., Sheppard, N., Williams, C., & Tyler, M. (2004). 
  Effect Of Season And Temperature On Mortality In Amphibians Due To Chytridiomycosis. 
  *Australian Veterinary Journal, 82*(7), 434–439. 
  https://doi.org/10.1111/j.1751-0813.2004.tb11137.x
  
Bradley, P. W., Brawner, M. D., Raffel, T. R., Rohr, J. R., Olson, D. H., Blaustein, A. R., & 
  Kerby, J. (2019). Shifts in temperature influence how Batrachochytrium dendrobatidis infects 
  amphibian larvae. *PloS One, 14*(9), e0222237–e0222237. 
  https://doi.org/10.1371/journal.pone.0222237

*Infection with Batrachochytrium dendrobatidis (Bd)*. (2020). Australia Government, Department 
  of Agriculture, Water and the Environment. 
  https://www.agriculture.gov.au/sites/default/files/documents/infection-batrachochytrium-dendr
  obatidis.pdf

Kriger, K. M., & Hero, J. M. (2007). Large-scale seasonal variation in the prevalence and 
  severity of chytridiomycosis. *Journal of Zoology (1987), 271*(3), 352–359. 
  https://doi.org/10.1111/j.1469-7998.2006.00220.x

Kriger, K. M., Pereoglou, F., & Hero, J. M. (2007). Latitudinal Variation in the Prevalence and
  Intensity of Chytrid (Batrachochytrium dendrobatidis) Infection in Eastern Australia. 
  *Conservation Biology, 21*(5), 1280–1290. https://doi.org/10.1111/j.1523-1739.2007.00777.x
  
Lips, K. R. (2016). Overview of chytrid emergence and impacts on amphibians. *Philosophical* 
  *Transactions of the Royal Society of London. Series B. Biological Sciences, 371*(1709), 
  20150465-. https://doi.org/10.1098/rstb.2015.0465
  
Piotrowski, J. S., Annis, S. L., & Longcore, J. E. (2004). Physiology of Batrachochytrium 
  dendrobatidis, a chytrid pathogen of amphibians. *Mycologia, 96*(1), 9–15. 
  https://doi.org/10.1080/15572536.2005.11832990
  
Raffel, T. R., Halstead, N. T., McMahon, T. A., Davis, A. K., & Rohr, J. R. (2015). Temperature
  variability and moisture synergistically interact to exacerbate an epizootic disease. 
  *Proceedings of the Royal Society. B, Biological Sciences, 282*(1801), 20142039–20142039. 
  https://doi.org/10.1098/rspb.2014.2039
  
Rauw, W. M. (2012). Immune response from a resource allocation perspective. *Frontiers in 
  Genetics, 3*, 267–267. https://doi.org/10.3389/fgene.2012.00267
  
Retallick, R. W. R., McCallum, H., Speare, R., & Georgina M. Mace. (2004). Endemic infection of
  the amphibian chytrid fungus in a frog community post-decline. *PLoS Biology, 2*(11), 
  e351–e351. https://doi.org/10.1371/journal.pbio.0020351

Sapsford, S. J., Alford, R. A., & Schwarzkopf, L. (2013). Elevation, Temperature, and Aquatic 
  Connectivity All Influence the Infection Dynamics of the Amphibian Chytrid Fungus in Adult 
  Frogs. *PloS One, 8*(12), e82425–e82425. https://doi.org/10.1371/journal.pone.0082425

Scheele, B. C., Pasmans, F., Skerratt, L. F., Berger, L., Martel, A., Beukema, W., Acevedo, A. 
  A., Burrowes, P. A., Carvalho, T., Catenazzi, A., De la Riva, I., Fisher, M. C., Flechas, S. 
  V., Foster, C. N., Frías-Álvarez, P., Garner, T. W. J., Gratwicke, B., Guayasamin, J. M.,   
  Hirschfeld, M., … Canessa, S. (2019). Amphibian fungal panzootic causes catastrophic and    
  ongoing loss of biodiversity. *Science, 363*(6434), 1459–1463. 
  https://doi.org/10.1126/science.aav0379

Sheets, C. N., Schmidt, D. R., Hurtado, P. J., Byrne, A. Q., Rosenblum, E. B., 
  Richards-Zawacki, C. L., & Voyles, J. (2021). Thermal Performance Curves of Multiple Isolates
  of Batrachochytrium dendrobatidis, a Lethal Pathogen of Amphibians. *Frontiers in Veterinary*
  *Science, 8*, 687084–687084. https://doi.org/10.3389/fvets.2021.687084

Skerratt, L. F., Berger, L., Speare, R., Cashins, S., McDonald, K. R., Phillott, A. D., Hines, 
  H. B., & Kenyon, N. (2007). Spread of Chytridiomycosis Has Caused the Rapid Global Decline 
  and Extinction of Frogs. *EcoHealth, 4*(2), 125–134. 
  https://doi.org/10.1007/s10393-007-0093-5

Stevenson, L. A., Alford, R. A., Bell, S. C., Roznik, E. A., Berger, L., & Pike, D. A. (2013). 
  Variation in Thermal Performance of a Widespread Pathogen, the Amphibian Chytrid Fungus 
  *Batrachochytrium dendrobatidis*. *PLoS ONE, 8*(9), e73830. 
  https://doi.org/10.1371/journal.pone.0073830

Sun, D., Ellepola, G., Herath, J., & Meegaskumbura, M. (2023). The two chytrid pathogens of 
  amphibians in Eurasia—climatic niches and future expansion. *BMC Ecology and Evolution,* 
  *23*(1), 26–26. https://doi.org/10.1186/s12862-023-02132-y
  
Tolles, J., & Luong, T. (2020). Modeling Epidemics With Compartmental Models. *JAMA: The 
  Journal of the American Medical Association, 323*(24), 2515–2516. 
  https://doi.org/10.1001/jama.2020.8420

Tytar, V., Nekrasova, O., Pupins, M., Skute, A., Kirjušina, M., Gravele, E., Mezaraupe, L., 
  Marushchak, O., Čeirāns, A., Kozynenko, I., & Kulikova, A. A. (2023). Modeling the 
  Distribution of the Chytrid Fungus Batrachochytrium dendrobatidis with Special Reference to 
  Ukraine. *Journal of Fungi, 9*(6), 607-. https://doi.org/10.3390/jof9060607
  
Voyles, J., Johnson, L. R., Rohr, J., Kelly, R., Barron, C., Miller, D., Minster, J., & 
  Rosenblum, E. B. (2017). Diversity in growth patterns among strains of the lethal fungal 
  pathogen Batrachochytrium dendrobatidis across extended thermal optima. *Oecologia, 184*(2), 
  363–373. https://doi.org/10.1007/s00442-017-3866-8