# Michaelis-Menten fits by nonlinear least squares, the Lineweaver-Burk alternative, and two ways to compare Km.
# R 4.3.3, base R (stats::nls).
d <- read.csv("kinetics.csv")                        # columns: experiment, enzyme, substrate_uM, well, rate

fit <- function(x) {                                 # nonlinear least squares: estimates, SEs, residual SS, df
  f <- nls(rate ~ vmax * substrate_uM / (km + substrate_uM), data = x,
           start = list(vmax = max(x$rate), km = median(x$substrate_uM)))
  cf <- summary(f)$coefficients
  list(p = cf[, 1], se = cf[, 2], ss = sum(resid(f)^2), df = nrow(x) - 2)
}
lineweaver_burk <- function(x) {                     # straight line through 1/rate against 1/substrate
  b <- coef(lm(I(1 / rate) ~ I(1 / substrate_uM), data = x))
  c(vmax = 1 / b[[1]], km = b[[2]] / b[[1]])
}

x <- subset(d, experiment == 1 & enzyme == "WT")
r <- fit(x); t <- qt(0.975, r$df); l <- lineweaver_burk(x)
cat(sprintf("Experiment 1, WT, nonlinear regression: Km %.2f uM (95%% CI %.2f to %.2f); Vmax %.1f (95%% CI %.1f to %.1f)\n",
            r$p[["km"]], r$p[["km"]] - t * r$se[["km"]], r$p[["km"]] + t * r$se[["km"]],
            r$p[["vmax"]], r$p[["vmax"]] - t * r$se[["vmax"]], r$p[["vmax"]] + t * r$se[["vmax"]]))
cat(sprintf("Experiment 1, WT, Lineweaver-Burk: Km %.2f uM; Vmax %.1f\n", l[["km"]], l[["vmax"]]))

cat("Km (uM) in each experiment, nonlinear regression [Lineweaver-Burk]:\n")
lr <- c()
for (e in sort(unique(d$experiment))) {
  k <- list()
  for (g in c("WT", "mutant")) {
    x <- subset(d, experiment == e & enzyme == g)
    k[[g]] <- fit(x)$p[["km"]]; k[[paste0(g, "_lb")]] <- lineweaver_burk(x)[["km"]]
  }
  lr <- c(lr, log(k$mutant / k$WT))
  cat(sprintf("  Experiment %d: WT %.2f [%.2f]; mutant %.2f [%.2f]\n", e, k$WT, k$WT_lb, k$mutant, k$mutant_lb))
}
tt <- t.test(lr); n <- length(lr)
cat(sprintf("Km ratio, mutant / WT, one ratio per experiment (paired t-test on log Km, n = %d): %.2f (95%% CI %.2f to %.2f), p = %.3f\n",
            n, exp(mean(lr)), exp(tt$conf.int[1]), exp(tt$conf.int[2]), tt$p.value))

# Global fit: all wells together, a separate Vmax for each preparation, one Km per enzyme (or one shared Km)
d$grp <- (d$experiment - 1) * 2 + (d$enzyme == "mutant") + 1; d$mut <- (d$enzyme == "mutant") + 1; ng <- max(d$grp)
full <- nls(rate ~ vmax[grp] * substrate_uM / (km[mut] + substrate_uM), data = d, start = list(vmax = rep(100, ng), km = c(10, 10)))
red <- nls(rate ~ vmax[grp] * substrate_uM / (km + substrate_uM), data = d, start = list(vmax = rep(100, ng), km = 10))
ss_full <- sum(resid(full)^2); ss_red <- sum(resid(red)^2); df_full <- nrow(d) - ng - 2
F <- (ss_red - ss_full) / (ss_full / df_full); pF <- pf(F, 1, df_full, lower.tail = FALSE)
cf <- summary(full)$coefficients; t <- qt(0.975, df_full)
kw <- cf["km1", 1]; km_ <- cf["km2", 1]; sw <- cf["km1", 2]; sm <- cf["km2", 2]
cat(sprintf("Global fit of all %d wells (a Vmax for each preparation, Km shared across experiments): Km WT %.1f (95%% CI %.1f to %.1f), mutant %.1f (95%% CI %.1f to %.1f); extra-sum-of-squares F(1, %d) = %.2f, p = %.1e\n",
            nrow(d), kw, kw - t * sw, kw + t * sw, km_, km_ - t * sm, km_ + t * sm, df_full, F, pF))