Temperature- and density-dependent mosquito dynamics model with seasonal variation in carrying capacity

# ---- MODEL PARAMETERS (for calculation and dynamic content if needed)
sigma_base <- 0.025                # Low-density larval death (σ)
E_base     <- 1/19                 # Maturation rate (E)
delta_base <- 0.1                  # Adult death rate (Δ)
gamma_base <- 4                    # Fecundity (Γ, used for fitting to R0m if required)
xi_base    <- 8                    # Extrinsic incubation rate (ξ)
K_s        <- 0.3                  # Seasonal carrying capacity amplitude
alpha      <- 1.5                  # Adult wild-type density per person (α)
N          <- 1000                 # Human population (arbitrary)
kappa_base <- 0.6                  # Biting rate (κ)
beta_mh    <- 1.0                  # Prob. mosquito → human (Bmh)
eta        <- 4                    # Human infectious period (η)
R0         <- 2.69                 # Reproduction number (R0m for mosquito only, in this model)
days       <- seq(0, 730, by = 0.1)
T_seasonal <- 25 + 4 * cos(2 * pi * days / 365)  # Example seasonal temperature

# ---- Print and check parameters dynamically (optional)
cat("Model parameters (from paper):\n")

## Model parameters (from paper):

cat(sprintf("σ = %.3f, E = %.5f, Δ = %.3f, Γ = %.1f, ξ = %.5f, K_s = %.2f, α = %.1f, N = %d, κ = %.1f, Bmh = %.1f, η = %.1f, R0 = %.2f\n",
            sigma_base, E_base, delta_base, gamma_base, xi_base, K_s, alpha, N, kappa_base, beta_mh, eta, R0))

## σ = 0.025, E = 0.05263, Δ = 0.100, Γ = 4.0, ξ = 8.00000, K_s = 0.30, α = 1.5, N = 1000, κ = 0.6, Bmh = 1.0, η = 4.0, R0 = 2.69

# ---- Eq. 37: Mean carrying capacity (K_bar)
diagnostic_block <- function(mean_delta, mean_sigma, gamma_val, E_base, xi_constant, alpha) {
  term <- ((E_base * gamma_val - mean_delta) / (mean_delta * mean_sigma) - 1)
  if (term <= 0) stop(sprintf("K_bar inner term negative: %.6f", term))
  K_bar <- (alpha * mean_delta / E_base) * (term)^(-1/xi_constant)
  if (!is.finite(K_bar) || K_bar <= 0) stop(sprintf("K_bar not positive: %.6g", K_bar))
  return(K_bar)
}

# ---- Eq. 39: Probability of transmission from human to mosquito (Bhm)
beta_hm_calc <- function(R0, beta_mh, kappa, alpha, eta, delta) {
  denom <- (beta_mh / (kappa^2 * alpha * eta)) * ((1 + delta * xi_base) / delta)
  beta_hm <- R0 / denom
  return(beta_hm)
}

# ---- ODIN ODE Model: as per equations
mosquito_model_universal <- odin::odin({
  pi <- 3.141592653589793
  delta <- interpolate(time, delta_vals, "linear")
  sigma <- interpolate(time, sigma_vals, "linear")
  E <- interpolate(time, E_vals, "linear")
  gamma <- interpolate(time, gamma_vals, "linear")
  kappa <- interpolate(time, kappa_vals, "linear")
  xi <- interpolate(time, xi_vals, "linear")
  time_interpolated <- interpolate(time, time, "linear")
  K <- K_bar * (1 + K_s * cos(2 * pi * time_interpolated / 365))
  mu_density <- sigma * (1 + L / (K * N))
  lambda <- beta_hm * kappa * MI / N
  output(lambda) <- TRUE
  deriv(L)  <- gamma * MS - (E + mu_density) * L
  deriv(MS) <- E * L - (delta + lambda) * MS
  deriv(ME) <- lambda * MS - (xi + delta) * ME
  deriv(MI) <- xi * ME - delta * MI
  initial(L)  <- L0
  initial(MS) <- MS0
  initial(ME) <- ME0
  initial(MI) <- MI0
  time[]          <- user()
  T[]             <- user()
  delta_vals[]    <- user()
  sigma_vals[]    <- user()
  E_vals[]        <- user()
  gamma_vals[]    <- user()
  kappa_vals[]    <- user()
  xi_vals[]       <- user()
  dim(time) <- user()
  dim(T)    <- user()
  dim(delta_vals) <- user()
  dim(sigma_vals) <- user()
  dim(E_vals) <- user()
  dim(gamma_vals) <- user()
  dim(kappa_vals) <- user()
  dim(xi_vals) <- user()
  beta_hm <- user()
  K_bar   <- user()
  K_s     <- user()
  N       <- user()
  L0  <- user()
  MS0 <- user()
  ME0 <- user()
  MI0 <- user()
})

# ---- Helper for post-processing model output for plotting
scenario_output <- function(out, days, scenario, curve_type, compartment_names) {
  plot_indices <- seq(1, length(days), by = 10)
  if (ncol(out) < 5) stop("ODE output does not have expected columns! Actual columns: ", paste(colnames(out), collapse=", "))
  pop_data <- as.data.frame(out[plot_indices, 2:5])
  colnames(pop_data) <- compartment_names
  pop_data$Day <- days[plot_indices]
  pop_data$Scenario <- scenario
  pop_data$CurveType <- curve_type
  pop_long <- pop_data %>%
    pivot_longer(cols = compartment_names, names_to = "Compartment", values_to = "Population")
  pop_long
}

# ---- Baseline (temp-independent) parameter set for the model
mean_delta <- delta_base
mean_sigma <- sigma_base
gamma_constant <- gamma_base
K_bar_val <- diagnostic_block(mean_delta, mean_sigma, gamma_constant, E_base, xi_base, alpha)
beta_hm_val <- beta_hm_calc(R0, beta_mh, kappa_base, alpha, eta, mean_delta)

baseline_params <- list(
  time = days,
  T = T_seasonal,
  delta_vals = rep(mean_delta, length(days)),
  sigma_vals = rep(mean_sigma, length(days)),
  E_vals = rep(E_base, length(days)),
  gamma_vals = rep(gamma_constant, length(days)),
  kappa_vals = rep(kappa_base, length(days)),
  xi_vals = rep(xi_base, length(days)),
  beta_hm = beta_hm_val,
  K_bar = K_bar_val,
  K_s = K_s,
  N = N,
  L0  = K_bar_val * 10,
  MS0 = K_bar_val * 10,
  ME0 = K_bar_val * 0.25,
  MI0 = 0
)
cat(sprintf("K_bar (Eq. 37) = %.6g\n", K_bar_val))

## K_bar (Eq. 37) = 1.7799