We defined triangular distributions for 18 models parameters, with bounds of ±20% about the baseline values. We used Latin hypercube sampling (LHS) to draw 1000 parameter combinations from these distributions, and simulated the model for each combination, to assess the model sensitivity.
The uninfected and infected red blood cell (RBC) populations in the circulation and spleen were compared to data collected from asymptomatic Papuan patients who underwent splenectomy. Circulating RBC counts were also compared to data collected from asymptomatic Papuan patients who had not experienced trauma that required splenectomy.
library(spleenrbc, warn.conflicts = FALSE)
library(dplyr, warn.conflicts = FALSE)
library(tidyr)
library(ggplot2)
library(gghalves)
suppressPackageStartupMessages(library(here))
output_dir <- here("outputs")
# Load Pf results.
pf_results <- readRDS(
file.path(output_dir, "sensitivity-analysis-pf-rbc-cints-all.rds")
) |>
mutate(day = time / 24)
# Load Pv results.
pv_results <- readRDS(
file.path(output_dir, "sensitivity-analysis-pv-rbc-cints-all.rds")
) |>
mutate(day = time / 24)
load_clinical_data <- function(species = c("Pf", "Pv")) {
species <- match.arg(species)
# Load HHS total RBC counts for asymptomatic Pf infections.
data(hhs_asymp, envir = environment())
df_asymp <- hhs_asymp |>
filter(species == !!species) |>
mutate(kind = "Uninfected RBCs: Circulation")
# Load uRBC and iRBC counts for splenectomised Pf patients.
data(splenectomised, envir = environment())
df_splen <- splenectomised |>
filter(infection == !!species) |>
pivot_longer(! c("patient", "infection")) |>
filter(! is.na(value)) |>
mutate(
cell_type = case_when(
grepl("^urbc", name) ~ "Uninfected RBCs",
grepl("^irbc", name) ~ "Infected RBCs"
),
location = case_when(
grepl("_spleen", name) ~ "Spleen",
grepl("_peri", name) ~ "Circulation"
),
kind = paste0(cell_type, ": ", location),
cell_type = NULL,
location = NULL
)
# Calculate the ratio of iRBC % in spleen vs circulation.
df_ratio <- df_splen |>
select(patient, infection, name, value) |>
pivot_wider(names_from = "name", values_from = "value") |>
mutate(
ratio = (irbc_spleen / urbc_spleen) / (irbc_peri / urbc_peri),
day = 152
) |>
select(day, ratio) |>
filter(is.finite(ratio))
list(asymp = df_asymp, splen = df_splen, ratio = df_ratio)
}
pf_clinical <- load_clinical_data("Pf")
pv_clinical <- load_clinical_data("Pv")
plot_results_vs_data <- function(df, clinical) {
# Plot uninfected and infected RBC populations against clinical data.
quantiles <- c(0.05, 0.25, 0.5, 0.75, 0.95)
labels <- c(
"U_t" = "Uninfected RBCs: Circulation",
"Ur_t" = "Uninfected RBCs: Spleen",
"I_t" = "Infected RBCs: Circulation",
"Ir_t" = "Infected RBCs: Spleen",
"iRBC_Ratio" = "iRBC: (Spleen % / Circulation %)"
)
df <- df |>
filter(name %in% names(labels)) |>
mutate(kind = labels[name])
df_splen <- bind_rows(
clinical$splen |>
select(value, kind),
clinical$ratio |> mutate(
value = ratio,
kind = "iRBC: (Spleen % / Circulation %)"
) |>
select(value, kind)
)
ggplot(df, aes()) +
geom_half_violin(
data = clinical$asymp,
mapping = aes(2.5, total_rbc),
linetype = "dashed",
draw_quantiles = quantiles,
side = "l",
width = 5
) +
geom_ribbon(
mapping = aes(day, ymin = ci_min, ymax = ci_max,
fill = cred_int, colour = cred_int),
linewidth = 0.5
) +
geom_half_violin(
data = df_splen,
mapping = aes(max(df$day) + 2, value),
position = position_identity(),
draw_quantiles = quantiles,
side = "r",
width = 5
) +
geom_jitter(
data = df_splen,
mapping = aes(max(df$day) + 7, value),
colour = "#af3f3f",
shape = 8,
width = 2
) +
scale_x_continuous(
"Time (days)",
breaks = scales::breaks_width(10)
) +
scale_y_log10("Number of cells") +
scale_colour_brewer(NULL) +
scale_fill_brewer(NULL) +
facet_wrap(~ kind, ncol = 1, scales = "free_y") +
guides(colour = guide_none())
}
plot_results_vs_data(pf_results, pf_clinical)
#> Warning in scale_y_log10("Number of cells"): log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> Warning: Removed 2 rows containing non-finite outside the scale range
#> (`stat_half_ydensity()`).