diff --git a/R/server-functions.R b/R/server-functions.R
index 6e4a93d..9fb8e92 100644
--- a/R/server-functions.R
+++ b/R/server-functions.R
@@ -210,3 +210,29 @@ plot_epi <- function(model, mark_max) {
                    )
   return(plot)
 }
+
+#' plot_reproductive_epi Function
+#' 
+#' This function generates a plot of the reproductive number over time
+#' @param model The model object
+#' 
+#' @return A plot displaying the reproductive number for the model over the 
+#' course of the simulation
+#' 
+#' @export
+plot_reproductive_epi <- function (model) {
+  # Calculating average rep. number for each unique source_exposure_date
+  rep_num <- get_reproductive_number(model)
+  average_rt <- aggregate(rt ~ source_exposure_date, data = rep_num, 
+                          FUN = mean)
+  # Plotting
+  reproductive_plot <- plot_ly(data = average_rt, x = ~source_exposure_date, 
+                               y = ~rt, type = 'scatter', 
+                               mode = 'lines+markers')
+  reproductive_plot <- 
+    reproductive_plot |> 
+      plotly::layout(title = "Reproductive Number",
+                     xaxis = list(title = 'Day (step)'),
+                     yaxis = list(title = 'Average Rep. Number'))
+  return(reproductive_plot)         
+}
\ No newline at end of file
diff --git a/models/shiny_seir.R b/models/shiny_seir.R
index 4bfddcb..699762b 100644
--- a/models/shiny_seir.R
+++ b/models/shiny_seir.R
@@ -31,12 +31,7 @@ shiny_seir <- function(input) {
   # Summary
   summary_seir <- function() summary(model_seir)
   # Reproductive Number
-  reproductive_seir <- function() 
-    plot_reproductive_number(
-      model_seir,
-      main = "SEIR Model Reproductive Number"
-    )
-
+  reproductive_seir <- function() plot_reproductive_epi(model_seir)
   # Table 
   table_seir <- function() {
 
diff --git a/models/shiny_seirconn.R b/models/shiny_seirconn.R
index ee81fa4..6df440d 100644
--- a/models/shiny_seirconn.R
+++ b/models/shiny_seirconn.R
@@ -21,15 +21,8 @@ shiny_seirconn <- function(input) {
   
   # Plot, summary, and reproductive number
   plot_seirconn <- function() plot_epi(model_seirconn)
-  
   summary_seirconn <- function() summary(model_seirconn)
-  
-  reproductive_seirconn <- function()
-    plot_reproductive_number(
-      model_seirconn,
-      main = "SEIRCONNECTED Model Reproductive Number"
-      )
-
+  reproductive_seirconn <- function() plot_reproductive_epi(model_seirconn)
   # Table 
   table_seirconn <- function() {
     df <- as.data.frame(get_hist_total(model_seirconn))
diff --git a/models/shiny_seirconnequity.R b/models/shiny_seirconnequity.R
index 0a4b1d7..37fcd41 100644
--- a/models/shiny_seirconnequity.R
+++ b/models/shiny_seirconnequity.R
@@ -139,13 +139,9 @@ shiny_seirconnequity <- function(input) {
   }
   
   summary_seirconnequity <- function() summary(model_seirconnequity)
+  reproductive_seirconnequity <- 
+    function() plot_reproductive_epi(model_seirconnequity)
   
-  reproductive_seirconnequity <- function()
-    plot_reproductive_number(
-      model_seirconnequity,
-      main = "SEIRCONNECTED Model Reproductive Number"
-      )
-
   # Table 
   table_seirconnequity <- function() {
     df <- as.data.frame(get_hist_total(model_seirconnequity))
diff --git a/models/shiny_seird.R b/models/shiny_seird.R
index 8be8fa2..c49a262 100644
--- a/models/shiny_seird.R
+++ b/models/shiny_seird.R
@@ -32,45 +32,31 @@ shiny_seird <- function(input) {
   # Summary
   summary_seird <- function() summary(model_seird)
   # Reproductive Number
-  reproductive_seird <- function() 
-    plot_reproductive_number(
-      model_seird,
-      main = "SEIRD Model Reproductive Number"
-    )
-
+  reproductive_seird <- function() plot_reproductive_epi(model_seird)
   # Table 
   table_seird <- function() {
-
     df <- as.data.frame(get_hist_total(model_seird))
-    
     # Subset to only include "infection" state
     infection_data <- df[df$state == "Infected", ]
-    
     # Row with the maximum count
     max_infection_row <- infection_data[which.max(infection_data$counts), ]
-    
     # Row number of the maximum count in the original data frame
     max_row_number <- which(
       df$date == max_infection_row$date & df$state == "Infected"
       )
-
     df[max_row_number,"counts"] <- sprintf(
       "<strong>%s</strong>", 
       df[max_row_number, "counts"]
       )
-
     # Making sure factor variables are ordered
     df$state <- factor(
       x      = df$state,
       levels = c("Susceptible", "Exposed", "Infected", "Removed")
       )
-    
     # Reshaping the data to wide format
     df <- reshape(df, idvar = "date", timevar = "state", direction = "wide")
-
     colnames(df) <- gsub(colnames(df), pattern = "counts.", replacement = "")
     df
-    
   }
 
   # Output list
diff --git a/models/shiny_sir.R b/models/shiny_sir.R
index 3b84592..e76abbe 100644
--- a/models/shiny_sir.R
+++ b/models/shiny_sir.R
@@ -29,12 +29,7 @@ shiny_sir <- function(input) {
   # Plot, summary and repnum
   plot_sir <- function() plot_epi(model_sir)
   summary_sir <- function() summary(model_sir)
-  reproductive_sir <- function()
-    plot_reproductive_number(
-      model_sir,
-      main = "SIR Model Reproductive Number"
-    )
-
+  reproductive_sir <- function() plot_reproductive_epi(model_sir)
   # Table 
   table_sir <- function() {
     df <- as.data.frame(get_hist_total(model_sir))
diff --git a/models/shiny_sirconn.R b/models/shiny_sirconn.R
index 70cd361..c8806e3 100644
--- a/models/shiny_sirconn.R
+++ b/models/shiny_sirconn.R
@@ -19,15 +19,8 @@ shiny_sirconn <- function(input) {
   
   # Plot, summary, and reproductive number
   plot_sirconn <- function() plot_epi(model_sirconn)
-  
   summary_sirconn <- function() summary(model_sirconn)
-  
-  reproductive_sirconn <- function()
-    plot_reproductive_number(
-      model_sirconn,
-      main = "SIRCONNECTED Model Reproductive Number"
-      )
-
+  reproductive_sirconn <- function() plot_reproductive_epi(model_sirconn)
   # Table 
   table_sirconn <- function() {
     df <- as.data.frame(get_hist_total(model_sirconn))
diff --git a/models/shiny_sird.R b/models/shiny_sird.R
index 751b80f..0a2cba7 100644
--- a/models/shiny_sird.R
+++ b/models/shiny_sird.R
@@ -30,12 +30,7 @@ shiny_sird <- function(input) {
   # Plot, summary, and repnum
   plot_sird <- function() plot_epi(model_sird)
   summary_sird <- function() summary(model_sird)
-  reproductive_sird <- function()
-    plot_reproductive_number(
-      model_sird,
-      main = "SIRD Model Reproductive Number"
-    )
-
+  reproductive_sird <- function() plot_reproductive_epi(model_sird)
   # Table 
   table_sird <- function() {
     df <- as.data.frame(get_hist_total(model_sird))
diff --git a/models/shiny_sis.R b/models/shiny_sis.R
index 2af15cb..91157f2 100644
--- a/models/shiny_sis.R
+++ b/models/shiny_sis.R
@@ -28,15 +28,8 @@ shiny_sis <- function(input) {
   
   # Plot, summary, and reproductive number
   plot_sis <- function() plot_epi(model_sis)
-  
   summary_sis <- function() summary(model_sis)
-  
-  reproductive_sis <- function()
-    plot_reproductive_number(
-      model_sis,
-      main = "SIS Model Reproductive Number"
-      )
-
+  reproductive_sis <- function() plot_reproductive_epi(model_sis)
   # Table 
   table_sis <- function() {
     df <- as.data.frame(get_hist_total(model_sis))
diff --git a/models/shiny_sisd.R b/models/shiny_sisd.R
index a767530..73fcb2a 100644
--- a/models/shiny_sisd.R
+++ b/models/shiny_sisd.R
@@ -30,12 +30,7 @@ shiny_sisd <- function(input) {
   # Plot, summary, and repnum
   plot_sisd <- function() plot_epi(model_sisd)
   summary_sisd <- function() summary(model_sisd)
-  reproductive_sisd <- function()
-    plot_reproductive_number(
-      model_sisd,
-      main = "SISD Model Reproductive Number"
-    )
-
+  reproductive_sisd <- function() plot_reproductive_epi(model_sisd)
   # Table 
   table_sisd <- function() {
     df <- as.data.frame(get_hist_total(model_sisd))
diff --git a/server.R b/server.R
index 85faae5..655403d 100644
--- a/server.R
+++ b/server.R
@@ -45,7 +45,7 @@ function(input, output) {
   output$model_plot <- renderPlotly({
     model_output()$epicurves_plot()
   })
-  output$model_reproductive_plot <- renderPlot({
+  output$model_reproductive_plot <- renderPlotly({
     model_output()$reproductive_plot()
   })
   output$model_summary <- renderPrint({
diff --git a/ui.R b/ui.R
index 1c23e0f..1c00eec 100644
--- a/ui.R
+++ b/ui.R
@@ -99,7 +99,7 @@ body <- dashboardBody(
   ),
   fluidRow(
     column(6, plotlyOutput("model_plot") %>% withSpinner(color="#009bff")),
-    column(6, plotOutput("model_reproductive_plot") %>% withSpinner(color="#009bff"))
+    column(6, plotlyOutput("model_reproductive_plot") %>% withSpinner(color="#009bff"))
   ),
   HTML("<br>"),
   fluidRow(