---
title: "Evidence for parameter values"
bibliography: articles/references.bib
csl: articles/journal-of-theoretical-biology.csl
link-citations: TRUE
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    code_folding: show
    toc: true
    number_sections: false
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{Evidence for parameter values}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 5,
  fig.asp = NULL,
  out.width = "95%"
)
```

This document collects sources of evidence for baseline parameter values.

```{r setup, class.source = 'fold-hide'}
library(spleenrbc, warn.conflicts = FALSE)
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(here))

output_dir <- here("outputs")

pf <- baseline_parameters("Pf")
pv <- baseline_parameters("Pv")

# Remove non-scalar derived parameters, such as the initial red blood cell
# populations.
remove_nonscalar_derived_parameters <- function(x) {
  x[sapply(x, function(x) length(x) == 1)]
}

# Define how to format parameter values of different magnitudes.
format_value <- function(x, nsmall = 4, digits = 5, min = 1e-4, max = 1e4,
                         delim = TRUE) {
  scientific <- x < min || x > max
  result <- format(x, digits = digits, nsmall = nsmall,
                   scientific = scientific)
  if (delim) {
    start <- "\\("
    end <- "\\)"
  } else {
    start <- ""
    end <- ""
  }
  if (scientific) {
    # Convert from scientific notation to "X x 10^Y".
    parts <- strsplit(result, "e")[[1]]
    digits <- parts[1]
    # Strip the plus sign from exponents.
    power <- sub("^\\+", "", parts[2])
    # Strip leading zeros from positive and negative exponents.
    power <- sub("^(-?)0*", "\\1", power)
    paste0(start, digits, "\\times 10^{", power, "}", end)
  } else {
    # Strip fractional digits for integer values.
    result <- sub("\\.0+$", "", result)
    paste0(start, result, end)
  }
}

# Convenience function for inline R expressions.
pretty <- function(x, digits = 3, nsmall = 0, delim = FALSE) {
  format_value(x, digits = digits, nsmall = nsmall, delim = delim)
}

# Return the baseline parameters for a species as a long data frame.
get_baseline_parameters <- function(species) {
  baseline_parameters(species = species) |>
    remove_nonscalar_derived_parameters() |>
    as.data.frame() |>
    dplyr::mutate(
      species = NULL,
      nospleen = NULL,
      Scenario = NULL
    ) |>
    tidyr::pivot_longer(
      everything(),
      names_to = "Parameter",
    ) |>
    dplyr::mutate(
      value = sapply(value, format_value)
    )
}

# Convert the baseline parameters for both species into data frames.
baseline <- dplyr::inner_join(
  get_baseline_parameters("Pf") |>
    dplyr::rename(Pf = value),
  get_baseline_parameters("Pv") |>
    dplyr::rename(Pv = value),
  by = "Parameter"
) |>
  # Sort by parameter name (case-insensitive).
  dplyr::arrange(tolower(Parameter)) |>
  # Highlight parameters whose values differ between species.
  dplyr::mutate(Parameter = dplyr::case_when(
    .default = Parameter,
    Pf != Pv ~ paste0("**", Parameter, "**")
  ))
```

```{r define-evidence, class.source = 'fold-hide'}
# For each parameter, add a link to the relevant section that describes the
# available evidence (if any).
baseline <- baseline |>
  dplyr::mutate(Evidence = dplyr::case_when(
    Parameter == "gamma" ~ "[Normoblast population in the bone marrow]",
    Parameter == "M_t_1" ~ "[Macrophage population in the spleen]",
    Parameter == "MaxFold" ~ "[Maximum fold increase in RBC production]",
    Parameter == "PMF" ~ "[Parasite multiplication]",
    Parameter == "rho0" ~ "[Reticulocyte release from bone marrow]",
    Parameter == "rho_slope" ~ "[Reticulocyte release from bone marrow]",
    Parameter == "rho_inflection" ~ "[Reticulocyte release from bone marrow]",
    Parameter == "kappa" ~ "[Reticulocyte release from bone marrow]",
    Parameter == "T_R" ~ "[Reticulocyte release from bone marrow]",
    Parameter == "T_R_min" ~ "[Reticulocyte release from bone marrow]",
    Parameter == "U_ss" ~ "[Steady-state RBC population]",
    .default = ""
  ))

knitr::kable(
  baseline,
  caption=paste(
    "Evidence for chosen parameter values.",
    "Parameters that differ between species are shown **in bold**."
  )
)
```

## Macrophage population in the spleen

Steven Kho:

> I cannot find a reference showing this exact figure.
> Our own estimates of total macrophages in the spleen based on counts by histology is \(8.7 \times 10^9\) cells, though our spleens are infected and enlarged, therefore may not represent baseline steady-state.
> I have done similar calculations from published histological counts in control spleens [@Urban05], and the number comes to \(1.2 \times 10^9\) cells.

## Maximum fold increase in RBC production

As per @Watson17:

> In extreme anaemia [RBC production] can be increased fivefold or more, for example.

Figure 3 of Appendix 1 shows the haematocrit response and reticulocyte count response for max-fold values of \(1..10\).

## Normoblast population in the bone marrow

Steven Kho:

> The reference we have previously used for bone marrow cellularity and content is [@Harris62] which states a mean of \(11.1 \times 10^9\) nucleated cells/kg body weight, of which 28.4% are normoblasts.
> Therefore, in an average 60kg Papuan male, this data suggests the bone marrow contains \(6.7 \times 10^{11}\) normoblasts, which is about 2 log-folds higher than the current value in the model.
> If I look at other studies cited in this paper (Table 5), values do fall in this range when you adjust to total bodyweight.

At steady-state we obtain \(\gamma = `r pretty(pf$gamma)`\), which results in a normoblast population of \(`r pretty(pf$gamma * pf$T_n)`\).
This is \(`r pretty(100 * pf$gamma * pf$T_n / 6.7e11 - 100)`\%\) higher than \(6.7 \times 10^{11}\).

### Adjusting the model normoblast population

Presumably we can't assume that the splenectomised Papuans had an average bodyweight of \(\approx 96\) kg.

We can adjust the normoblast population at homeostasis by adjusting model parameters.
However, adjusting the baseline rate parameter \(\rho_0\) and scaling parameter \(\kappa\) for the reticulocyte release rate results in a reduction of \(< 1\%\).
Here are the results obtained by adjusting other parameters:

- \(\rho\): can only increase the normoblast population by decreasing the release rate.
- \(\lambda_u\): minimal change for rate parameter between \(10^{-10}\) and \(10^{10}\).
- \(\delta_u\): minimal change for rate parameter between \(10^{-3}\) and \(1\).
- \(\delta'_u\): minimal change for rate parameter between \(10^5\) and \(10^5\).
- \(\delta_u^{\mathrm{min}}\): a 100-fold reduction decreases the normoblast population to \(\approx 8.5 \times 10^{11}\).

**But** the decrease in \(\delta_u^{\mathrm{min}}\) also causes the circulating RBC population to remain stable during infections.
  With the current baseline parameters, the final circulating RBC count is \(\approx 73.7\%\) of the steady-state count.
  Decreasing \(\delta_u^{\mathrm{min}}\) 100-fold results in a final circulating RBC count that is \(\approx 96.6\%\) of the steady-state count.

Note that the median total RBC count in asymptomatic participants in the HHS differs very little between uninfected (\(1.75 \times 10^{13}\)) and infected (\(1.73 \times 10^{13}\) for Pf, \(1.66 \times 10^{13}\) for Pv).

Steven Kho:

> We reported in the AJH paper [@Kho23] that asymptomatic infection resulted in a 8.1% loss of circulating RBCs, which translates to 91.9% of steady-state.
> Perhaps the rate of removal of RBCs (except young and old) in the spleen could be tweaked to a smaller fold-change to address this?

If we scale \(\delta_u^{\mathrm{min}}\) by \(0.15\), the steady-state normoblast population is \(8.64 \times 10^{11}\), and the final circulating RBC count is \(\approx 8.3\%\) for Pf and \(\approx 9.7\%\) for Pv.

## Parasite multiplication

From @Simpson02:

> The mean population estimate of "PMR" was approximately 8, and was highly dependent on the P. falciparum "strain".
> PMR also varied significantly between patients with a 90% prediction interval varying from 5.5 to 12.3-fold.

## Reticulocyte release from bone marrow

As per the caption of Figure 2 in @Koepke86:

> With increasing anaemia (and erythropoietin production) the maturation time of the erythroid marrow normoblasts and reticulocytes progressively shortens from a normal 3.5 days to 1.5 days or less.
> Conversely the reticulocytes in the peripheral blood persist for a longer time when the patient is anaemic.

The manuscript text also states:

> The maturation time of the reticulocyte in the peripheral blood is taken as 1 day when the PCV is 0.45±0.5, 1.5 days when the PCV is 0.35±0.05, 2 days when the PCV is 0.25±0.05, and 3 days when the PCV is 0.15±0.05 (Hillman & Finch 1967).

Accordingly, the residence time in the bone marrow decreases from 3.5 days when the PCV is 0.45±0.5, down to 1.5 day when the PCV is 0.15±0.05.

The following function provides a reasonable characterisation, while maintaining a release age of 3.5 days in the absence of anaemia:

\begin{align}
  t_r &= \begin{cases}
    T_R & \text{when } \mathbf{U_c} > U_{ss}
    \\
    T_R^{\min} + (T_R - T_R^{\min}) \cdot \left(1 +
      \exp\left[- s \cdot (U_{\mathrm{frac}} - i)\right]\right)^{-1}
    & \text{when } \mathbf{U_c} \le U_{ss}
  \end{cases}
  \\
  U_{\mathrm{frac}} &= \frac{\mathbf{U_c}}{U_{ss}} \\
  s &= 10 \quad \mathrm{(slope)}
  \\
  i &= 0.5 \quad \mathrm{(inflection)}
\end{align}

```{r plot-retic-release, class.source = 'fold-hide'}
orig_retic_release_age <- function(u_frac) {
  T_R <- 3.5 * 24
  T_R_min <- 24

  pkg_data <- new.env()
  utils::data("rbc_steady_state", package = "spleenrbc", envir = pkg_data)
  u_ss <- pkg_data$rbc_steady_state

  u_rbc <- pmin(u_frac * u_ss, u_ss)
  scale <- log(u_ss) - log(u_rbc)
  t_release <- T_R_min + (T_R - T_R_min) * exp(-100 * scale)
  t_release
}

new_retic_release_age <- function(u_frac, inflection, slope) {
  T_R <- 3.5 * 24
  T_R_min <- 24
  u_frac <- pmin(u_frac, 1)

  scale <- T_R - T_R_min
  T_R_min + scale / (1 + exp(- slope * (u_frac - inflection)))
}

plot_retic_release_age <- function() {
  u_frac <- seq(from = 0, to = 1.0, by = 1e-3)
  slope <- 10
  inflection <- 0.5
  orig_age <- orig_retic_release_age(u_frac)
  new_age <- new_retic_release_age(u_frac, inflection, slope)

  df <- data.frame(
    u_pcnt = rep(100 * u_frac, 2),
    age = c(orig_age, new_age),
    equation = rep(c("Original", "New"), each = length(u_frac))
  )

  df_koepke <- data.frame(
    hematocrit = c(45, 35, 25, 15),
    age = 24 * c(3.5, 3.0, 2.5, 1.5)
  ) |> dplyr::mutate(
    u_pcnt = 100 * hematocrit / 45,
    equation = "Koepke & Koepke (1986)"
  )

  ggplot(df, aes(u_pcnt, age / 24)) +
    geom_line(mapping = aes(colour = equation)) +
    geom_point(data = df_koepke) +
    expand_limits(y = 0) +
    scale_x_continuous(
      "uRBC Population (% of steady-state)",
      breaks = c(0, 33, 66.7, 100),
      labels = c("0%", "33%", "67%", "100%")
    ) +
    ylab("Minimum age of release (days)") +
    scale_colour_brewer(
      "Equation: ",
      palette = "Dark2"
    ) +
    theme(legend.position = "top")
}

plot_retic_release_age()
```

## Steady-state RBC population

We used the median RBC count from Papuans with no malaria infection, who had not reported fever in the past 24 hours, as reported by @Pava16.

## References {-}