# Illustration of hook effect
# Thanks to Sarrita Adams, rexvlucyletby2023.com ("Insulin Science").
#
# The "capture antibody" coats the wells of a tray
# We add a liquid in which the substance whose concentration we want
# to measure, in this case C-peptide, is dissolved, as well as
# "detection antibodies" and "signal antibodies".
# I can leave the signal antibodies out of the story, since
# the important action concerns capture, C-peptide, and detection particles.
# I will give them the names A, B and C.
# The idea is that A and B particles bind to form what I will call AB particles,
# then C particles bind to them forming ABC particles. At the end of the reactions
# we are going to measure the quantity of ABC and the idea is that it should be
# proportional to the original concentration of B.
#
# However there is another binding which competes for B, namely B and C can
# bind to form BC.
#
# What happens depends on the initial concentrations of A, B and C, and on
# the rates of the three reactions. I'll denote the rates by a.b, ab.c and b.c.
#
# At the end of the reactions we have final concentrations of six substances
# which I denote A, AB, ABC, BC, and C
#
# I use the R library "deSolve", and within it, the function "ode".
# In other words, we will solve a system of six ordinary differential
# equations.
#
# I will do this for a range of different values of the initial
# concentration of the "C-peptide" B. Notice how the final concentration
# of ABC initially increases linearly with the C-peptide concentration,
# but then turns around, and at very high concentrations becomes negligible.
# The "b.c" reaction won from the "a.b" reaction and as a result hardly any ABC
# was generated


library(deSolve)
nrep <- 21
concentrations <- seq(0, 2000, 100)
results <- rep(0, 21)
u0 <- c(A = 1000, AB = 0, ABC = 0, B = 80, BC = 0, C = 1000)
p0 <-c(a.b = 1e-04 , ab.c = 1e-05, b.c = 1e-03)
times <- seq(0, 500, by = 0.01)
update <- function(t, u, p){##
  a.b  <- p["a.b"]##
  ab.c <- p["ab.c"]##
  b.c  <- p["b.c"]##
  A   <- u["A"]##
  AB  <- u["AB"]##
  ABC <- u["ABC"]##
  B   <- u["B"]##
  BC  <- u["BC"]##
  C   <- u["C"]##
  dA <-  -a.b * A * B##
  dAB <-  a.b * A * B  -ab.c * AB * C##
  dABC <- ab.c * AB * C##
  dB <- -a.b * A * B -b.c * B * C##
  dBC <- b.c * B * C##
  dC <- - ab.c * AB * C - b.c * B * C##
  list(c(dA, dAB, dABC, dB, dBC, dC))
}
for (i in 1:nrep){#
 B0 <- concentrations[i]#
 u0["B"] <- B0#
 out <- ode(y  = u0, times = times, func = update, parms = p0)#
 results[i] <- out[50001, "ABC"]#
 }
plot(concentrations, results, type = "b")