---
title: "ParathyroidHormone (Ahn 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ParathyroidHormone (Ahn 2014)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(nlmixr2lib)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r readModelDb("Ahn_2014_parathyroidHormone")()$meta$reference`
* Description: `r readModelDb("Ahn_2014_parathyroidHormone")()$meta$description`
* Article: <https://doi.org/10.4196/kjpp.2014.18.3.217> (open access, Korean Journal of Physiology & Pharmacology)

This is a semi-mechanistic indirect-response PD model. Calcium intake (oral) is
absorbed into an unobserved Ca pool that inhibits parathyroid hormone (PTH)
secretion via an Emax model with Emax fixed to 1. A parallel homeostatic
indirect-response model describes observed plasma ionized Ca with no coupling
to the calcium dose -- the paper's central insight is that the dose does not
visibly perturb observed Ca (because of Ca homeostasis), but the PTH response
to the absorbed Ca is a sensitive marker of calcium absorption.

The model is K-PD-style: the depot and the unobserved Ca pool carry no
physiological mass scale; the depot is normalized to 1 unit at t=0 (Ahn 2014
Methods p. 218), and the relative bioavailability of CaCO3 vs thermal spring
water (Ahn 2014 Table 2 'Relative F1' = 1.98) acts as a multiplicative shift on
the depot bioavailability. Validation therefore follows the
`endogenous-validation.md` recipe (steady-state, perturbation-recovery,
mass-balance, replicate the published VPCs).

# Population

24 healthy Korean adults (22 male / 2 female; age 21-39 years, median 26;
weight 55.1-79.3 kg, mean 68.5) from a single randomized parallel trial
conducted at Seoul St. Mary's Hospital (Catholic University of Korea). Subjects
received calcium as either Geumjin thermal spring water (240 mL containing 400
mg elemental Ca; n=12) or calcium carbonate tablets (500 mg CaCO3 = 2 x 200 mg
elemental Ca, with 240 mL normal saline or 340 mL purified water; n=12, with
the two sub-arms pooled in the final analysis). Days 1 and 7 received the full
morning dose; Days 2-6 split the daily dose into twice-daily administrations.
Blood samples for ionized Ca and PTH were collected pre-dose and at 0.5, 1,
1.5, 2, 3, 4, 6, and 8 hours post-dose on Days 1 and 7. Demographics are in
Ahn 2014 Table 1.

```{r metadata}
mod_fn <- readModelDb("Ahn_2014_parathyroidHormone")
mod    <- mod_fn()
str(mod$meta$population)
str(mod$meta$covariateData)
```

# Source trace

Per-parameter origin is recorded as in-file comments next to each `ini()` entry
in `inst/modeldb/endogenous/Ahn_2014_parathyroidHormone.R`. The table below
collects them in one place for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka = log(0.796)` | ka = 0.796 1/hr | Ahn 2014 Table 2 (95% CI 0.438-1.15) |
| `lkin_ca = log(3.41)` | kin_ca = 3.41 mmol/L/hr | Ahn 2014 Table 2 (95% CI 3.39-3.43) |
| `lkout_ca = fixed(log(2.86))` | kout_ca = 2.86 1/hr (FIXED) | Ahn 2014 Table 2 + Results p. 220 (fixed to Abraham et al. 2009 literature value 2.862 1/hr) |
| `lkin_pth = log(21.6)` | kin_pth = 21.6 pg/mL/hr | Ahn 2014 Table 2 (95% CI 12.4-30.8) |
| `lkout_pth = log(0.849)` | kout_pth = 0.849 1/hr | Ahn 2014 Table 2 (95% CI 0.513-1.19) |
| `lec50 = log(0.158)` | EC50 = 0.158 mmol/L | Ahn 2014 Table 2 (95% CI 0.0924-0.224) |
| `emax = fixed(1)` | Emax = 1 (FIXED) | Ahn 2014 Methods p. 218 + Results p. 220 |
| `lfdepot = fixed(log(1))` | F_thermal = 1 (FIXED reference) | Ahn 2014 Methods p. 218 |
| `e_form_caco3_fdepot = log(1.98)` | Relative F1 = 1.98 for CaCO3 | Ahn 2014 Table 2 (RSE 24%; 95% CI 1.06-2.90) |
| `etalka ~ 0.406` | IIV variance on ka (~71% CV) | Ahn 2014 Table 2 (95% CI 0.161-0.651) |
| `etalkin_ca ~ 0.000229` | IIV variance on kin_ca (~1.5% CV) | Ahn 2014 Table 2 (published CI '0.000814, 0.000377' appears order-inverted; bootstrap CI 0.0001-0.0004 confirms magnitude) |
| `etalkin_pth ~ 0.0453` | IIV variance on kin_pth (~21.4% CV) | Ahn 2014 Table 2 (Results p. 220 narrative 'random individual difference of 21.4%') |
| `addSd = 0.0258` | Additive SD on observed Ca, mmol/L | Ahn 2014 Table 2 (SD_ca; 95% CI 0.0212-0.0304) |
| `propSd_PTH = 0.21` | Proportional SD on PTH, fraction | Ahn 2014 Table 2 (CV_pth = 21.0%; 95% CI 18.8-23.2%) |
| `d/dt(depot) = -ka * depot` | n/a | Ahn 2014 Methods p. 218 (schematic Fig. 1) |
| `d/dt(ca_unobs) = ka*depot - kout_ca*ca_unobs` | n/a | Ahn 2014 Methods p. 218 (eq. for unobserved Ca) |
| `d/dt(ca) = kin_ca - kout_ca*ca` | n/a | Ahn 2014 Methods p. 218 (eq. for observed Ca) |
| `d/dt(pth) = kin_pth*(1 - Emax*ca_unobs/(EC50+ca_unobs)) - kout_pth*pth` | n/a | Ahn 2014 Methods p. 218 (eq. for PTH with Emax inhibition) |
| `f(depot) = exp(lfdepot + e_form_caco3_fdepot * FORM_CACO3)` | F=1 (thermal) or 1.98 (CaCO3) | Ahn 2014 Methods p. 218 + Table 2 |
| Steady-state IC `ca_unobs(0) = 0` | 0 | Implicit: no calcium dose has been absorbed at t=0 (Ahn 2014 Methods p. 218 narrative: '[Ca2+_unobs] starts from 0 before dose'); also Fig. 1 |
| Steady-state IC `ca(0) = kin_ca/kout_ca` | ~1.19 mmol/L | Ahn 2014 Methods p. 218 assumption (no circadian rhythm) -- drug-free steady state of observed Ca |
| Steady-state IC `pth(0) = kin_pth/kout_pth` | ~25.4 pg/mL | Ahn 2014 Methods p. 218 assumption (no circadian rhythm) -- drug-free steady state of PTH at t=0 when ca_unobs=0 |

## Units (per ODE term)

The packaged model is K-PD-style: the depot amount and the unobserved Ca pool
carry no physiological mass scale. ka transfers the normalized depot signal
into the unobserved Ca pool, which is then matched dimensionally against EC50
(mmol/L) for the Emax inhibition. The paper treats the unobserved Ca pool
implicitly in mmol/L (so it pairs with EC50 in the Emax denominator) without
specifying a strict mass-balance scaling for the depot -> unobserved Ca
transition.

| ODE term | Units (paper's bookkeeping) | Notes |
|----------|-----------------------------|-------|
| `d/dt(depot)` | unit/hr | Normalized depot amount (dose at t=0 set to 1 with bioavailability F multiplying); `ka * depot` is the transfer rate. |
| `d/dt(ca_unobs)` | mmol/L/hr | `ka * depot` is interpreted as a concentration-rate input (paper's K-PD convention); `kout_ca * ca_unobs` is mmol/L/hr. |
| `d/dt(ca)` | mmol/L/hr | `kin_ca` has units mmol/L/hr; `kout_ca * ca` is mmol/L/hr. Independent of dose. |
| `d/dt(pth)` | pg/mL/hr | `kin_pth` has units pg/mL/hr; the bracketed inhibition factor is dimensionless; `kout_pth * pth` is pg/mL/hr. |
| `inhib = emax * ca_unobs / (ec50 + ca_unobs)` | dimensionless | Both ca_unobs and ec50 in mmol/L. |

# Parameter table (paper vs. file)

```{r param-table}
data.frame(
  parameter   = c("ka (1/hr)", "kin_ca (mmol/L/hr)", "kout_ca (1/hr, FIXED)",
                  "kin_pth (pg/mL/hr)", "kout_pth (1/hr)", "EC50 (mmol/L)",
                  "Emax (FIXED)", "Relative F1 (CaCO3 vs thermal)",
                  "IIV(ka) variance", "IIV(kin_ca) variance",
                  "IIV(kin_pth) variance",
                  "Additive SD on Ca (mmol/L)", "Proportional SD on PTH (CV)"),
  paper       = c("0.796 (95% CI 0.438-1.15)",
                  "3.41 (95% CI 3.39-3.43)",
                  "2.86 (literature, Abraham 2009)",
                  "21.6 (95% CI 12.4-30.8)",
                  "0.849 (95% CI 0.513-1.19)",
                  "0.158 (95% CI 0.0924-0.224)",
                  "1 (set during model development)",
                  "1.98 (RSE 24%; 95% CI 1.06-2.90)",
                  "0.406 (95% CI 0.161-0.651)",
                  "0.000229 (CI in paper appears typo'd; bootstrap 0.0001-0.0004)",
                  "0.0453 (95% CI 0.0260-0.0646; 21.4% CV)",
                  "0.0258 (95% CI 0.0212-0.0304)",
                  "21.0% (95% CI 18.8-23.2%)"),
  packaged    = c(round(exp(-0.2281561), 3),
                  round(exp(1.22671229), 3),
                  round(exp(1.05082162), 3),
                  round(exp(3.07269331), 3),
                  round(exp(-0.16369609), 3),
                  round(exp(-1.84516024), 3),
                  1,
                  round(exp(0.68309684), 3),
                  0.406, 0.000229, 0.0453, 0.0258, 0.21)
)
```

# Steady-state check

Solve the model without any calcium intake (no dose), zero IIV, for 24 hours.
The state should remain at the seeded baseline to numerical tolerance:
ca = kin_ca / kout_ca = 3.41 / 2.86 ~= 1.192 mmol/L, and pth = kin_pth /
kout_pth = 21.6 / 0.849 ~= 25.44 pg/mL.

```{r steady-state}
mod_typ <- mod |> rxode2::zeroRe()

ev_ss <- data.frame(
  id         = 1L,
  time       = seq(0, 24, by = 0.5),
  amt        = 0,
  evid       = 0L,
  cmt        = "Cc",
  FORM_CACO3 = 0
)

sim_ss <- rxode2::rxSolve(mod_typ, events = ev_ss, returnType = "data.frame")
range(sim_ss$Cc)
range(sim_ss$PTH)
```

```{r steady-state-assert}
stopifnot(diff(range(sim_ss$Cc))  < 1e-6)
stopifnot(diff(range(sim_ss$PTH)) < 1e-6)
stopifnot(abs(sim_ss$Cc[1]  - 3.41 / 2.86)   < 1e-3)
stopifnot(abs(sim_ss$PTH[1] - 21.6 / 0.849)  < 1e-2)
```

# Typical-value PTH and Ca trajectories after a single dose

Simulate the typical-value (zero-IIV) response to a single calcium dose under
each treatment arm. Per the paper's K-PD convention the depot is dosed with
`amt = 1` at t=0; bioavailability `F = 1` for thermal spring water and
`F = 1.98` for CaCO3.

```{r typical-dose}
make_typical_events <- function(form_caco3, times = seq(0, 8, by = 0.1)) {
  rbind(
    data.frame(id = 1L, time = 0, amt = 1, evid = 1L, cmt = "depot",
               FORM_CACO3 = form_caco3),
    data.frame(id = 1L, time = times, amt = 0, evid = 0L, cmt = "Cc",
               FORM_CACO3 = form_caco3)
  )
}

sim_typical <- bind_rows(
  rxode2::rxSolve(mod_typ, events = make_typical_events(0), returnType = "data.frame") |>
    mutate(treatment = "Thermal spring water (F=1)"),
  rxode2::rxSolve(mod_typ, events = make_typical_events(1), returnType = "data.frame") |>
    mutate(treatment = "Calcium carbonate (F=1.98)")
)
```

```{r typical-pth-figure}
ggplot(sim_typical, aes(time, PTH, colour = treatment)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 21.6 / 0.849, linetype = "dashed", colour = "grey50") +
  labs(x = "Time (hr)", y = "PTH (pg/mL)",
       title = "Typical-value PTH response after a single calcium dose",
       caption = "Dashed line: drug-free baseline PTH = kin_pth/kout_pth = 25.4 pg/mL. Compare to Ahn 2014 Figure 5.") +
  theme_minimal()
```

```{r typical-ca-figure}
ggplot(sim_typical, aes(time, Cc, colour = treatment)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 3.41 / 2.86, linetype = "dashed", colour = "grey50") +
  labs(x = "Time (hr)", y = "Observed Ca (mmol/L)",
       title = "Typical-value observed Ca trajectory (no coupling to dose by design)",
       caption = "Per the paper's K-PD design, observed Ca is decoupled from the calcium dose and stays at the homeostatic baseline ~1.19 mmol/L. Reproduces Ahn 2014 Figure 4 (no notable dose effect).") +
  theme_minimal()
```

The typical-value PTH curve shows the expected indirect-response shape: a
nadir at ~1-2 hours post-dose (when ca_unobs peaks), followed by recovery
toward baseline as ca_unobs is eliminated. The CaCO3 arm shows a deeper nadir
because its higher bioavailability drives a larger unobserved-Ca exposure.

# Visual predictive check (replicates Ahn 2014 Figures 4 and 5)

Reproduce the VPCs from Ahn 2014 Figures 4 (ionized Ca) and 5 (PTH) at the
"per treatment" stratification (lower panels of those figures). 250 virtual
subjects per treatment arm at the Day-1 sampling design.

```{r vpc-cohort}
set.seed(20140214L)
n_per_arm <- 250L

obs_times <- c(0, 0.5, 1, 1.5, 2, 3, 4, 6, 8)

make_cohort <- function(n, form_caco3, id_offset = 0L) {
  ids <- id_offset + seq_len(n)
  dose_rows <- data.frame(id = ids, time = 0, amt = 1, evid = 1L,
                          cmt = "depot", FORM_CACO3 = form_caco3)
  obs_rows  <- expand.grid(id = ids, time = obs_times) |>
    transform(amt = 0, evid = 0L, cmt = "Cc", FORM_CACO3 = form_caco3)
  rbind(dose_rows, obs_rows[order(obs_rows$id, obs_rows$time), ])
}

events_vpc <- bind_rows(
  make_cohort(n_per_arm, form_caco3 = 0, id_offset =        0L) |>
    mutate(treatment = "Thermal spring water"),
  make_cohort(n_per_arm, form_caco3 = 1, id_offset = n_per_arm) |>
    mutate(treatment = "Calcium carbonate")
)
stopifnot(!anyDuplicated(unique(events_vpc[, c("id", "time", "evid")])))

sim_vpc <- rxode2::rxSolve(mod, events = events_vpc,
                           keep = c("treatment", "FORM_CACO3"),
                           returnType = "data.frame")
```

```{r vpc-pth}
vpc_pth <- sim_vpc |>
  filter(time > 0) |>
  group_by(time, treatment) |>
  summarise(
    Q05 = quantile(PTH, 0.05, na.rm = TRUE),
    Q50 = quantile(PTH, 0.50, na.rm = TRUE),
    Q95 = quantile(PTH, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_pth, aes(time, Q50, group = treatment, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.9) +
  labs(x = "Time (hr)", y = "PTH (pg/mL)",
       title = "VPC of PTH by treatment (n=250 per arm)",
       caption = "Replicates the per-treatment lower panel of Ahn 2014 Figure 5. Median +/- 5/95th percentile band.") +
  theme_minimal()
```

```{r vpc-ca}
vpc_ca <- sim_vpc |>
  filter(time > 0) |>
  group_by(time, treatment) |>
  summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(vpc_ca, aes(time, Q50, group = treatment, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.9) +
  labs(x = "Time (hr)", y = "Observed Ca (mmol/L)",
       title = "VPC of observed Ca by treatment (n=250 per arm)",
       caption = "Replicates the per-treatment lower panel of Ahn 2014 Figure 4. Per the paper's K-PD design the observed-Ca distribution is dose-independent.") +
  theme_minimal()
```

The PTH VPC reproduces the published shape: a nadir at ~1-2 hours after the
calcium dose, deeper in the CaCO3 arm, with the median returning toward
baseline by 6-8 hours; the 5/95 percentile band widths match the ~10-40 pg/mL
spread visible in Ahn 2014 Figure 5 lower panel.

The observed-Ca VPC is intentionally featureless across treatments: per the
paper's K-PD design (Methods p. 218 assumption 2), the observed Ca compartment
is decoupled from the calcium dose, so both arms collapse to the same
homeostatic-baseline distribution centered at ~1.19 mmol/L. The residual
spread reflects the additive SD = 0.0258 mmol/L plus the small IIV on
kin_ca (~1.5% CV).

# PTH AUC-equivalent: relative bioavailability validation

The paper's central quantitative claim is that calcium carbonate has ~1.98x
the bioavailability of thermal spring water (Table 2 'Relative F1', 95% CI
1.06-2.90; p < 0.05). Reproduce this claim by computing the area between the
typical PTH trajectory and the drug-free baseline for each treatment arm.

```{r relative-bioavailability}
auc_pth_drop <- function(form_caco3, t_end = 8) {
  ev <- make_typical_events(form_caco3, times = seq(0, t_end, by = 0.01))
  s <- rxode2::rxSolve(mod_typ, events = ev, returnType = "data.frame")
  baseline <- 21.6 / 0.849
  # Area between baseline and the PTH trajectory (positive when PTH is suppressed).
  drop <- pmax(0, baseline - s$PTH)
  sum((drop[-1] + drop[-length(drop)]) / 2 * diff(s$time))
}

auc_thermal <- auc_pth_drop(0)
auc_caco3   <- auc_pth_drop(1)

data.frame(
  treatment           = c("Thermal spring water", "Calcium carbonate"),
  auc_pth_drop_pg_h_per_mL = c(auc_thermal, auc_caco3),
  ratio_vs_thermal    = c(1, auc_caco3 / auc_thermal)
)
```

The ratio of PTH-suppression AUC between the two arms reflects the encoded
1.98x bioavailability multiplier; the AUC ratio is close to 1.98 but not
exactly equal because the Emax inhibition is nonlinear in `ca_unobs`. The
direction (CaCO3 > thermal) and approximate magnitude (~2x) match the paper's
published finding.

# Mass-balance / flux check (algebraic)

At t=0 (no dose, ca_unobs = 0), check that the seeded steady-state baselines
make both indirect-response equations evaluate to zero:

```{r mass-balance}
kin_ca   <- 3.41
kout_ca  <- 2.86
kin_pth  <- 21.6
kout_pth <- 0.849
ec50     <- 0.158
emax     <- 1

ca_ss        <- kin_ca / kout_ca
pth_ss       <- kin_pth / kout_pth
ca_unobs_ss  <- 0
inhib_ss     <- emax * ca_unobs_ss / (ec50 + ca_unobs_ss)

dca       <- kin_ca - kout_ca * ca_ss
dpth      <- kin_pth * (1 - inhib_ss) - kout_pth * pth_ss
dca_unobs <- 0 - kout_ca * ca_unobs_ss

c(ca_ss = ca_ss, pth_ss = pth_ss, dca = dca, dpth = dpth, dca_unobs = dca_unobs)

stopifnot(abs(dca) < 1e-9, abs(dpth) < 1e-9, abs(dca_unobs) < 1e-9)
```

# Assumptions and deviations

* **Compartment naming deviation -- `ca`, `ca_unobs`, `pth`.** The model has
  four states: a depot (canonical), and three paper-specific compartments for
  observed Ca, unobserved Ca, and PTH. The canonical nlmixr2lib compartment
  list (`depot`, `central`, `peripheral1`, ...) does not include endogenous-
  biomarker compartments for ionized calcium or parathyroid hormone, so this
  model uses the paper-faithful names. `checkModelConventions()` emits a WARN
  for each of the three; the warning is expected and intentional. This
  matches the pattern in `Aksenov_2018_uricAcid.R` (`serum`, `urine`) and
  `igg_kim_2006.R` (`igg`).
* **K-PD-style depot.** The depot amount and the unobserved-Ca pool carry no
  physiological mass scale. The user dosess `amt = 1` (Ahn 2014 Methods p. 218
  normalization: '[Ca2+_depot] at t=0 was assumed to be 1'); bioavailability
  scales by the `e_form_caco3_fdepot * FORM_CACO3` exponent of the depot.
  `checkModelConventions()` emits a dimensional-incompatibility WARN between
  the dosing-unit declaration ('unit (normalized)') and the concentration
  numerator (`mmol`); this is expected for K-PD models and is preserved by
  design.
* **No occasion (IOV) construct in the model file.** Ahn 2014 reports an
  inter-occasion variability term on kin_pth of 0.00969 (~9.84% CV between
  Day 1 and Day 7). The packaged model encodes only the IIV component
  (`etalkin_pth ~ 0.0453`, ~21.4% CV); the IOV is documented as a comment in
  `inst/modeldb/endogenous/Ahn_2014_parathyroidHormone.R` `ini()`. For
  occasion-aware simulation (e.g., reproducing the Day 1 vs Day 7 stratified
  panels of Ahn 2014 Figures 4 and 5 upper rows), the user can add an
  occasion-indexed eta on `lkin_pth` following the standard nlmixr2 IOV
  pattern; the variance to use is 0.00969.
* **Published 95% CI on IIV(kin_ca) appears typo'd in Ahn 2014 Table 2.** The
  table reports IIV(kin_ca) = 0.000229 with '95% CI 0.000814, 0.000377', which
  is order-inverted (the lower CI is greater than the upper). The bootstrap
  median (0.00021) and bootstrap 95th percentile range (0.0001-0.0004) confirm
  that the point estimate 0.000229 is correct; the magnitude in the model
  file matches the bootstrap. The CI is reproduced verbatim in source-trace
  comments without any correction (the file reproduces the paper).
* **No drug-coupled observed Ca by design.** Ahn 2014 Methods p. 218
  assumption 2 ('the drop in the PTH response is only due to the increased
  but unobserved Ca2+, as a negative feedback') decouples observed Ca from
  the calcium dose. The observed-Ca compartment has a fixed `kin_ca` (no dose
  term in its ODE) and stays at the homeostatic baseline ~1.19 mmol/L in
  expectation. This is the paper's central K-PD-style assumption, not an
  implementation simplification.
* **Two CaCO3 sub-arms pooled.** The trial randomized 6 subjects to CaCO3 +
  240 mL normal saline and 6 to CaCO3 + 340 mL purified water. The published
  analysis (Methods p. 218) pooled these 12 subjects into a single 'calcium
  carbonate' arm; the packaged model follows this convention with a single
  `FORM_CACO3 = 1` indicator.
* **Cohort generalisability.** The Pottgiesser-equivalent cohort here is
  predominantly male healthy Korean adults (22 M / 2 F) with normal calcium
  and vitamin-D status. Subjects with renal impairment, vitamin-D deficiency
  or supplementation, hyperparathyroidism / hypoparathyroidism, malabsorption
  syndromes, or active calcium-modifying medications are out of scope. The
  Discussion (p. 220-221) explicitly flags that the mild homeostatic
  condition limits extrapolation to extreme physiological perturbations and
  that more sophisticated feedback models (Abraham et al. 2009 reference [6];
  Shrestha et al. 2010 reference [8]) are needed for hypocalcemic-clamp or
  EGTA-chelation regimes.
