#packages used
library('ggplot2')
library('deSolve')
library('patchwork')

Abstract

Acute paralysis virus infects honey bees colonies by using mites as a vector. This study addresses the disease dynamics by focusing on three critical parameters: transmission rate, healthy removal rate, and infectious removal rate. By isolating the healthy removal rate from the infectious removal rate, we are able to determine the relative effects of each parameter. My modelling study quantifies experiential apiary methods such as cold-snapping, pesticide-use, and population supplements. Further studies on the interplay between infectious and healthy removal rates may be required.

Introduction

Honey bees ( Apis Mellifera ) are one of nature’s greatest eusocial insects. Though individual bees tend and guard the hive, tight quarters in the honeycomb attract intruders. One such species is the Varroa destructor mite. The adult female Varroa mite feeds on honey bee haemolymph by piercing through and attaching itself to the bee’s intersegmental membrane (ISM). A parasitized bee has decreased fitness. Besides a reduction in fitness, the Varroa mite also acts as a vector for acute paralysis virus (APV). APV kills larvae and pupae before it can mature to an adult bee [1].

Transmission & Disease Information

Direct APV transmission between infected and non-infected bees is minimal [2]. APV transmission only occurs when a virus-carrying mite pierces the honey bee ISM. A virus-free mite can become infected if it pierces an infected bee. Moreover, if a parasitized bee dies, a virus-carrying mite can detach and infect another bee.

Characteristics of APV indicate that the mite vector is required for transmission. APV cannot proliferate without the vector.

Pathology

APV-carrying mites do not show any adverse reactions. This positive-sense RNA virus is 9,491 nucleotides long and contains two open reading frames [3]. Transcription of APV RNA causes nerve cell damage resulting in complete paralysis. Infected bees are unable to fly, lose hairs, and tremble rapidly.

Methods

Despite the complex nature of mite, bee, and virus dynamics, this model can be simplified into three basic parameters: β is the infection rate of un-infected bees, α is the removal rate of infected bees, and γ is the removal rate of healthy bees.

Sources for each parameter were derived from multiple studies, with both removal rates (α and γ) taken from Fukuda & Sekiguchi’s study on seasonal change and life span [4]. The study was conducted over eight months from April (spring) to November (wintering). Every ten days, two honey bee ( Apis Mellifera ) colonies were examined by removing and tagging 200-300 young worker bees. Tagged bees were summed for either infected or healthy and the colonies’ removal rate was calculated.

The infection rate was determined using Bailey’s 1964 study on APV viruses in Austria and Switzerland. In this study, 20 bees from two separate colonies were extracted and placed in a cage. The bees were placed into different temperatures to simulate seasons. Each bee was injected with APV viral sample to determine the probability of infection with exposure to virus [5].

Parameters for removal rate included ranges for removal counts in different seasons; however, infection rate was determined by changing incubation temperature to simulate infection in different seasons. Thus, the infection rate may not be fully representative of field conditions. Moreover, the injection of identical RNA viral doses may not accurately portray stochasticity in mite-to-bee RNA exchange.

Below are parameter values extracted:

## 
## 
## |value name |  spring| summer|  autumn|   winter|
## |:----------|-------:|------:|-------:|--------:|
## |beta       | 0.19840|  0.146| 0.19000| 0.033840|
## |alpha      | 0.20000|  0.200| 0.20000| 0.005300|
## |gamma      | 0.02272|  0.040| 0.02272| 0.005263|

Results

Disease Dynamics

The below uses ODEs, rather than the discrete Euler approximations.

Spring dynamics base + plot code, includes raw data 


#spring dynamics

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values <- c(
  beta = values$spring[1],
  alpha = values$spring[2],
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values
)

simulation_spring <- as.data.frame(simulation_spring)
simulation_spring
##    time             S        I        R
## 1     1  9.990000e+02   1.0000   0.0000
## 2     2  3.009746e-10 806.5782 193.4218
## 3     3 -2.031009e-11 645.5357 354.4643
## 4     4  4.167399e-10 516.6471 483.3529
## 5     5 -1.235364e-09 413.4929 586.5071
## 6     6 -1.231309e-09 330.9347 669.0653
## 7     7  1.540449e-09 264.8600 735.1400
## 8     8 -1.808350e-08 211.9780 788.0220
## 9     9  7.628691e-08 169.6542 830.3458
## 10   10 -2.652082e-07 135.7810 864.2190
#plot spring

p_1 <- ggplot(data = simulation_spring, aes(x = time)) +
  geom_line(aes(y = simulation_spring$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Spring Bee SIR Model') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
Summer dynamics base + plot code, includes raw data 


#summer dynamics

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

summer_values <- c(
  beta = values$summer[1],
  alpha = values$summer[2],
  gamma = values$summer[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_summer <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- summer_values
)

simulation_summer <- as.data.frame(simulation_summer)
simulation_summer
##    time             S        I        R
## 1     1  9.990000e+02   1.0000   0.0000
## 2     2  4.205067e-09 795.6297 204.3703
## 3     3 -7.643415e-09 625.8643 374.1357
## 4     4  5.501336e-09 492.3221 507.6779
## 5     5 -8.597457e-09 387.2744 612.7256
## 6     6 -1.079606e-07 304.6409 695.3591
## 7     7  1.499479e-07 239.6392 760.3608
## 8     8 -3.007435e-07 188.5071 811.4929
## 9     9  1.635670e-06 148.2850 851.7150
## 10   10 -8.390663e-06 116.6453 883.3547
#plot summer

p_2 <- ggplot(data = simulation_summer, aes(x = time)) +
  geom_line(aes(y = simulation_summer$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_summer$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_summer$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Summer Bee SIR Model') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
Autumn dynamics base + plot code, includes raw data 


#autumn dynamics

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}
  
autumn_values <- c(
  beta = values$autumn[1],
  alpha = values$autumn[2],
  gamma = values$autumn[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_autumn <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- autumn_values
)

simulation_autumn <- as.data.frame(simulation_autumn)
simulation_autumn
##    time             S        I        R
## 1     1  9.990000e+02   1.0000   0.0000
## 2     2  2.945765e-10 806.8556 193.1444
## 3     3 -1.179727e-10 645.7577 354.2423
## 4     4  3.319920e-10 516.8248 483.1752
## 5     5 -1.194587e-09 413.6351 586.3649
## 6     6 -1.144043e-09 331.0485 668.9515
## 7     7  1.228912e-09 264.9510 735.0490
## 8     8 -1.721266e-08 212.0508 787.9492
## 9     9  7.144888e-08 169.7126 830.2874
## 10   10 -2.393565e-07 135.8277 864.1723
#plot autumn

p_3 <- ggplot(data = simulation_autumn, aes(x = time)) +
  geom_line(aes(y = simulation_autumn$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_autumn$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_autumn$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Autumn Bee SIR Model') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
Winter dynamics base + plot code, includes raw data 


#winter dynamics

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

winter_values <- c(
  beta = values$winter[1],
  alpha = values$winter[2],
  gamma = values$winter[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)
t <- 1:10

simulation_winter <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- winter_values
)

simulation_winter <- as.data.frame(simulation_winter)
simulation_winter
##    time             S        I         R
## 1     1  9.990000e+02   1.0000  0.000000
## 2     2 -3.194731e-08 991.6293  8.370713
## 3     3 -2.207709e-11 981.2098 18.790183
## 4     4 -6.653963e-11 970.8998 29.100153
## 5     5  1.129948e-11 960.6982 39.301790
## 6     6 -6.638201e-12 950.6038 49.396236
## 7     7  2.530708e-13 940.6154 59.384615
## 8     8  6.502164e-13 930.7320 69.268043
## 9     9  7.078220e-13 920.9524 79.047622
## 10   10 -3.097723e-13 911.2756 88.724444
#plot winter

p_4 <- ggplot(data = simulation_winter, aes(x = time)) +
  geom_line(aes(y = simulation_winter$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_winter$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_winter$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Winter Bee SIR Model') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
#simulation plots

all_plots <- (p_1 + p_2) / (p_3 + p_4)

all_plots + plot_annotation(
  title = 'HoneyBee SIR Simulations from Spring to Winter') +
  plot_layout(guides = 'collect')

Caption: SIR simulations for Honey Bee Populations. Susceptible counts decreased too rapidly to be shown accurately on the winter plots. Simulations were only run for 10 days with a starting bee population of 1,000 (999 susceptible and 1 infected). Y-axis was capped at 1,000 individuals.

Using only the simplified model parameters, we can see that infected bee numbers have a linear and proportional relationship to susceptible bee counts. This relationship is arrested at around day 2 when the infected numbers peak. Spring, summer, and autumn infected counts peaked at around 800 bees; however, winter infected counts peaked at 992 bees. This winter result was expected as the rate of honey bee removal (α and γ) was very low. Since these values are extracted from a more complex model, parameters such as “mite (vector) birth rate”, “carrying capacity of mites given bee population”, and “tolerance of bees to infection” do not come into play.

The initial size of the population does not affect the outcome.

Parameter Uncertainties

#uncertainty values 1 (u1 - halved beta)

c1_u1 <- c("beta", "alpha", "gamma")
c2_u1 <- c((0.1984/2), 0.2, 0.02272)
c3_u1 <- c("NA", "NA", "NA")
c4_u1 <- c("NA", "NA", "NA")
c5_u1 <- c((0.03384/2), 0.0053, 0.005263)
values_u1 <- data.frame(c1_u1, c2_u1, c3_u1, c4_u1, c5_u1)
names(values_u1) <- c("value name", "spring", "summer", "autumn", "winter")
print(knitr::kable(values_u1))
## 
## 
## |value name |  spring|summer |autumn |   winter|
## |:----------|-------:|:------|:------|--------:|
## |beta       | 0.09920|NA     |NA     | 0.016920|
## |alpha      | 0.20000|NA     |NA     | 0.005300|
## |gamma      | 0.02272|NA     |NA     | 0.005263|
#uncertainty values 2 (u2 - tripled gamma)

c1_u2 <- c("beta", "alpha", "gamma")
c2_u2 <- c(0.1984, 0.2, (0.02272*3))
c3_u2 <- c("NA", "NA", "NA")
c4_u2 <- c("NA", "NA", "NA")
c5_u2 <- c(0.03384, 0.0053, (0.005263*3))
values_u2 <- data.frame(c1_u2, c2_u2, c3_u2, c4_u2, c5_u2)
names(values_u2) <- c("value name", "spring", "summer", "autumn", "winter")
print(knitr::kable(values_u2))
## 
## 
## |value name |  spring|summer |autumn |   winter|
## |:----------|-------:|:------|:------|--------:|
## |beta       | 0.19840|NA     |NA     | 0.033840|
## |alpha      | 0.20000|NA     |NA     | 0.005300|
## |gamma      | 0.06816|NA     |NA     | 0.015789|
##uncertainty values 3 (u3 - tripled alpha) - not included in plot analysis

c1_u3 <- c("beta", "alpha", "gamma")
c2_u3 <- c(0.1984, (0.2*3), 0.02272)
c3_u3 <- c("NA", "NA", "NA")
c4_u3 <- c("NA", "NA", "NA")
c5_u3 <- c(0.03384, (0.0053*3), 0.005263)
values_u3 <- data.frame(c1_u3, c2_u3, c3_u3, c4_u3, c5_u3)
names(values_u3) <- c("value name", "spring", "summer", "autumn", "winter")
Simulation and plot code for uncertainties 


model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_u1 <- c(
  beta = values_u1$spring[1],
  alpha = values_u1$spring[2],
  gamma = values_u1$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_u1 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_u1
)

simulation_spring_u1 <- as.data.frame(simulation_spring_u1)
p_1_u1 <- ggplot(data = simulation_spring_u1, aes(x = time)) +
  geom_line(aes(y = simulation_spring_u1$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring_u1$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring_u1$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Spring Uncertainties SIR', subtitle = 'Transmission Rate Halved') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

winter_values_u1 <- c(
  beta = values_u1$winter[1],
  alpha = values_u1$winter[2],
  gamma = values_u1$winter[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_winter_u1 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- winter_values_u1
)

simulation_winter_u1 <- as.data.frame(simulation_winter_u1)
p_4_u1 <- ggplot(data = simulation_winter_u1, aes(x = time)) +
  geom_line(aes(y = simulation_winter_u1$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_winter_u1$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_winter_u1$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Winter Uncertainties SIR', subtitle = 'Transmission Rate Halved') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_u2 <- c(
  beta =  values_u2$spring[1],
  alpha = values_u2$spring[2],
  gamma = values_u2$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_u2 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_u2
)

simulation_spring_u2 <- as.data.frame(simulation_spring_u2)
p_1_u2 <- ggplot(data = simulation_spring_u2, aes(x = time)) +
  geom_line(aes(y = simulation_spring_u2$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring_u2$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring_u2$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Spring Uncertainties SIR', subtitle = 'Healthy Removal Rate Tripled') +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

winter_values_u2 <- c(
  beta = values_u2$winter[1],
  alpha = values_u2$winter[2],
  gamma = values_u2$winter[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_winter_u2 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- winter_values_u2
)

simulation_winter_u2 <- as.data.frame(simulation_winter_u2)
p_4_u2 <- ggplot(data = simulation_winter_u2, aes(x = time)) +
  geom_line(aes(y = simulation_winter_u2$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_winter_u2$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_winter_u2$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 1000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Winter Uncertainties SIR', subtitle = 'Healthy Removal Rate Tripled') +
  scale_color_identity(name = 'Legend',
                       breaks = c( 'red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')
#uncertainty plots. original plot codes hidden, see above drop-down for full source

all_plots <- (p_1_u1 + p_4_u1) / (p_1_u2 + p_4_u2)

all_plots + plot_annotation(
  title = 'Honey Bee Uncertainty Simulations Spring & Winter') +
  plot_layout(guides = 'collect')

Caption: Plots showing uncertainty for spring and winter simulations. These simulations were selected because they contained the most different simulations numbers. Summer and autumn simulations were quite similar to spring. Two conditions were tested: halved transmission rate and tripled healthy removal rate.

By halving transmission rate and then tripling the healthy removal rate, we introduced uncertainties to the parameters. Comparing uncertainty simulation outcomes to the base simulations, there is little difference in the end stage.

Since little changes to the simulation outcome were found, it is clear that direct viral transmission does not play a key role in the transmission of APV in honey bees. Without other parameters mentioned above, transmission dynamics cannot be determined.

Moreover, in all simulations above, honey bee removal outputs combined healthy (γ) and infected (α) removal rates. We are unable to tell if the bees die due to recovery or infection. I wish to add more complexity to determine the difference between healthy removal rate and infected removal rate. In the “control” section, I will separately test the magnitude of healthy and infected removal rates.

Host Density Threshold for Persistence

Threshold values code 


#threshold value code

NT_spring <- (values$spring[2] +  values$spring[3]) / values$spring[1]
NT_summer <- (values$summer[2] +  values$summer[3]) / values$summer[1]
NT_autumn <- (values$autumn[2] +  values$autumn[3]) / values$autumn[1]
NT_winter <- (values$winter[2] +  values$winter[3]) / values$winter[1]

NT_spring_u1 <- (values_u1$spring[2] +  values_u1$spring[3]) / values_u1$spring[1]
NT_winter_u1 <- (values_u1$winter[2] +  values_u1$winter[3]) / values_u1$winter[1]

NT_spring_u2 <- (values_u2$spring[2] +  values_u2$spring[3]) / values_u2$spring[1]
NT_winter_u2 <- (values_u2$winter[2] +  values_u2$winter[3]) / values_u2$winter[1]

NT <- data.frame(
  seasons = c("spring", "summer", "autumn", "winter"),
  base_sim = c(NT_spring, NT_summer, NT_autumn, NT_winter),
  u1_sim = c(NT_spring_u1, "NA", "NA", NT_winter_u1),
  u2_sim = c(NT_spring_u2, "NA", "NA", NT_winter_u2))

Below are the threshold values for spring to winter. Summer and autumn uncertainties were not tested.

## 
## 
## |seasons |  base_sim|u1_sim            |u2_sim            |
## |:-------|---------:|:-----------------|:-----------------|
## |spring  | 1.1225806|2.24516129032258  |1.35161290322581  |
## |summer  | 1.6438356|NA                |NA                |
## |autumn  | 1.1722105|NA                |NA                |
## |winter  | 0.3121454|0.624290780141844 |0.623197399527187 |
Threshold u1 winter simulation code 


#threshold u1 spring simulation


model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_u1_thres <- c(
  beta = values_u1$spring[1],
  alpha = values_u1$spring[2],
  gamma = values_u1$spring[3]
)

initial_model_values_thres <- c(
  S = 1,
  I = 1,
  R = 0
)

t <- 1:25

simulation_spring_u1_thres <- ode(
  y <- initial_model_values_thres,
  times <- t,
  func <- model,
  parms <- spring_values_u1_thres
)

simulation_spring_u1_thres <- as.data.frame(simulation_spring_u1_thres)
Threshold u1 winter plot code 


#u1 spring threshold plot

p_1_u1_thres <- ggplot(data = simulation_spring_u1_thres, aes(x = time)) +
  geom_line(aes(y = simulation_spring_u1_thres$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring_u1_thres$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring_u1_thres$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 2)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = 'Spring U1 Threshold SIR', subtitle = 'Transmission rate halved; Initial susceptible number set below threshold value') +
  scale_color_identity(name = 'Legend',
                       breaks = c( 'red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')

Caption: Spring simulation with transmission rate halved and honey bee population at 2 (one infected and one susceptible). We can see that the disease was eradicated after 20 days: an equilibrium between removed and susceptible counts was reached. As time approaches infinity, infectious counts approach zero. Moreover, susceptible counts plateaued at 1.5 from below. Whether that signifies healthy or infected removed, is unknown.

Given that transmission spread does not depend on the honey bee population size, the host density threshold makes sense. If there was a way to incorporate mite population size, this threshold would be more informative. Knowing the threshold values, I ran a simple threshold check with spring uncertainty parameters of halved transmission and host (honey bee) population of less than the threshold (<2.245). Simulation with these parameters yielded population recovery from disease.

One point to note is the threshold values for winter. Given the slower infection rate, I expected winter threshold to be higher than spring NT values; however, this was not the case. Rather, the ratio of removal rate to transmission rate was lower for winter. This difference is likely due to the infected removal rate (α) being much lower than the other seasons.

Healthy and Infected Interplay

This section is dedicated to determining the interplay between healthy removal rate (γ) and infected removal rate (α). I take isolated each value in both spring and winter conditions to see what effect each removal rate has on overall population numbers.

Spring healthy isolated code 


#isolated spring healthy - ish

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_ish <- c(
  beta = values$spring[1],
  alpha = 0,
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_ish <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_ish
)

simulation_spring_ish <- as.data.frame(simulation_spring_ish)
simulation_spring_ish
##    time             S        I         R
## 1     1  9.990000e+02   1.0000   0.00000
## 2     2 -7.561997e-11 978.3098  21.69015
## 3     3  4.249141e-14 956.3333  43.66675
## 4     4 -8.333141e-11 934.8503  65.14966
## 5     5 -3.178309e-11 913.8500  86.14999
## 6     6  1.404765e-10 893.3214 106.67859
## 7     7  6.430825e-11 873.2540 126.74602
## 8     8 -7.174328e-11 853.6373 146.36267
## 9     9  3.483648e-11 834.4614 165.53864
## 10   10  9.266143e-11 815.7161 184.28385
Spring infected isolated code 


#isolated spring infected - isi

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_isi <- c(
  beta = values$spring[1],
  alpha = values$spring[2],
  gamma = 0
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_isi <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_isi
)

simulation_spring_isi <- as.data.frame(simulation_spring_isi)
simulation_spring_isi
##    time             S        I        R
## 1     1  9.990000e+02   1.0000   0.0000
## 2     2  6.591155e-09 824.4590 175.5410
## 3     3  2.139497e-10 675.0099 324.9901
## 4     4 -1.617062e-09 552.6512 447.3488
## 5     5 -2.436153e-09 452.4728 547.5272
## 6     6 -6.506169e-08 370.4534 629.5466
## 7     7  4.312902e-08 303.3018 696.6982
## 8     8 -3.056763e-07 248.3227 751.6773
## 9     9  1.174338e-06 203.3094 796.6906
## 10   10 -2.174048e-06 166.4558 833.5442
Winter healthy isolated code 


#isolated winter healthy - iwh

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

winter_values_iwh <- c(
  beta = values$winter[1],
  alpha = values$winter[2],
  gamma = values$winter[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_winter_iwh <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- winter_values_iwh
)

simulation_winter_iwh <- as.data.frame(simulation_winter_iwh)
simulation_winter_iwh
##    time             S        I         R
## 1     1  9.990000e+02   1.0000  0.000000
## 2     2 -3.194731e-08 991.6293  8.370713
## 3     3 -2.207709e-11 981.2098 18.790183
## 4     4 -6.653963e-11 970.8998 29.100153
## 5     5  1.129948e-11 960.6982 39.301790
## 6     6 -6.638201e-12 950.6038 49.396236
## 7     7  2.530708e-13 940.6154 59.384615
## 8     8  6.502164e-13 930.7320 69.268043
## 9     9  7.078220e-13 920.9524 79.047622
## 10   10 -3.097723e-13 911.2756 88.724444
Winter infected isolated code 


#isolated winter healthy - iwi

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

winter_values_iwi <- c(
  beta = values$winter[1],
  alpha = values$winter[2],
  gamma = 0
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_winter_iwi <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- winter_values_iwi
)

simulation_winter_iwi <- as.data.frame(simulation_winter_iwi)
simulation_winter_iwi
##    time             S        I         R
## 1     1  9.990000e+02   1.0000  0.000000
## 2     2 -3.175463e-08 995.7910  4.209019
## 3     3  1.086919e-11 990.5272  9.472754
## 4     4 -7.794198e-10 985.2913 14.708662
## 5     5  1.276028e-10 980.0831 19.916890
## 6     6 -2.436265e-11 974.9024 25.097589
## 7     7 -4.422557e-12 969.7491 30.250903
## 8     8  4.217208e-12 964.6230 35.376977
## 9     9  3.458022e-13 959.5240 40.475955
## 10   10 -5.140981e-13 954.4520 45.547979

Caption: Comparison of base plots with isolated healthy (γ) and infected removal rate (α) simulations.

By comparing spring simulations in panels I and III, we can see that γ has comparatively little effect on the population end-stage as α. This is to say that the effect of α is magnitudes greater than γ in the spring. For plots II and IV, isolation of either variable yields no significant difference from the base plot; although, the γ simulation seems to have a higher removal increase. This result indicates infectious removal rate (α) has a higher effect in the spring and both removal rates are equally weighted in the winter.

Controlling Spread

Since this population has a low threshold value, culling of any bee populations does not produce any benefit. Thus, methods to control APV spread lie exclusively with external environmental controls and mite population controls.

Control Option 1

The first control option involves temperature control. As illustrated in the winter simulations, winter colonies have a lower transmission rate (β). If it is possible to simulate winter-like conditions without reducing healthy removal rate, then the temperature may work to slow down APV spread.

#1 Temperature control code main 


#1 temperature control

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c1 <- c(
  beta = 0.005,
  alpha = values$spring[2],
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_c1 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c1
)

simulation_spring_c1 <- as.data.frame(simulation_spring_c1)
simulation_spring_c1
##    time            S        I          R
## 1     1 9.990000e+02   1.0000   0.000000
## 2     2 8.896223e+02 105.2125   5.165233
## 3     3 7.259226e+01 810.6180 116.789708
## 4     4 1.492672e+00 708.6956 289.811714
## 5     5 6.198170e-02 568.4103 431.527700
## 6     6 4.849381e-03 454.9698 545.025393
## 7     7 6.310486e-04 364.1336 635.865733
## 8     8 1.233262e-04 291.4308 708.569122
## 9     9 3.339493e-05 233.2434 766.756531
## 10   10 1.173150e-05 186.6738 813.326199
#1 Temperature control code - gamma only 


#1 temperature control gamma only

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c1_g <- c(
  beta = 0.005,
  alpha = 0,
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_c1_g <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c1_g
)

simulation_spring_c1_g <- as.data.frame(simulation_spring_c1_g)
simulation_spring_c1_g
##    time             S        I           R
## 1     1  9.990000e+02   1.0000   0.0000000
## 2     2  8.727230e+02 126.6629   0.6140615
## 3     3  4.590574e+01 940.0980  13.9962786
## 4     4  3.674414e-01 963.6989  35.9336971
## 5     5  3.129703e-03 942.4082  57.5886383
## 6     6  2.970599e-05 921.2412  78.7588058
## 7     7  3.108021e-07 900.5466  99.4534292
## 8     8 -4.299554e-08 880.3168 119.6831619
## 9     9 -8.374236e-07 860.5415 139.4584616
## 10   10 -2.595145e-07 841.2105 158.7895250
#1 Temperature control code - alpha only 


#1 temperature control alpha only

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c1_a <- c(
  beta = 0.005,
  alpha = values$spring[2],
  gamma = 0
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_c1_a <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c1_a
)

simulation_spring_c1_a <- as.data.frame(simulation_spring_c1_a)
simulation_spring_c1_a
##    time            S        I          R
## 1     1 9.990000e+02   1.0000   0.000000
## 2     2 8.878171e+02 107.4633   4.719567
## 3     3 6.913449e+01 824.0375 106.827994
## 4     4 1.296093e+00 732.8078 265.896130
## 5     5 4.661370e-02 601.0487 398.904687
## 6     6 3.058647e-03 492.1348 507.862112
## 7     7 3.288421e-04 402.9283 597.071370
## 8     8 5.294856e-05 329.8900 670.109898
## 9     9 1.189711e-05 270.0915 729.908529
## 10   10 3.498988e-06 221.1326 778.867421
(p_1_c1 +  guide_area() + plot_layout(guides = 'collect') + p_1_c1_g + p_1_c1_a) + plot_annotation(
  title = 'Spring Control 1 SIR Model', tag_levels = 'I', subtitle = bquote(~beta~"is 0.005 for all plots"))

Caption: Control 1 SIR simulations done with infection rate (β) at 0.005. This means the infection rate per bee per day is 0.5%. Plot I shows population numbers with combined γ and α. Plot II and III show population numbers with isolated γ and α to see which removal rate has greater effect. At least 3/4 of the removed population is infected, whereas 1/4 is healthy.

Here we assumed a temperature decrease lowers the infection rate (β) and does not affect removal rate of both healthy and infected. We can see that in a population of 999 susceptible and 1 infected, the infection spread within colony is slower than the base plot. However, this is only a short-term solution as removed bees still dominate the end-stage (with 3/4 dominance). Healthy-removed bee counts are statistically higher than the base model, but still remain low. Thus, although the easiest to implement, control option 1 is not the best long-term solution.

Control Option 2

Next, it would be beneficial to see the interaction between an elevated removal rate and the rest of the SIR system. By treating the bees with pesticides, mite counts in honey bee colonies may decrease. In this simple SIR model, the effects of pesticide treatment are likely to influence the number of healthy-removed bees. Thus, increasing healthy bee removal rate (γ) can be a proxy for effective pesticide treatment.

#2 Healthy control code 


#2 healthy rate control


model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c2 <- c(
  beta = values$spring[1],
  alpha = values$spring[2],
  gamma = 0.8
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_c2 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c2
)

simulation_spring_c2 <- as.data.frame(simulation_spring_c2)
simulation_spring_c2
##    time            S           I        R
## 1     1 9.990000e+02   1.0000000   0.0000
## 2     2 1.935216e-09 380.9968743 619.0031
## 3     3 2.126486e-08 140.1609375 859.8391
## 4     4 7.269967e-06  51.5625641 948.4374
## 5     5 3.531056e-08  18.9688070 981.0312
## 6     6 3.348738e-09   6.9782436 993.0218
## 7     7 1.375074e-09   2.5671563 997.4328
## 8     8 9.961003e-10   0.9444031 999.0556
## 9     9 8.848233e-10   0.3474258 999.6526
## 10   10 8.470967e-10   0.1278106 999.8722
#2 Healthy control code - gamma only 


#2 healthy control gamma only

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c2_g <- c(
  beta = values$spring[1],
  alpha = 0,
  gamma = 0.8
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_c2_g <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c2_g
)

simulation_spring_c2_g <- as.data.frame(simulation_spring_c2_g)
simulation_spring_c2_g
##    time             S           I        R
## 1     1  9.990000e+02   1.0000000   0.0000
## 2     2 -1.374256e-10 462.0847891 537.9152
## 3     3  1.274545e-09 207.6284595 792.3715
## 4     4 -2.246713e-07  93.2939242 906.7061
## 5     5 -3.456050e-06  41.9198464 958.0802
## 6     6  5.455879e-08  18.8357855 981.1642
## 7     7  3.721657e-09   8.4634718 991.5365
## 8     8  4.511859e-11   3.8028822 996.1971
## 9     9  1.680552e-11   1.7087437 998.2913
## 10   10  5.877575e-12   0.7677871 999.2322
#2 Healthy control code - alpha only 


#2 healthy control alpha only

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c2_a <- c(
  beta = values$spring[1],
  alpha = values$spring[2],
  gamma = 0
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:10

simulation_spring_c2_a <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c2_a
)

simulation_spring_c2_a <- as.data.frame(simulation_spring_c2_a)
simulation_spring_c2_a
##    time             S        I        R
## 1     1  9.990000e+02   1.0000   0.0000
## 2     2  6.591155e-09 824.4590 175.5410
## 3     3  2.139497e-10 675.0099 324.9901
## 4     4 -1.617062e-09 552.6512 447.3488
## 5     5 -2.436153e-09 452.4728 547.5272
## 6     6 -6.506169e-08 370.4534 629.5466
## 7     7  4.312902e-08 303.3018 696.6982
## 8     8 -3.056763e-07 248.3227 751.6773
## 9     9  1.174338e-06 203.3094 796.6906
## 10   10 -2.174048e-06 166.4558 833.5442
(p_1_c2 +  guide_area() + plot_layout(guides = 'collect') + p_1_c2_g + p_1_c2_a) + plot_annotation(
  title = 'Spring Control 2 SIR Model', tag_levels = 'I', subtitle = bquote(~gamma~"is 0.8 for all plots, except for plot III where"~gamma~"= 0"))

Caption: Control 2 SIR simulations done with healthy removal rate (γ) is 0.8. This means that the healthy removal rate per bee per day is 80%. Plot I shows population numbers with combined removal rates (γ and α). Plots II and III show population numbers with isolated γ and α. Plot II has high correlation with the control 2 simulation, indicating γ is the greatest contributor to population levels.

Here, we see that increasing γ values as a proxy for a pesticide treatment killing mites allowed for successful healthy bee population recovery. Infected populations only reached a maximum of 500 bees, as seen in plot II. Of note is the interaction between γ and α values. In plot III, infectious removal populations are still high. This discrepancy may indicate some type of competition between both removal parameters. Thus, with these simplistic SIR ODEs, it is impossible to pinpoint the competitive effect between both γ and α parameters. In practice, this form of competition may be pathological: with 80% chance of survival and 20% chance of infected death.

Pesticide treatments may be effective in killing mites, but may be harmful to local insect biodiversity [4]. Mites have also evolved pesticide-resistant strains.

Control Option 3

Even though synthetic pesticides may fail, honey bee colonies have acquired a natural defense: bigger is better. An increased population size leads to stronger colonies and lower infection [6]. In this final control method, I supplement existing simulations with more susceptible bees at intervals of 2, 4, 6, and 8. At each addition, infectious removal rate (α) and infection rate (β) is reduced by 3 times (multiplied by 1/3).

#3 Population supplement code 


#3 population supplements section 1

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c3_sec1 <- c(
  beta = values$spring[1],
  alpha = values$spring[2],
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = 999,
  I = 1,
  R = 0
)

t <- 1:2

simulation_spring_c3_sec1 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c3_sec1
)

simulation_spring_c3_sec1 <- as.data.frame(simulation_spring_c3_sec1)
simulation_spring_c3_sec1
##   time            S        I        R
## 1    1 9.990000e+02   1.0000   0.0000
## 2    2 3.009746e-10 806.5782 193.4218
#3 population supplements section 2

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c3_sec2 <- c(
  beta = (values$spring[1]/3),
  alpha = (values$spring[2]/3),
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = (3.009746e-10+999),
  I = 806.5782,
  R = 193.4218
)

t <- 1:2

simulation_spring_c3_sec2 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c3_sec2
)

simulation_spring_c3_sec2 <- as.data.frame(simulation_spring_c3_sec2)
simulation_spring_c3_sec2
##   time             S         I        R
## 1    1  9.990000e+02  806.5782 193.4218
## 2    2 -3.751952e-11 1652.1841 346.8159
#3 population supplements section 3

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c3_sec3 <- c(
  beta = ((values$spring[1]/3)/3),
  alpha = ((values$spring[2]/3)/3),
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = (-3.751952e-11+999),
  I = 1652.1841 ,
  R = 346.8159
)

t <- 1:2

simulation_spring_c3_sec3 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c3_sec3
)

simulation_spring_c3_sec3 <- as.data.frame(simulation_spring_c3_sec3)
simulation_spring_c3_sec3
##   time             S        I        R
## 1    1  9.990000e+02 1652.184 346.8159
## 2    2 -2.026017e-11 2535.594 462.4056
#3 population supplements section 4

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c3_sec4 <- c(
  beta = (((values$spring[1]/3)/3)/3),
  alpha = (((values$spring[2]/3)/3)/3),
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = (-2.026017e-11+999),
  I = 2535.594,
  R = 462.4056
)

t <- 1:16

simulation_spring_c3_sec4 <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c3_sec4
)

simulation_spring_c3_sec4 <- as.data.frame(simulation_spring_c3_sec4)
simulation_spring_c3_sec4
##    time             S        I         R
## 1     1  9.990000e+02 2535.594  462.4056
## 2     2 -6.620656e-09 3431.019  565.9810
## 3     3 -4.239525e-10 3329.192  667.8072
## 4     4  7.642100e-11 3230.388  766.6112
## 5     5 -1.626751e-10 3134.517  862.4830
## 6     6 -9.173561e-10 3041.490  955.5093
## 7     7 -3.662721e-10 2951.225 1045.7750
## 8     8  1.276861e-09 2863.638 1133.3619
## 9     9  8.480513e-10 2778.650 1218.3493
## 10   10 -7.206897e-10 2696.185 1300.8144
## 11   11 -3.754103e-11 2616.168 1380.8320
## 12   12  1.313975e-09 2538.525 1458.4749
## 13   13 -3.175579e-10 2463.186 1533.8135
## 14   14 -2.167234e-09 2390.083 1606.9162
## 15   15  7.616497e-10 2319.150 1677.8494
## 16   16  3.937006e-09 2250.322 1746.6773
#3 population supplements section 4 gamma only

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c3_sec4_g <- c(
  beta = (((values$spring[1]/3)/3)/3),
  alpha = 0,
  gamma = values$spring[3]
)

initial_model_values <- c(
  S = (-2.026017e-11+999),
  I = 2535.594,
  R = 462.4056
)

t <- 1:16

simulation_spring_c3_sec4_g <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c3_sec4_g
)

simulation_spring_c3_sec4_g <- as.data.frame(simulation_spring_c3_sec4_g)
simulation_spring_c3_sec4_g
##    time             S        I         R
## 1     1  9.990000e+02 2535.594  462.4056
## 2     2 -8.583033e-09 3456.199  540.8004
## 3     3 -2.067960e-10 3378.559  618.4405
## 4     4 -3.852282e-10 3302.663  694.3363
## 5     5  7.402297e-11 3228.473  768.5266
## 6     6 -9.655607e-12 3155.949  841.0504
## 7     7  1.062724e-11 3085.055  911.9451
## 8     8 -3.626664e-12 3015.752  981.2472
## 9     9 -8.593835e-12 2948.007 1048.9925
## 10   10 -1.586068e-11 2881.783 1115.2161
## 11   11 -9.766868e-12 2817.048 1179.9521
## 12   12 -4.309888e-12 2753.766 1243.2338
## 13   13  1.437491e-12 2691.906 1305.0940
## 14   14  4.597058e-13 2631.435 1365.5646
## 15   15 -1.057171e-12 2572.323 1424.6767
## 16   16  1.060250e-12 2514.539 1482.4610
#3 population supplements section 4 alpha only

model <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- -beta * I * S
    dI <- (beta * S - (alpha + gamma)) * I
    dR <- (alpha + gamma) * I
    return(list(c(dS, dI, dR)))
  })
}

spring_values_c3_sec4_a <- c(
  beta = (((values$spring[1]/3)/3)/3),
  alpha = (((values$spring[2]/3)/3)/3),
  gamma = 0
)

initial_model_values <- c(
  S = (-2.026017e-11+999),
  I = 2535.594,
  R = 462.4056
)

t <- 1:16

simulation_spring_c3_sec4_a <- ode(
  y <- initial_model_values,
  times <- t,
  func <- model,
  parms <- spring_values_c3_sec4_a
)

simulation_spring_c3_sec4_a <- as.data.frame(simulation_spring_c3_sec4_a)
simulation_spring_c3_sec4_a
##    time             S        I        R
## 1     1  9.990000e+02 2535.594 462.4056
## 2     2 -1.239763e-08 3508.841 488.1585
## 3     3  9.770069e-10 3482.946 514.0539
## 4     4 -1.765102e-09 3457.241 539.7583
## 5     5 -8.875659e-10 3431.727 565.2728
## 6     6  2.751459e-10 3406.401 590.5991
## 7     7  1.365386e-10 3381.261 615.7384
## 8     8 -7.704743e-11 3356.307 640.6923
## 9     9  6.703784e-14 3331.538 665.4619
## 10   10  3.118957e-11 3306.951 690.0488
## 11   11 -9.155411e-12 3282.545 714.4542
## 12   12 -1.015392e-11 3258.320 738.6796
## 13   13  6.949838e-12 3234.274 762.7261
## 14   14  2.689648e-12 3210.404 786.5952
## 15   15 -4.240148e-12 3186.712 810.2881
## 16   16  8.972769e-13 3163.193 833.8061
simulation_spring_c3_joined <- rbind(simulation_spring_c3_sec1, simulation_spring_c3_sec2, simulation_spring_c3_sec3, simulation_spring_c3_sec4)

simulation_spring_c3_joined$time <- (1:22)
simulation_spring_c3_joined
##    time             S         I         R
## 1     1  9.990000e+02    1.0000    0.0000
## 2     2  3.009746e-10  806.5782  193.4218
## 3     3  9.990000e+02  806.5782  193.4218
## 4     4 -3.751952e-11 1652.1841  346.8159
## 5     5  9.990000e+02 1652.1841  346.8159
## 6     6 -2.026017e-11 2535.5944  462.4056
## 7     7  9.990000e+02 2535.5940  462.4056
## 8     8 -6.620656e-09 3431.0186  565.9810
## 9     9 -4.239525e-10 3329.1924  667.8072
## 10   10  7.642100e-11 3230.3884  766.6112
## 11   11 -1.626751e-10 3134.5166  862.4830
## 12   12 -9.173561e-10 3041.4903  955.5093
## 13   13 -3.662721e-10 2951.2246 1045.7750
## 14   14  1.276861e-09 2863.6377 1133.3619
## 15   15  8.480513e-10 2778.6503 1218.3493
## 16   16 -7.206897e-10 2696.1852 1300.8144
## 17   17 -3.754103e-11 2616.1676 1380.8320
## 18   18  1.313975e-09 2538.5247 1458.4749
## 19   19 -3.175579e-10 2463.1861 1533.8135
## 20   20 -2.167234e-09 2390.0834 1606.9162
## 21   21  7.616497e-10 2319.1502 1677.8494
## 22   22  3.937006e-09 2250.3223 1746.6773
simulation_spring_c3_joined_a <- rbind(simulation_spring_c3_sec1, simulation_spring_c3_sec2, simulation_spring_c3_sec3, simulation_spring_c3_sec4_a)

simulation_spring_c3_joined_a$time <- (1:22)
simulation_spring_c3_joined_a
##    time             S         I        R
## 1     1  9.990000e+02    1.0000   0.0000
## 2     2  3.009746e-10  806.5782 193.4218
## 3     3  9.990000e+02  806.5782 193.4218
## 4     4 -3.751952e-11 1652.1841 346.8159
## 5     5  9.990000e+02 1652.1841 346.8159
## 6     6 -2.026017e-11 2535.5944 462.4056
## 7     7  9.990000e+02 2535.5940 462.4056
## 8     8 -1.239763e-08 3508.8411 488.1585
## 9     9  9.770069e-10 3482.9457 514.0539
## 10   10 -1.765102e-09 3457.2413 539.7583
## 11   11 -8.875659e-10 3431.7268 565.2728
## 12   12  2.751459e-10 3406.4005 590.5991
## 13   13  1.365386e-10 3381.2612 615.7384
## 14   14 -7.704743e-11 3356.3073 640.6923
## 15   15  6.703784e-14 3331.5377 665.4619
## 16   16  3.118957e-11 3306.9508 690.0488
## 17   17 -9.155411e-12 3282.5454 714.4542
## 18   18 -1.015392e-11 3258.3200 738.6796
## 19   19  6.949838e-12 3234.2735 762.7261
## 20   20  2.689648e-12 3210.4044 786.5952
## 21   21 -4.240148e-12 3186.7115 810.2881
## 22   22  8.972769e-13 3163.1935 833.8061
simulation_spring_c3_joined_g <- rbind(simulation_spring_c3_sec1, simulation_spring_c3_sec2, simulation_spring_c3_sec3, simulation_spring_c3_sec4_g)

simulation_spring_c3_joined_g$time <- (1:22)
simulation_spring_c3_joined_g
##    time             S         I         R
## 1     1  9.990000e+02    1.0000    0.0000
## 2     2  3.009746e-10  806.5782  193.4218
## 3     3  9.990000e+02  806.5782  193.4218
## 4     4 -3.751952e-11 1652.1841  346.8159
## 5     5  9.990000e+02 1652.1841  346.8159
## 6     6 -2.026017e-11 2535.5944  462.4056
## 7     7  9.990000e+02 2535.5940  462.4056
## 8     8 -8.583033e-09 3456.1992  540.8004
## 9     9 -2.067960e-10 3378.5591  618.4405
## 10   10 -3.852282e-10 3302.6633  694.3363
## 11   11  7.402297e-11 3228.4730  768.5266
## 12   12 -9.655607e-12 3155.9492  841.0504
## 13   13  1.062724e-11 3085.0545  911.9451
## 14   14 -3.626664e-12 3015.7524  981.2472
## 15   15 -8.593835e-12 2948.0071 1048.9925
## 16   16 -1.586068e-11 2881.7835 1115.2161
## 17   17 -9.766868e-12 2817.0475 1179.9521
## 18   18 -4.309888e-12 2753.7658 1243.2338
## 19   19  1.437491e-12 2691.9056 1305.0940
## 20   20  4.597058e-13 2631.4350 1365.5646
## 21   21 -1.057171e-12 2572.3229 1424.6767
## 22   22  1.060250e-12 2514.5386 1482.4610
p_1_c3 <- ggplot(data = simulation_spring_c3_joined, aes(x = time)) +
  geom_line(aes(y = simulation_spring_c3_joined$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring_c3_joined$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring_c3_joined$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 4000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = bquote("Both"~alpha~"and"~gamma~"are included")) +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')

p_1_c3_a <- ggplot(data = simulation_spring_c3_joined_a, aes(x = time)) +
  geom_line(aes(y = simulation_spring_c3_joined_a$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring_c3_joined_a$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring_c3_joined_a$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 4000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = bquote("Only"~alpha~"is included")) +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')

p_1_c3_g <- ggplot(data = simulation_spring_c3_joined_g, aes(x = time)) +
  geom_line(aes(y = simulation_spring_c3_joined_g$S, color = 'red'), size = 1.5) +
  geom_line(aes(y = simulation_spring_c3_joined_g$I, color = 'blue'), size = 1.5) +
  geom_line(aes(y = simulation_spring_c3_joined_g$R, color = 'green'), size = 1.5) +
  scale_y_continuous(limits = c(0, 4000)) +
  theme_classic() +
  labs(x = 'Time (Days)', y = 'Number of Bees', title = bquote("Only"~gamma~"is included")) +
  scale_color_identity(name = 'Legend',
                       breaks = c('red', 'blue', 'green'),
                       labels = c('Susceptible', 'Infectious', 'Removed'),
                       guide = 'legend')

Caption: Control 3 SIR simulations done by starting with base simulations until timestep 2. Next, 999 susceptible bees are added to the system. At the meantime, α and β are reduced by 3 times. This process is repeated until timestep 8, after which the simulation was left to run for another 16 days. Plots I and II show the same process, but with γ and α isolated.

In control 3 simulations, we can see that healthy removal rates formulate a significant proportion at the end-stage. This result assumes that the supplement of additional population counts lowers the transmission rate (β) and the infectious removal rate (α) proportionally. This assumption is based on current research determining higher honey bee population counts results in a more robust colony. Thus, with these assumptions, I have portrayed population supplements in infected colonies as a viable long-term solution to mite infestations and APV disease.

Discussion

Descriptions of expected dynamics in each control option are found below their respective plot. I do not wish to repeat them here. However, I will comment on the effectiveness of each control method.

The first control method involved temperature. Temperatures only lower the transmission rate (β). This type of cold-shocking only slows down mite-transmitted APV. Although not effective in the long run, this method can be used effectively as damage control before other solutions are found.

My second control method involved pesticides. Pesticides increase the healthy removal rate (γ) by acting as a proxy for successful mite reduction. This method has the most short-term effectiveness as most removed bees are projected to be healthy.

My third control method involves population supplements. This method is based on the assumption that larger colony sizes lead to lower transmission and infectious death. This method is the most successful in the long-term, with a balance of healthy and infectious removed populations shown on the plots. The mechanisms of this process are unknown and deserve further study (perhaps with more complex models).

In apiary programs, all three methods have been used (perhaps even all in tandem). The effectiveness of each solution seems to correspond well with my simple SIR simulations, with cold-shock slowing APV spread, pesticides having short-term benefit, and population supplements having the greatest long-term benefit.

In my own experience in beekeeping, mites have always been found inside hives. Low mite populations do not pose any risk to robust colonies; only after a long wintering period (stress) do mite counts rise and treatment is required. Interpreting the winter base plot, we conclude that winter conditions lead to the slowest transmission, but has the lowest host density threshold. This means that if transmission rates were to rise, colonies in the winter would suffer the quickest infectious death. This begs the question: how long can a cold-snap be before the effectiveness of slowing transmission causes undue stress on the honey bee colony?

Through our simple SIR simulations with actual honey bee APV transmission, infectious removal, and healthy removal rates, we are able to correlate expected populations changes with real-life inferences. Further work on the interplay between infectious and health removal rates is necessary to produce accurate results.

References

Please find the references here: References on uToronto Sharepoint (Click here)

All Code

Please find entire R Markdown document here: Entire R Markdown doc