---
title: "The contribution of host cell-directed vs. parasite-directed immunity to the disease and dynamics of malaria infections" 
subtitle: "Computational data analysis"
author: "Aaron A. King and Nina Wale"
output:
  html_document:
    code_folding: show
    df_print: paged
    highlight: tango
    number_sections: no
    toc: yes
    toc_depth: 4
    toc_float:
      collapsed: yes
      smooth_scroll: yes
  pdf_document:
    toc: yes
    toc_depth: 3
bibliography: ms.bib
csl: pnas.csl
params:
  prefix: decomp
---

```{css echo=FALSE,purl=FALSE}
hr {
	border-width: 3px;
	border-color: #00274c;
}
h2 {
	color: #00274c;
}
h4 {
	font-style: italic;
}
blockquote {
	background-color: #ffeca3;
	border-style: solid;
	border-width: 2px;
	border-color: #00274c;
	padding: 1em;
}
```
-------------------------------

This is the code supplement to the paper ["The contribution of host cell-directed vs. parasite-directed immunity to the disease and dynamics of malaria infections"](https://doi.org/10.1073/pnas.1908147116).


```{r knitr_opts,cache=FALSE,purl=FALSE,include=FALSE}
library(knitr)
opts_chunk$set(
  progress=TRUE,
  prompt=FALSE,tidy=FALSE,highlight=TRUE,
  purl=FALSE,
  strip.white=TRUE,
  warning=FALSE,
  message=FALSE,
  error=FALSE,
  echo=TRUE,
  cache=TRUE,
  cache.extra=rand_seed,
  dev="png",
  dev.args=c(bg="white"),
  dpi=150,
  results='markup',
  fig.show='asis',
  size='small',
  fig.align='left',
  fig.path=paste0("figure/",params$prefix,"-"),
  cache.path=paste0("cache/",params$prefix,"-"))
```
```{r setup,message=FALSE,cache=FALSE,purl=TRUE}
library(plyr)
library(tidyverse)
library(magrittr)
library(tibble)
library(stringi)
library(pomp)
library(panelPomp)
library(foreach)
library(iterators)
library(doRNG)
library(aakmisc) ## available at https://kingaa.github.io/
stopifnot(packageVersion("pomp")>="2.2")
stopifnot(packageVersion("panelPomp")>="0.9.1")
stopifnot(packageVersion("aakmisc")>="0.26.2")
options(
  stringsAsFactors=FALSE,
  keep.source=TRUE,
  encoding="UTF-8"
)
set.seed(407958184)
```
```{r plot_pkgs,message=FALSE,cache=FALSE,purl=FALSE}
library(ggplot2)
library(ggpmisc)
library(RColorBrewer)
library(plotly)
library(cowplot)
library(scales)
theme_set(theme_bw())
```
-------------------------------

## Experimental data

The following codes import and pre-process the data.
Experimental design and protocols are described in the [paper](https://doi.org/10.1073/pnas.1908147116).
The data themselves are available in this archive, in the file named `data.csv`.
The data consist of the following fields, measured daily on each of 15 mice.

- `day`: day post-infection
- `box`: mice are kept in a box, according to their treatment. 
  This value is the unique identifier of each box.
- `mouse`: unique identifier for each individual mouse in a box
- `ama_density`: Density of parasites per microliter of blood; 
  calculated from CT.
- `rbc_density`: total density of red blood cells, of all stages of maturity, as derived from a coulter counter.
- `ter119_density`: density of cells in the erythroid family, as derived by flow cytometry
- `cd71_density`: density of reticulocytes, as derived by flow cytometry

```{r import_data,purl=TRUE}
read_csv("data.csv",
         col_types="iiinnnn"
) %>%
  mutate(
    mouseid=sprintf("%02d-%02d",box,mouse),
    box=sprintf("%02d",box),
    mouse=sprintf("%02d",mouse),
    paba=as.factor(box),
    rbc_density=rbc_density/1000
  ) %>%
  mutate(
    paba=recode(
      paba,
      "01"="0.05","02"="0.005","03"="0.0005","04"="0","05"="control"
    )
  ) %>%
  select(
    day,
    Pd=ama_density,
    RBC=rbc_density,
    Ter119=ter119_density,
    CD71=cd71_density,
    mouseid,
    paba,box,mouse
  ) %>%
  arrange(mouseid,day) %>%
  mutate(
    paba=as.character(paba),
    Ter119=ifelse(Ter119==0,NA,Ter119),
    CD71=ifelse(CD71==0,NA,CD71)
  ) %>%
  mutate(
    Eryth=(1-CD71/Ter119)*RBC,
    Retic=CD71/Ter119*RBC
  ) -> flow
```

```{r plot_data,fig.height=10,purl=FALSE}
flow %>%
  gather(variable,value,-day,-mouseid,-paba,-mouse,-box) %>% 
  ggplot(aes(x=day,y=value,group=mouseid,color=paba))+
  geom_line()+
  facet_wrap(~variable,scales="free_y",ncol=1)+
  scale_y_log10()+
  scale_colour_manual(
    name="pABA concentration (%)", 
    breaks=c("0.05", "0.005", "0.0005", "0", "control"), 
    values=c("turquoise3","purple", "darkgreen",  "darkorange2", "darkgrey"))

flow %>%
  gather(variable,value,-day,-mouseid,-paba,-mouse,-box) %>%
  ggplot(aes(x=day,y=value,group=mouseid,color=mouseid))+
  geom_line()+
  facet_grid(variable~paba,scales="free_y")+
  scale_y_log10()
```

-------------------------------

## Estimation of inoculum and burst sizes

We estimate the dose (inoculum size) and burst-size parameter, $\beta$, by regressing the log of parasite growth on day of infection, restricted to the first 4 days.
The slope is our estimate of $\beta$; 
the intercept is the logarithm of the inoculum size.

We first fit a model with a common slope and then a model with a slope for each pABA treatment.
Model comparison indicates that the latter is well supported by the data.
Inspection of the coefficients for the various treatments reveals that they increase with pABA concentration.

```{r beta-dose-plot,purl=FALSE}
flow %>%
  select(day,Pd,mouseid,paba) %>%
  filter(!paba=="control")%>%
  ggplot(aes(x=day,y=Pd,group=mouseid,color=paba))+
  geom_line()+
  scale_y_log10()+
  scale_colour_manual(
    name="pABA concentration (%)", 
    breaks=c("0.05", "0.005", "0.0005", "0"),
    values=c("turquoise3","purple", "darkgreen",  "darkorange2"))

flow %>% 
  filter(day<=4) %>%
  filter(!paba=="control")%>%
  select(day,Pd,mouseid,paba) %>%
  ggplot(aes(x=day,y=Pd,group=mouseid,color=paba))+
  geom_line()+
  scale_y_log10()+
  scale_colour_manual(
    name="pABA concentration (%)",
    breaks=c("0.05", "0.005", "0.0005", "0"),
    values=c("turquoise3","purple", "darkgreen",  "darkorange2"))
```

```{r beta-dose,purl=TRUE}
flow %>%
  filter(day<=4) %>%
  lm(log(Pd)~day+mouseid-1,data=.) -> fit1
summary(fit1)

flow %>%
  filter(day<=4) %>%
  lm(log(Pd)~box:day+mouseid-1,data=.) -> fit2
summary(fit2)

anova(fit1,fit2)
```

> **From the text:**
As expected, the initial growth rate of infections increased with the concentration of pABA that they were administered.


```{r beta-dose2}
coef(fit2) %>% 
  bind_rows() %>%
  gather(var,val) %>%
  mutate(
    var=stri_replace_all_regex(var,"mouseid(\\d{2})-(\\d{2})","dose[$1-$2]"),
    var=stri_replace_all_regex(var,"box(\\d{2}):day","Beta[$1]"),
    val=exp(val)
  ) -> theta1

expand.grid(
  box=sprintf("%02d",1:5),
  mouse=sprintf("%02d",1:3)
) %>%
  mutate(
    mouseid=paste0(box,"-",mouse),
    betavar=sprintf("Beta[%s]",box),
    dosevar=sprintf("dose[%s]",mouseid)
  ) %>%
  left_join(theta1,by=c("betavar"="var")) %>%
  rename(Beta=val) %>%
  left_join(theta1,by=c("dosevar"="var")) %>%
  rename(dose=val) %>%
  select(-betavar,-dosevar) %>%
  arrange(box,mouse) %>%
  mutate(
    Beta=coalesce(Beta,0),
    dose=coalesce(dose,0)
  ) -> theta

theta
```

-------------------------------

## Model implementation

The model equations are:
$$
\begin{aligned}
K_t &= \left(R_t+E_t\right)\,\left(1-\exp{\left(-\frac{M_t}{R_t+E_t}\right)}\right)\\
M_{t+1} &= \beta\,K_t\,\exp{\left(-\frac{W_t+N_t}{R_t+E_t}\right)}\\
E_{t+1} &= \left(R_t+E_t\right)\,\exp{\left(-\frac{M_t+N_t}{R_t+E_t}\right)}
\end{aligned}
$$

The state variables $K$, $M$, and $E$ are linked to the observed variables <tt>Pd</tt>, <tt>RBC</tt>, and <tt>Retic</tt> via a lognormal measurement model.

The following codes implement the model using the **pomp** package [@pomp;@King2016].
Documentation on this package is available at [the **pomp** website](https://kingaa.github.io/pomp/).

```{r build_pomp,cache=FALSE,purl=TRUE}
library(doParallel)
registerDoParallel()

theta %>%
  mutate(
    sigmaPd = 2, 
    sigmaRBC = 0.1, 
    sigmaRetic = 0.3, 
    sigmaW = 1, 
    sigmaN = 0.5, 
    sigmaR = 0.5,
    E_0 = 8.0e6, 
    R_0 = 3e5, 
    W_0 = 8.8e4, 
    N_0 = 8e5
  ) -> theta

foreach (m = iter(theta,"row"),.inorder=TRUE,.combine=c) %dopar% {
  
  flow %>%
    filter(mouseid==m$mouseid) %>%
    select(day,Pd,RBC,Retic) %>%
    mutate(Retic=if_else(day %in% c(0,14),NA_real_,Retic)) %>%
    pomp(
      params=select(m,-mouseid,-box,-mouse) %>% unlist(),
      times="day",
      t0=0,
      rmeasure=Csnippet("
  Retic = rlnorm(log(1+R),sigmaRetic)-1;
  RBC = rlnorm(log(1+E+R),sigmaRBC)-1;
  Pd = rlnorm(log(1+K),sigmaPd)-1;"),
  dmeasure=Csnippet("
  double l1, l2, l3;
  l1 = (R_FINITE(Retic)) ? dlnorm(1+Retic,log(1+R),sigmaRetic,1) : 0;
  l2 = (R_FINITE(RBC)) ? dlnorm(1+RBC,log(1+E+R),sigmaRBC,1) : 0;
  l3 = (R_FINITE(Pd) && Pd>0) ? dlnorm(1+Pd,log(1+K),sigmaPd,1) : 0;
  lik = (give_log) ? l1+l2+l3 : exp(l1+l2+l3);"),
  rprocess=discrete_time(
    step.fun=Csnippet("
  double Mold = M;
  M = Beta*K*exp(-(W+N)/(R+E));
  E = (R+E)*exp(-(Mold+N)/(R+E));
  N = rlnorm(log(N),sigmaN);
  W = rlnorm(log(W),sigmaW);
  R = rlnorm(log(R),sigmaR);
  K = (R+E>0) ? (R+E)*(1-exp(-M/(R+E))): 0;
  "),
  delta.t=1
  ),
  partrans=parameter_trans(
    log=c("sigmaW","sigmaR","sigmaN",
          "sigmaPd","sigmaRBC","sigmaRetic",
          "N_0","W_0","E_0","R_0")
  ),
  rinit=Csnippet("
  E = E_0;
  R = R_0;
  N = N_0;
  W = W_0;
  M = 0;
  K = dose;"),
  statenames=c("E","R","W","N","M","K"),
  paramnames=c(
    "Beta","dose",
    "sigmaPd","sigmaRBC","sigmaRetic",
    "sigmaW","sigmaN","sigmaR",
    "E_0","R_0","W_0","N_0"
  )
    )
} %>% 
  set_names(theta$mouseid) -> pos
```

-------------------------------

## Parameter estimation


```{r mpi_setup,purl=TRUE,echo=FALSE}
if (file.exists("CLUSTER")) {
  scan("CLUSTER",what=integer(0)) -> ncpu
  library(doMPI)
  cl <- startMPIcluster(ncpu,verbose=TRUE,logdir="/tmp")
  registerDoMPI(cl)
}
```

-------------------------------

### Initial iterated filtering

To narrow the focus of our search for the maximum likelihood parameter estimate (MLE), we begin by running the iterated filtering algorithm [@Ionides:2015dw] on the data from each individual mouse.
For each mouse, we run 50 independent `mif2` computations, initiating each at a distinct point.
These start points are distributed within a large hyperrectangle in the 10-dimensional parameter space.

```{r shared_params,purl=TRUE}
shared_params <- c(
  "E_0","N_0","R_0",
  "sigmaN","sigmaPd","sigmaR",
  "sigmaRBC","sigmaRetic","sigmaW",
  "W_0"
)
```

```{r mf1,cache=FALSE,purl=TRUE}
bake(file="m5mf1.rds",{
  
  theta %>%
    select(-box,-mouse,-mouseid,-Beta,-dose) %>%
    gather(var,val) %>%
    group_by(var) %>%
    summarize(min=min(val),max=max(val)) %>% 
    ungroup() %>%
    gather(bound,val,-var) %>%
    spread(var,val) %>%
    column_to_rownames("bound") %>%
    as.matrix() -> box
  
  profileDesign(
    mouseid=theta$mouseid,
    lower=0.5*box["min",shared_params],
    upper=1.5*box["max",shared_params],
    nprof=50
  ) -> starts
  
  registerDoRNG(2022051968)
  
  foreach (start=iter(starts,"row"),
           .inorder=FALSE,.combine=bind_rows,
           .packages=c("pomp","magrittr","plyr","dplyr")
  ) %dopar%
    {
      tryCatch({
        pos[[start$mouseid]] -> m
        start %>% select(-mouseid) %>% unlist() -> theta
        coef(m,names(theta)) <- unname(theta)
        start <- cbind(start["mouseid"],as.data.frame(as.list(coef(m))))
        m %>%
          mif2(Nmif=50,Np=10000,tol=exp(-60),
               cooling.fraction.50=0.1,cooling.type="geometric",
               rw.sd=rw.sd(
                 sigmaPd=0.02,sigmaRBC=0.02,sigmaRetic=0.02,
                 sigmaW=0.02,sigmaN=0.02,sigmaR=0.02,
                 E_0=ivp(0.05),R_0=ivp(0.05),W_0=ivp(0.05),N_0=ivp(0.05)
               )
          ) %>%
          mif2(Nmif=100) -> m
        raply(5,m %>% pfilter(Np=50000,tol=exp(-60)) %>% logLik()) %>% 
          logmeanexp(se=TRUE) -> pfs
        m %>% coef() %>% rbind() %>% as.data.frame() -> fin
        fin$loglik <- pfs[[1]]
        fin$loglik.se <- pfs[[2]]
        fin$mouseid <- start$mouseid
        fin$msg <- "mf1"
        start$msg <- "start"
        bind_rows(start,fin)
      },
      error=function (e) {
        fin <- start
        fin$msg <- conditionMessage(e)
        full_join(start,fin,by=names(start))
      })
    } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  
  res
}) -> mf1
```

```{r plotfuns,purl=FALSE}
library(GGally)

sigmaplot1 <- function (data, ..., ll = TRUE) {
  
  if (missing(..1)) {
    ae <- aes(color=msg)
  } else {
    ae <- aes(...)
  }
  
  cols <- c("sigmaW","sigmaN","sigmaR","sigmaRBC","sigmaRetic","sigmaPd")
  collabs <- c("sigma[W]","sigma[N]","sigma[R]","sigma[RBC]","sigma[Re]","sigma[Pd]")
  if (ll) {
    cols <- c(cols,"loglik")
    collabs <- c(collabs,"log(L)")
  }
  
  data %>%
    ggpairs(mapping=ae,upper=NULL,legend=1,
            lower=list(continuous=wrap("points",alpha=0.5,size=0.2)),
            diag=list(continuous=wrap("densityDiag",alpha=0.5)),
            columns=cols,columnLabels=collabs,labeller="label_parsed")+
    theme(legend.position="bottom",axis.text.x=element_text(angle=-90))
}

icplot1 <- function (data, ..., ll = TRUE) {
  
  if (missing(..1)) {
    ae <- aes(color=msg)
  } else {
    ae <- aes(...)
  }
  
  cols <- c("E_0","R_0","W_0","N_0")
  collabs <- c("E[0]","R[0]","W[0]","N[0]")
  if (ll) {
    cols <- c(cols,"loglik")
    collabs <- c(collabs,"log(L)")
  }
  
  data %>%
    ggpairs(mapping=ae,upper=NULL,legend=1,
            lower=list(continuous=wrap("points",alpha=0.5,size=0.2)),
            diag=list(continuous=wrap("densityDiag",alpha=0.5)),
            columns=cols,columnLabels=collabs,labeller="label_parsed")+
    theme(legend.position="bottom",axis.text.x=element_text(angle=-90))
}

```

We display the results of this initial search.

```{r mf1-plot1,fig.height=6,fig.width=6,purl=FALSE}
mf1 %>%
  group_by(mouseid) %>%
  mutate(loglik=loglik-max(loglik,na.rm=TRUE)) %>%
  ungroup() %>%
  filter(is.na(loglik) | loglik>-10) %>%
  sigmaplot1()+
  labs(title="individual mice")

mf1 %>%
  group_by(mouseid) %>%
  mutate(loglik=loglik-max(loglik,na.rm=TRUE)) %>%
  ungroup() %>%
  filter(is.na(loglik) | loglik>-10) %>%
  icplot1()+
  labs(title="individual mice")

mf1 %>%
  group_by(mouseid) %>%
  mutate(loglik=loglik-max(loglik,na.rm=TRUE)) %>%
  ungroup() %>%
  filter(!is.na(loglik) | loglik>-10) %>%
  sigmaplot1(color=mouseid)+
  labs(title="best likelihoods")

mf1 %>%
  group_by(mouseid) %>%
  mutate(loglik=loglik-max(loglik,na.rm=TRUE)) %>%
  ungroup() %>%
  filter(!is.na(loglik) | loglik>-10) %>%
  separate(mouseid,into=c("box","mouse"),remove=FALSE) %>%
  filter(box != "05") %>%
  sigmaplot1(color=mouseid)+
  labs(title="best likelihoods, excluding controls")

```

-------------------------------

### Panel iterated filtering

At the next stage of parameter estimation, we pool the data across mice.
In particular, we treat the full data set as a panel of data and use panel iterated filtering [@Breto2019], as implemented in the [**panelPomp** package](https://cbreto.github.io/panelPomp/), to estimate parameters.
We run 250 independent panel iterated filtering computations, initiating each at a different point within the hyper-rectangle spanned by the estimates from the previous stage.

Naturally, we exclude the controls from this exercise.

```{r exclude}
control_mice <- grep("^05-",names(pos),perl=TRUE)
```

```{r mf2,cache=FALSE,purl=TRUE}
bake(file="m5mf2.rds",{
  
  mf1 %>%
    separate(mouseid,into=c("box","mouse"),remove=FALSE) %>%
    ## remove controls and mouse that died early
    filter(box != "05" & mouseid != "01-01") %>%
    group_by(mouseid) %>%
    filter(!is.na(loglik) | loglik>max(loglik,na.rm=TRUE)-10) %>%
    ungroup() %>%
    select(starts_with("sigma"),ends_with("_0")) %>%
    gather(var,val) %>%
    group_by(var) %>%
    summarize(min=min(val),max=max(val)) %>% 
    ungroup() %>%
    gather(bound,val,-var) %>%
    spread(var,val) %>%
    column_to_rownames("bound") %>%
    as.matrix() -> box
  
  sobolDesign(
    lower=box["min",shared_params],
    upper=box["max",shared_params],
    nseq=250
  ) -> starts

  registerDoRNG(1510289681)
  
  foreach (start=iter(starts,"row"),
           .inorder=FALSE,.combine=bind_rows,
           .packages=c("panelPomp","magrittr","plyr","dplyr"),
           .options.mpi=list(chunkSize=1)) %dopar%
    {
      tryCatch({
        panelPomp(pos[-control_mice],shared=unlist(start)) -> m
        m %>% coef() %>% rbind() %>% as.data.frame() -> start
        m %>%
          mif2(Nmif=300,Np=20000,tol=exp(-180),
               cooling.fraction.50=0.75,cooling.type="geometric",
               rw.sd=rw.sd(
                 sigmaPd=0.01,sigmaRBC=0.01,sigmaRetic=0.01,
                 sigmaW=0.01,sigmaN=0.01,sigmaR=0.01,
                 E_0=ivp(0.02),R_0=ivp(0.02),W_0=ivp(0.02),N_0=ivp(0.02)
               )
          ) -> m
        m %>% coef() %>% rbind() %>% as.data.frame() -> fin
        fin$msg <- "mf2"
        start$msg <- "mf1"
        bind_rows(start,fin)
      },
      error=function (e) {
        fin <- start
        fin$msg <- conditionMessage(e)
        full_join(start,fin,by=names(start))
      })
    } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  
  res
}) -> mf2
```

The above required `r signif(attr(mf2,"system.time")[3]/3600,2)` hr on `r attr(mf2,"ncpu")` cpus.

```{r mf2-plot1,purl=FALSE}
mf2 %>% 
  sigmaplot1(ll=FALSE)+
  labs(title="panel IF2, round 1")

mf2 %>% 
  icplot1(ll=FALSE)+
  labs(title="panel IF2, round 1")
```

Having examined the results from the first round of panel iterated filtering, we perform a second round.
This is initiated at the endpoint of the previous calculation.

```{r mf3,cache=FALSE,purl=TRUE}
bake(file="m5mf3.rds",{
  
  mf2 %>%
    filter(msg=="mf2") %>%
    arrange(sigmaW,sigmaN,sigmaR) %>%
    select(-msg) -> starts
  
  registerDoRNG(599108918)
  
  foreach (start=iter(starts,"row"),
           .inorder=FALSE,.combine=bind_rows,
           .packages=c("panelPomp","magrittr","plyr","dplyr"),
           .options.mpi=list(chunkSize=1)) %dopar%
    {
      tryCatch({
        panelPomp(pos[-control_mice],shared=unlist(start)) -> m
        m %>% coef() %>% rbind() %>% as.data.frame() -> start
        m %>%
          mif2(Nmif=100,Np=20000,tol=exp(-180),
               cooling.fraction.50=0.1,cooling.type="geometric",
               rw.sd=rw.sd(
                 sigmaPd=0.01,sigmaRBC=0.01,sigmaRetic=0.01,
                 sigmaW=0.01,sigmaN=0.01,sigmaR=0.01,
                 E_0=ivp(0.02),R_0=ivp(0.02),W_0=ivp(0.02),N_0=ivp(0.02)
               )
          ) -> m
        m %>% coef() %>% rbind() %>% as.data.frame() -> fin
        fin$msg <- "mf3"
        start$msg <- "mf2"
        bind_rows(start,fin)
      },
      error=function (e) {
        fin <- start
        fin$msg <- conditionMessage(e)
        full_join(start,fin,by=names(start))
      })
    } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  
  res
}) -> mf3
```

The above required `r signif(attr(mf3,"system.time")[3]/3600,2)` hr on `r attr(mf3,"ncpu")` cpus.


```{r mf3-plot1,purl=FALSE}
mf3 %>% sigmaplot1(ll=FALSE)+
  labs(title="panel IF2, round 2")

mf3 %>% icplot1(ll=FALSE)+
  labs(title="panel IF2, round 2")

```

At the stopping points, we compute the likelihood via the particle filter.

```{r pf3,cache=FALSE,purl=TRUE}
bake(file="m5pf3.rds",{
  
  bind_rows(mf2,mf3) %>%
    distinct() %>%
    arrange(sigmaW,sigmaN,sigmaR) -> starts
  
  registerDoRNG(865465456)
  
  foreach (start=iter(starts,"row"),
           .inorder=FALSE,.combine=bind_rows,
           .packages=c("panelPomp","magrittr","plyr","dplyr"),
           .options.mpi=list(chunkSize=1)) %dopar%
    {
      tryCatch({
        start %>% select(-msg) %>% unlist() -> p
        panelPomp(pos[-control_mice],shared=p) -> m
        replicate(
          n=10,
          m %>% pfilter(Np=150000,tol=exp(-180))
        ) %>% 
          sapply(logLik) %>%
          logmeanexp(se=TRUE) -> lls
        start$loglik <- lls[1]
        start$loglik.se <- lls[2]
        start
      },
      error=function (e) {
        start$msg <- conditionMessage(e)
        start
      })
    } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  
  res
}) -> pf3
```

The above required `r signif(attr(pf3,"system.time")[3]/3600,2)` hr on `r attr(pf3,"ncpu")` cpus.

```{r pf3-plot1,purl=FALSE}
pf3 %>% 
  sigmaplot1()+
  labs(title="panel iterated filtering results")

pf3 %>% 
  icplot1()+
  labs(title="panel iterated filtering results")

pf3 %>%
  filter(loglik+2*loglik.se>max(loglik-2*loglik.se)-100) %>%
  sigmaplot1()+
  labs(title="best likelihoods from panel iterated filtering")

pf3 %>% 
  filter(loglik+2*loglik.se>max(loglik-2*loglik.se)-100) %>%
  icplot1()+
  labs(title="best likelihoods from panel iterated filtering")
```

```{r pf3-2}
pf3 %>%
  group_by(msg) %>%
  filter(loglik+2*loglik.se>max(loglik-2*loglik.se)-100) %>%
  summarize(n=length(loglik),max=max(loglik),
            mean=mean(loglik),sd=sd(loglik),meanSE=mean(loglik.se))
```

-------------------------------

### Profile likelihood

To further refine the estimates, we profile over over $\sigma_W$.
This eliminates one direction of search, but requires considerably more effort.

```{r mf4,purl=TRUE,cache=FALSE}
bake(file="m5mf4.rds",{
  
  shared_params <- c(
    "E_0","N_0","R_0",
    "sigmaN","sigmaPd","sigmaR",
    "sigmaRBC","sigmaRetic",
    "W_0"
  )
  
  pf3 %>%
    as_tibble() %>%
    filter(loglik>max(loglik)-100) %>%
    select(starts_with("sigma"),ends_with("_0")) %>%
    gather(var,val) %>%
    group_by(var) %>%
    summarize(min=min(val),max=max(val)) %>%
    ungroup() %>%
    gather(bound,val,-var) %>%
    spread(var,val) %>%
    column_to_rownames("bound") %>%
    as.matrix() -> box 
  
  profileDesign(
    sigmaW=seq(from=0.2,to=box["max","sigmaW"],length=100),
    lower=box["min",shared_params],
    upper=box["max",shared_params],
    nprof=10
  ) -> starts
  
  registerDoRNG(933950321)
  
  foreach (
    start=iter(starts,"row"),
    .combine=bind_rows,
    .packages=c("plyr","dplyr","magrittr","panelPomp")
  ) %dopar% {
    
    start %>% unlist() -> p
    panelPomp(pos[-control_mice],shared=p) -> m
    m %>%
      mif2(Nmif=100,Np=10000,tol=exp(-180),
           rw.sd=rw.sd(
             sigmaPd=0.02,sigmaRBC=0.02,sigmaRetic=0.02,
             sigmaW=0,sigmaN=0.02,sigmaR=0.02,
             E_0=ivp(0.02),R_0=ivp(0.02),W_0=ivp(0.02),N_0=ivp(0.02)
           ),
           cooling.type="geometric",cooling.fraction.50=0.1) %>%
      mif2(Np=20000) %>%
      mif2(Np=40000) -> m
    m %>% coef() %>% rbind() %>% as.data.frame() %>%
      cbind(data.frame(msg="mf4"))
  } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  res
  
}) -> mf4
```

The above required `r signif(attr(mf4,"system.time")[3]/3600,2)` hr on `r attr(mf4,"ncpu")` cpus.

As before, we compute the likelihood at each stopping point.

```{r pf4,purl=TRUE,cache=FALSE}
bake(file="m5pf4.rds",{
  
  mf4 %>%
    arrange(sigmaW,sigmaN,sigmaR) %>%
    mutate(id=row_number(msg)) -> starts
  
  registerDoRNG(537031349)
  
  foreach (
    start=iter(starts,"row"),
    .combine=bind_rows,
    .packages=c("plyr","dplyr","tibble","magrittr","panelPomp")
  ) %dopar% {
    
    start %>% select(-msg) %>% unlist() -> p
    panelPomp(pos[-control_mice],shared=p) -> m
    replicate(n=10,
              m %>% pfilter(Np=150000,tol=exp(-180)) %>% logLik()
    ) -> lls
    start %>% 
      remove_rownames() %>%
      cbind(data.frame(loglik=lls))
    
  } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  res
  
}) -> pf4
```

The above required `r signif(attr(pf4,"system.time")[3]/3600,2)` hr on `r attr(pf4,"ncpu")` cpus.

```{r pf4-plot1,purl=FALSE}
pf4 %>%
  bind_rows(pf3) %>%
  ggplot(aes(x=sigmaW,y=loglik,color=msg,group=1))+
  geom_point(alpha=0.5)+
  lims(y=c(-11500,NA))

pf4 %>%
  group_by(id) %>%
  mutate(mean=mean(loglik)) %>%
  ungroup() %>%
  mutate(cw=cut(sigmaW,breaks=seq(0,max(sigmaW)+0.02,by=0.1))) %>%
  group_by(cw) %>%
  filter(mean>max(mean)-10) %>%
  ungroup() %>%
  select(-cw,-mean) %>%
  ggplot(aes(x=sigmaW,y=loglik,color=msg,group=1))+
  geom_point(alpha=0.5)+
  geom_smooth(method="loess")

pf4 %>%
  group_by(id) %>%
  mutate(mean=mean(loglik)) %>%
  ungroup() %>%
  mutate(cw=cut(sigmaW,breaks=seq(0,max(sigmaW)+0.02,by=0.1))) %>%
  group_by(cw) %>%
  filter(mean>max(mean)-10) %>%
  ungroup() %>%
  select(-cw,-mean) %>%
  ggplot(aes(x=sigmaW,y=loglik,color=msg,group=1))+
  geom_point(alpha=0.5)+
  geom_smooth(method="loess")+
  lims(y=c(-100,NA)+max(pf4$loglik))

pf4 %>%
  group_by(id) %>%
  mutate(loglik=mean(loglik)) %>%
  ungroup() %>%
  mutate(cw=cut(sigmaW,breaks=seq(0,max(sigmaW)+0.02,by=0.1))) %>%
  group_by(cw) %>%
  filter(loglik==max(loglik)) %>%
  ungroup() %>%
  select(-cw) %>%
  sigmaplot1()+
  labs(title="profile over sigmaW")

pf4 %>%
  group_by(id) %>%
  mutate(mean=mean(loglik)) %>%
  ungroup() %>%
  mutate(cw=cut(sigmaW,breaks=seq(0,max(sigmaW)+0.02,by=0.1))) %>%
  group_by(cw) %>%
  filter(mean>max(mean)-10) %>%
  ungroup() %>%
  select(-cw,-mean) %>%
  icplot1()+
  labs(title="profile over sigmaW")
```

-------------------------------

### Maximum likelihood estimate

As a point estimate, we simply take the highest likelihood found.

```{r mle}
pf4 %>%
  select(loglik,starts_with("sigma"),ends_with("_0")) %>%
  filter(loglik==max(loglik)) -> mle

mle %>%
  gather(var,val,-loglik) %>%
  mutate(val=signif(val,3)) %>%
  spread(var,val) %>%
  select(loglik,starts_with("sigma"),ends_with("_0"))
```

-------------------------------

## Smoothing

The following codes draw samples from the smoothed distibution of the state variables.

```{r sm1,cache=FALSE,purl=TRUE}
bake(file="m5sm1.rds",{
  
  registerDoRNG(1510516029)
  
  foreach (i=1:2000, 
           .inorder=FALSE,.combine=bind_rows,
           .packages=c("panelPomp","magrittr","reshape2","plyr","tibble","tidyr","dplyr")
  ) %dopar% 
    {
      
      mle %>% select(-loglik) %>% unlist() -> p
      panelPomp(pos,shared=p) -> m
      
      m %>%
        pfilter(Np=1e5,filter.traj=TRUE,tol=exp(-180)) %>%
        as("list") %>%
        llply(filter.traj) %>%
        melt() %>%
        rename(mouseid=L1) %>%
        mutate(rep=i)
      
    } -> res
  attr(res,"ncpu") <- getDoParWorkers()
  res
  
}) -> sm1
```

The above required `r signif(attr(sm1,"system.time")[3]/3600,2)` hr on `r attr(sm1,"ncpu")` cpus.

In the [Supporting Information](https://doi.org/10.1073/pnas.1908147116), the following quantities are defined.
$$
\begin{aligned}
S^{M}_t &= \exp{\left(-\frac{M_{t}}{R_{t}+E_{t}}\right)}\\
S^{N}_t &= \exp{\left(-\frac{N_{t}}{R_{t}+E_{t}}\right)}\\
S^{W}_t &= \exp{\left(-\frac{W_{t}}{R_{t}+E_{t}}\right)}\\
Q^{ps}_{t} &= \left(1-S^{M}_{t}\right)\,S^{W}_{t}\,\,S^{N}_{t}\\
Q^{pn}_{t} &= \frac{N_{t}}{N_{t}+W_{t}}\,\left(1-S^{M}_{t}\right)\,\left(1-S^{W}_{t}\,S^{N}_{t}\right)\\
Q^{pw}_{t} &= \frac{W_{t}}{N_{t}+W_{t}}\,\left(1-S^{M}_{t}\right)\,\left(1-S^{W}_{t}\,S^{N}_{t}\right)\\
Q^{un}_{t} &= S^{M}_{t}\,\left(1-S^{N}_{t}\right)\\ 
Q^{us}_{t} &= S^{M}_{t}\,S^{N}_{t}\\
\lambda^r_{t} &= \beta\,\frac{R_{t}}{M_{t}}\,Q^{ps}_{t}\\
\lambda^e_{t} &= \beta\,\frac{E_{t}}{M_{t}}\,Q^{ps}_{t}\\
\lambda^n_{t} &= \beta\,\frac{R_{t}+E_{t}}{M_{t}}\,Q^{pn}_{t}\\
\lambda^w_{t} &= \beta\,\frac{R_{t}+E_{t}}{M_{t}}\,Q^{pw}_{t}\\
\lambda^u_{t} &= \beta-\lambda^r_{t}-\lambda^e_{t}-\lambda^n_{t}-\lambda^w_{t}\\
\end{aligned}
$$

We compute median and 95% confidence intervals of these from the smoothing-distribution samples as follows.

```{r jfc,purl=FALSE}
mle %>% unlist() -> p
pos %>% panelPomp(shared=p) %>% coef() %>% 
  rbind() %>% as_tibble() -> theta
theta %>%
  select(starts_with("Beta"),starts_with("dose")) %>%
  gather(variable, value) %>%
  tidyr::extract(variable, into=c("variable","mouseid"), 
                 regex="([[:alnum:]]+)\\[(.+?)\\]") %>%
  spread(variable,value) -> bdf

sm1 %>%
  as_tibble() %>%
  spread(variable,value) %>%
  left_join(bdf,by=c("mouseid")) %>% 
  group_by(mouseid) %>%
  mutate(
    SM=exp(-M/(R+E)),
    SN=exp(-N/(R+E)),
    SW=exp(-W/(R+E)),
    Qps=(1-SM)*SW*SN,
    Qpn=N/(N+W)*(1-SM)*(1-SW*SN),
    Qpw=W/(N+W)*(1-SM)*(1-SW*SN),
    Qun=SM*(1-SN),
    Qus=SM*SN,
    lambda_r=Beta*R/M*Qps,
    lambda_e=Beta*E/M*Qps,
    lambda_n=Beta*(R+E)/M*Qpn,
    lambda_w=Beta*(R+E)/M*Qpw,
    lambda_u=Beta-lambda_r-lambda_e-lambda_n-lambda_w,
    loss=(R+E)*(1-Qus),
    rbc=E+R
  ) %>%
  gather(variable,value,-rep,-time,-mouseid) %>%
  group_by(mouseid,time,variable) %>%
  summarize(
    lo=quantile(value,prob=0.05,na.rm=TRUE),
    hi=quantile(value,prob=0.95,na.rm=TRUE),
    med=mean(value,na.rm=TRUE)
  ) %>%
  ungroup() %>%
  right_join(
    flow %>%
      select(mouseid,box,mouse) %>%
      distinct(),
    by="mouseid"
  ) -> jfc
```


-------------------------------

## Figures

Here, we provide the codes that generate the figures of the paper.

-------------------------------

### Figure 2

```{r fig2,purl=FALSE,fig.height = 3.2,fig.width = 3.2}
flow %>%
  select(time=day, P=Pd, RBC, R=Retic, mouseid) %>%
  mutate(rbc=RBC) %>%
  select(-RBC) %>%
  gather(variable, med, -mouseid, -time) %>%
  mutate(type="data") %>%
  filter(!variable=="R" | !time%in%c(0,14)) %>%
  filter(time<22) -> datplot

jfc %>%
  filter(variable %in% c("rbc","K","R")) %>%
  mutate(
    variable=recode(variable,K="P",rbc="rbc"),
    type="model",
    lo=trnc(lo,c(1,1e8)),
    med=trnc(med,c(1,1e8)),
    hi=trnc(hi,c(1,1e8))) %>%
  select(mouseid, time, variable, med, hi, lo, type) %>%
  full_join(datplot,by=c("mouseid","time","variable","med","type")) %>%
  filter(mouseid%in%c("02-01","03-01","04-03") & time<22) -> x

data_labels <- c(
  "P" =  "Parasitized RBCs",
  "rbc"  = "Total RBCs",
  "R" = "Reticulocytes")

mouse_labels <- c(
  "02-01" = "mouse 1",
  "03-01" = "mouse 2",
  "04-03" = "mouse 3")

x %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+   
  scale_x_continuous(limits=c(0,30), breaks=seq(0,30,5))+
  scale_color_manual(name="",labels=c("data","model"),values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","orange"), 
                    guide=FALSE)+
  facet_grid(
    variable~mouseid, 
    scales="free_y", 
    labeller=labeller(variable=data_labels))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title= element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size=6, family="Helvetica", angle=0),
    strip.text.x = element_blank(),
    strip.text.y = element_text(size=6, family="Helvetica", angle=0, hjust = 0),
    strip.background = element_blank(),
    legend.text=element_text(size=6, family="Helvetica"),
    legend.background = element_rect(fill=alpha('white', 0)),
    legend.key = element_rect(fill = NA),
    legend.position = c(1, 0),
    legend.justification = c(1, 0),
    legend.title=element_blank()
  )+
  guides(colour = guide_legend(override.aes = list(shape = NA))) -> data_l


x %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+  
  scale_x_continuous(limits=c(0,30), breaks=seq(0,30,5))+
  scale_color_manual(
    name="",labels=c("data","model"),
    values=c("black","orange"))+
  scale_fill_manual(
    name="",values=c("black","orange"), guide=FALSE)+
  facet_grid(
    variable~mouseid,
    scales="free_y",
    labeller=labeller(variable=data_labels, mouseid=mouse_labels))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title= element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size=6, family="Helvetica", angle=0),
    strip.text.y = element_text(size=9, family="Helvetica",angle=0),
    strip.text.x = element_blank(),
    strip.background = element_blank(),
    legend.text=element_text(size=6, family="Helvetica"),
    legend.background = element_rect(fill=alpha('white', 0)),
    legend.key = element_rect(fill = NA),
    legend.position = c(1, 0),
    legend.justification = c(1, 0),
    legend.title=element_blank()
  )+
  guides(colour = guide_legend(override.aes = list(shape = NA))) -> data_mousel

jfc %>%
  filter(
    variable %in% c("W","R","N","K", "rbc") & 
      mouseid%in%c("02-01","03-01","04-03")) %>%
  mutate(variable=factor(variable, levels=c("K", "rbc", "R","N", "W"))) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=variable,color=variable,fill=variable))+
  geom_line(aes(linetype=variable))+
  geom_ribbon(alpha=0.2,color=NA)+
  labs(x="day",y="")+
  scale_linetype_manual(
    values=c("dashed", rep("solid",4)),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing", "Targeted killing"))+
  scale_color_manual(
    values=c("grey","grey","pink","royal blue","darkorchid4"),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing", "Targeted killing"))+
  scale_fill_manual(
    values=c("grey","grey","pink","royal blue","darkorchid4"),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing", "Targeted killing"), 
    guide="none")+ 
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)),
                limits=c(1e4,1e8),oob=trnc)+
  scale_x_continuous(limits=c(0,30), breaks=seq(0,30,5))+
  facet_grid(~mouseid)+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    panel.border = element_rect(size=0.25),
    strip.text.x = element_blank(),
    strip.text.y = element_text(size=6, family="Helvetica", angle=0, hjust = 0),
    strip.background = element_blank(),
    legend.text=element_text(size=6, family="Helvetica"),
    legend.background = element_rect(fill=alpha('white', 0)),
    axis.text.y = element_text(size=7, family="Helvetica", angle=0),
    legend.justification = c(1, 0),
    legend.title=element_blank()
  ) -> traj_l

data.frame(
  variable = c("P", "P", "R", "R", "rbc", "rbc"), 
  time = 0, 
  med = c(1,1e7,1e5,1e7,1e6,1e7),
  lo=NA,
  hi=NA,
  color=NA, 
  type="data"
) -> blank_data
## To define limits of free scales, use these data to make a geom_blank 
## following http://chrischizinski.github.io/rstats/2014/08/24/using_geom_blank/

x %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line(size=0.25)+
  geom_point(data=filter(x, type=="data"), size=0.5)+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+  
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  scale_color_manual(name="",labels=c("data","model"),values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","orange"), 
                    guide=FALSE)+
  facet_grid(variable~mouseid, scales="free_y",
             labeller=labeller(variable=data_labels))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title= element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size=6, family="Helvetica", angle=0),
    strip.text = element_blank(),
    strip.background = element_blank(),
    legend.text=element_text(size=6, family="Helvetica"),
    legend.background = element_rect(fill=alpha('white', 0)),
    legend.key = element_rect(fill = NA),
    legend.key.size = unit(0.5, "lines"),
    legend.position = c(1, 0),
    legend.justification = c(1,0),
    legend.title=element_blank(),
    panel.border = element_rect(size=0.25)
  )+
  geom_blank(data = blank_data, aes(x = time, y = med))+
  guides(colour = guide_legend(override.aes = list(shape = NA))) -> data

jfc %>%
  filter(variable %in% c("W","R","N","K", "rbc") & 
           mouseid%in%c("02-01","03-01","04-03")) %>%
  mutate(variable=factor(variable, levels=c("K", "rbc", "R","N", "W"))) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=variable,color=variable,fill=variable))+
  geom_line(aes(linetype=variable))+
  geom_ribbon(alpha=0.2,color=NA)+
  scale_linetype_manual(
    values=c("dashed", rep("solid",4)),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing","Parasite-directed killing"))+
  scale_color_manual(
    values=c("grey","grey","pink","royal blue","darkorchid4"),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing","Parasite-directed killing"))+
  scale_fill_manual(
    values=c("grey","grey","pink","royal blue","darkorchid4"),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing","Parasite-directed killing"), 
    guide="none")+ 
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x), 
    labels=trans_format('log10',math_format(10^.x)),
    limits=c(1e4,1e8),oob=trnc)+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  facet_grid(~mouseid)+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.text.x = element_blank(),
    strip.text.y = element_text(size=6, family="Helvetica", angle=0),
    strip.background = element_blank(),
    legend.text=element_text(size=6, family="Helvetica"),
    legend.background = element_rect(fill=alpha('white', 0)),
    axis.text = element_text(size=6, family="Helvetica", angle=0),
    legend.justification = c(1, 1),
    legend.title=element_blank(),
    legend.position= "none",
    axis.title=element_blank(),
    panel.border = element_rect(size=0.25)
  ) -> traj

library(grid)
print(data,vp=viewport(y=1/3+1/3,height=2/3))
print(traj,vp=viewport(y=1/6,height=1/3))
```

The figure in the main text has slight differences:
some strictly stylistic elements were edited.

-------------------------------

### Figure 3

```{r fig3,fig.height = 7, fig.width = 9.1, purl=FALSE}
mouse_labels <- c(
  "02-01" = "mouse 1",
  "03-01" = "mouse 2",
  "04-03" = "mouse 3"
)

jfc %>%
  filter(
    variable == "Beta",
    mouseid%in%c("02-01","03-01","04-03"), 
    time==0
  ) -> maxrepro

jfc %>%
  filter(
    variable %in% c("lambda_w", "lambda_n","lambda_r","lambda_e", "lambda_u"),
    mouseid%in%c("02-01","03-01","04-03")
  ) %>%
  mutate(
    variable=factor(
      variable, 
      levels=c("lambda_u","lambda_w","lambda_n","lambda_r","lambda_e"))
  ) %>% 
  ggplot(aes(x=time,y=med,group=variable,fill=variable))+
  geom_area(color=NA)+
  geom_hline(data=maxrepro, aes(yintercept=med), linetype="dashed")+
  geom_hline(aes(yintercept=1))+
  geom_vline(aes(xintercept=6), color="white")+
  labs(x="day",y= "offspring/ parasite")+
  scale_fill_manual(
    name="",
    values=c("chocolate4", "midnight blue","steelblue1","mistyrose", "light pink"), 
    labels= c("unrealized reproduction", "destroyed by targeted killing", 
              "destroyed by indiscriminate killing",
              "emerged from reticulocytes", "emerged from erythrocytes"))+
  facet_grid(~mouseid, labeller=labeller(mouseid=mouse_labels))+
  theme(
    strip.text.x = element_text(size=14, family="Helvetica"),
    strip.text.y= element_blank(),
    strip.background = element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.text.x=element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=14),
    legend.justification=c(0,1),
    legend.text = element_text(size=10),
    legend.title = element_blank()
  )+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,by=5)) -> fig3p

jfc %>%
  filter(
    variable %in% c("Qpn", "Qpw", "Qun", "Qps", "Qus"),
    mouseid%in%c("02-01","03-01","04-03")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  gather(variable,value,-c(mouseid, time, box, mouse)) %>%
  mutate(
    variable=factor(variable,levels=c("Qun","Qps","Qpn", "Qpw","Qus"))
  ) -> fig3x

fig3x %>%
  ggplot()+
  geom_area(aes(x=time, y=value, fill=variable, group=variable))+
  geom_vline(aes(xintercept=6), color="white")+
  labs(x="day",y="fraction of RBCs")+ 
  scale_fill_manual(
    name="",
    values=c("seagreen","royal blue", "steelblue1",
             "midnight blue", "blanchedalmond"),
    labels=c("unparasitized & destroyed by indiscriminate killing", 
             "parasitized & destroyed by emerging parasites",
             "parasitized & destroyed by indiscriminate killing", 
             "parasitized & destroyed by targeted killing",
             "survived"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  facet_grid(~mouseid)+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.title = element_text(size=14),
    legend.text = element_text(size=10),
    legend.justification=c(0,1),
    legend.title = element_blank()
  )+
  guides(fill = guide_legend(reverse=F, title=NULL)) -> fig3h

jfc %>%
  filter( mouseid%in%c("02-01","03-01","04-03")) %>%
  select(-c(lo,hi)) %>%
  spread(variable,med) %>%
  mutate(
    leadR=lead(R),
    calc=lead(R-loss),
    turnover=pmin(leadR,loss),
    surplus=pmax(leadR-turnover,0),
    deficit=pmax(loss-turnover,0)
  ) -> a

a %>%
  select(time,mouseid,turnover,surplus,deficit) %>%
  gather(variable,value,-time,-mouseid) %>%
  mutate(variable=factor(variable,levels=c("surplus","deficit","turnover"))) %>%
  ggplot(aes(x=time,y=value,fill=variable,color=variable))+
  geom_bar(stat="identity")+
  scale_fill_manual(values=c(surplus="black",deficit="red",turnover="white"))+
  scale_color_manual(values=c(surplus=NA,deficit=NA,turnover="beige"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(labels=aakmisc::scinot)+
  facet_wrap(~mouseid)+
  labs(fill="",color="",x="",y="RBC/ \u03BCl blood")+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    axis.text.x=element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=14),
    legend.text=element_text(size=10),
    legend.justification=c(0,1),
    legend.title = element_blank()
  ) -> fig3rep

plot_grid(
  fig3p,fig3rep,fig3h, 
  ncol = 1, 
  align="v", 
  axis="lr", 
  rel_heights=c(1.1,0.7,1),
  labels=c("A","B","C")
) -> fig3

fig3
```

-------------------------------

### Figure 4

```{r fig4,purl=FALSE, fig.width = 7}
mouse_labels <- c(
  "02-01" = "mouse 1",
  "03-01" = "mouse 2",
  "04-03" = "mouse 3")

jfc %>%
  filter(mouseid%in%c("02-01","03-01","04-03") & time<22) %>%
  select(-c(lo,hi)) %>%
  spread(variable,med) %>%
  ddply(.(mouseid), mutate, 
        landmark_a=time[which.max(K)],
        landmark_b=time[which.max(N)]) %>%
  select("mouseid","time","K","R","E","N","W",starts_with("landmark_")) %>%
  mutate(landmark_a=case_when(
    mouseid=="04-03" ~ 7L,
    TRUE ~ landmark_a)) -> ia

ggplot(data=ia, aes(x=K,y=W,group=mouseid, color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed",
                        length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(ia, time==0), aes(fill=mouseid), shape=21, color="black")+
  geom_point(data=filter(ia, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(ia, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(ia,time==landmark_a), label="a", color="black",
            nudge_x=0.2, nudge_y=0.1)+
  geom_text(data=filter(ia,time==landmark_b), label="b", color="black",
            nudge_x=-0.2)+
  labs(x="Parasite density (K)",y="Targeted killing (W)")+
  theme(
    axis.title=element_text(size=10, family="Helvetica"), 
    axis.text= element_text(size=8, family="Helvetica"),
    legend.position="none",
    strip.text=element_blank(),
    plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm"))+
  scale_color_manual(values=rep("goldenrod4",3))+
  scale_fill_manual(values=rep("goldenrod4",3))+
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)), limits=c(10^5.3,10^7.5))+
  scale_x_log10(
    limits=c(1e2,1e7),breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)), oob=trnc)+
  facet_wrap(~mouseid,nrow=1) -> kw

ggplot(data=ia, aes(x=R,y=W,group=mouseid, color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed",
                        length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(ia, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(ia, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(ia,time==landmark_a), label="a",
            color="black", nudge_y=0.1)+
  geom_text(data=filter(ia,time==landmark_b), label="b", 
            color="black", nudge_x=0.09)+
  geom_point(data=filter(ia, time==0), aes(fill=mouseid), shape=21, color="black")+
  labs(y="Targeted killing (W)",x="Reticulocyte supply (R)")+
  theme(axis.title=element_text(size=10, family="Helvetica"), 
        axis.text= element_text(size=8, family="Helvetica"),
        legend.position="none",
        strip.text=element_blank(),
        plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm"))+
  scale_color_manual(values=rep("goldenrod4",3))+
  scale_fill_manual(values=rep("goldenrod4",3))+
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)), 
    limits=c(10^5.3,10^7.5))+
  scale_x_log10(
    limits=c(10^5,10^6.75),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)))+
  facet_wrap(~mouseid,nrow=1) -> rw

ggplot(data=ia, aes(x=R,y=N,group=mouseid, color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(ia, time==0), aes(fill=mouseid), shape=21, color="black")+
  geom_point(data=filter(ia, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(ia, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(ia,time==landmark_a), label="a", color="black", 
            nudge_x=-0.075)+
  geom_text(data=filter(ia,time==landmark_b), label="b", color="black",
            nudge_x=0.075)+
  labs(x="Reticulocyte supply (R)",y="Indiscriminate killing (N)")+
  theme(
    axis.title=element_text(size=10, family="Helvetica"), 
    axis.text= element_text(size=8, family="Helvetica"),
    legend.position="none",
    strip.text=element_blank(),
    plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm"))+
  scale_color_manual(values=rep("goldenrod4",3))+
  scale_fill_manual(values=rep("goldenrod4",3))+
  scale_y_log10(
    limits=c(10^5,10^6.9),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)),
    oob=trnc)+
  scale_x_log10(
    limits=c(1e5,10^6.75),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)))+
  facet_wrap(~mouseid,nrow=1) -> rn

ggplot()+
  geom_segment(aes(x = 0, y = 4, xend = 5.7, yend = 4, colour = "up"), size=2)+
  geom_segment(aes(x = 7.3, y = 4, xend = 9.7, yend = 4, colour = "down",
                   linetype="all"), size=2)+
  geom_segment(aes(x = 12.3, y = 4, xend = 15, yend = 4, colour = "up",
                   linetype="some"),size=2)+
  geom_segment(aes(x = 15, y = 4, xend = 21, yend = 4, colour = "down"), size=2)+
  
  geom_segment(aes(x = 0, y = 3, xend = 5.7, yend = 3, colour = "up"), size=2)+
  geom_segment(aes(x = 7.3, y = 3, xend = 9.7, yend = 3, colour = "down"), size=2)+
  geom_segment(aes(x = 12.3, y = 3, xend = 21, yend = 3, colour = "up"),  size=2)+
  
  geom_segment(aes(x = 0, y = 2, xend = 5.7, yend = 2, colour = "up"), size=2)+
  geom_segment(aes(x = 7.3, y = 2, xend = 9.7, yend = 2, colour = "up"), size=2)+
  geom_segment(aes(x = 12.3, y = 2, xend = 21, yend = 2, colour = "down"), size=2)+
  
  geom_segment(aes(x = 0, y=1, xend = 5.7, yend = 1, colour = "down"), size=2)+
  geom_segment(aes(x = 7.3, y=1, xend = 9.7, yend = 1, colour = "up"), size=2)+
  geom_segment(aes(x = 12.3, y=1, xend = 21, yend = 1, colour = "down"), size=2)+
  geom_text(aes(x=6.5,y=c(1,2,3,4),label="a"),size=5)+
  geom_text(aes(x=11,y=c(1,2,3,4),label="b"),size=5)+
  scale_color_manual("",values=c("blue","red"), breaks=c("up","down"), 
                     labels = c("increasing", "decreasing"))+
  scale_linetype_manual("",values=c("solid","dotted"), breaks=c("all","some"),
                        guide= "none")+
  scale_y_continuous(
    name="",
    labels=c("Reticulocyte\nsupply (R)","Indiscriminate\nkilling (N)",
             "Targeted\nkilling (W)","Parasite\ndensity (K)"), 
    limits=c(0.5,4.25), 
    breaks=c(1,2,3,4))+
  scale_x_continuous(
    name="Time (days)", 
    breaks=c(0,6.5,11,21), 
    labels=c("0","6 - 7","10 - 12","21"))+
  theme(
    axis.text.y = element_text(hjust=0,size=10, color="black"),
    axis.text.x = element_text(size=10, color="black"),
    axis.ticks.x=element_blank()
  ) -> schematic

ggplot(data=ia, aes(x=K,y=W,group=mouseid, color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(ia, time==0), aes(fill=mouseid), shape=21, color="black")+
  geom_point(data=filter(ia, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(ia, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(ia,time==landmark_a), label="a", color="black",
            nudge_x=0.2, nudge_y=0.1)+
  geom_text(data=filter(ia,time==landmark_b), label="b", color="black", 
            nudge_x=-0.1)+
  labs(x="Parasite density (K)",y="Targeted killing (W)")+
  theme(axis.title=element_text(size=10, family="Helvetica"), 
        axis.text= element_text(size=8, family="Helvetica"),
        legend.position="none",
        strip.text.y = element_blank(),
        strip.text.x = element_text(size=10, family="Helvetica"),
        strip.background = element_blank(),
        plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm"))+
  scale_color_manual(values=rep("goldenrod4",3))+
  scale_fill_manual(values=rep("goldenrod4",3))+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels=trans_format('log10',math_format(10^.x)), 
                limits=c(10^5.3,10^7.5))+
  scale_x_log10(limits=c(1e2,1e7),
                breaks = trans_breaks("log10", function(x) 10^x),
                labels=trans_format('log10',math_format(10^.x)), oob=trnc)+
  facet_grid(~mouseid,labeller=labeller(mouseid=mouse_labels)) -> swirl_facet

plot_grid(
  kw,rw, rn, schematic, 
  ncol = 1, 
  align="v", 
  axis="lr", 
  rel_heights=c(1,1,1,1),
  labels=c("A","B", "C", "D")) -> fig4_combo

plot_grid(
  swirl_facet,rw, rn, schematic, 
  ncol = 1, 
  align="v", 
  axis="lr", 
  rel_heights=c(1,1,1,1)) -> fig4_facetlabels

fig4_combo
```

The figure in the main text has some slight differences:
some strictly stylistic elements were edited.

-------------------------------

### Figure S1

```{r figS1,purl=FALSE, fig.width=6,fig.height=9}
flow %>%
  select(time=day, P=Pd, RBC, R=Retic, mouseid, box) %>%
  mutate(rbc=RBC) %>%
  select(-RBC) %>%
  gather(variable, med, -mouseid, -time, -box) %>%
  mutate(type="data")  %>%
  filter(!variable=="R" | !time%in%c(0,14)) %>%
  filter(time<22) -> datplot

jfc %>%
  filter(variable %in% c("rbc","K","R")) %>%
  mutate(
    variable=recode(variable,K="P",rbc="rbc"),
    type="model",
    lo=trnc(lo,c(1,1e8)),
    med=trnc(med,c(1,1e8)),
    hi=trnc(hi,c(1,1e8))) %>%
  select(mouseid, time, variable, med, hi, lo, type) %>%
  full_join(datplot,by=c("mouseid","time","variable","med","type")) -> x

data.frame(
  variable = c("P", "P", "R", "R", "rbc", "rbc"), time = 0, 
  med = c(1,1e7,1e5,1e7,1e6,1e7),
  lo=NA, hi=NA, color =NA, type="data"
) -> blank_data

data_labels <- c(
  "P" =  "P",
  "rbc"  = "RBCs",
  "R" = "R")

x %>%
  filter(mouseid%in%c("01-01", "01-02", "01-03")) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  geom_point(
    data=filter(
      x, 
      type=="data" & mouseid%in%c("01-01", "01-02", "01-03")))+
  labs(x="day",y="")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  scale_color_manual(name="",labels=c("data","model"),
                     values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","dark orange"), guide=FALSE)+
  facet_grid(variable~mouseid, scales="free_y")+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    plot.background=element_rect(color="#993404"),
    strip.text = element_blank(),
    axis.title= element_blank(),
    axis.text.x= element_blank(),
    legend.title=element_blank(),
    legend.position="none"
  )+
  geom_text(data = filter(x, mouseid == "01-01" & variable=="P"),
            aes(x= Inf, y=Inf, label="†"), hjust = 1.5, vjust=1.5)+
  geom_text(data = filter(x, mouseid == "01-02" & variable=="P"),
            aes(x= Inf, y=Inf, label="V"), hjust = 1.5, vjust=1.5)+
  geom_blank(data = blank_data, aes(x = time, y = med)) -> highdm

#0.005
x %>%
  filter(mouseid%in%c("02-01", "02-02", "02-03")) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,  group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  geom_point(data=filter(x, type=="data" & mouseid%in%c("02-01", "02-02", "02-03")))+
  labs(x="day",y="")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  scale_color_manual(name="",labels=c("data","model"),
                     values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","dark orange"), guide=FALSE)+
  facet_grid(
    variable~mouseid, scales="free_y", 
    labeller=labeller(variable=data_labels))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    plot.background=element_rect(color="#d95f0e"),
    strip.text.x = element_blank(),
    strip.text.y= element_text(family="Helvetica", angle=0),
    strip.background = element_blank(),
    axis.title= element_blank(),
    axis.text = element_blank(),
    legend.title=element_blank(),
    legend.position="none"
  )+
  geom_text(data = filter(x, mouseid == "02-03" & variable=="P"),
            aes(x= Inf, y=Inf, label="V"), hjust = 1.5, vjust=1.5)+
  geom_blank(data = blank_data, aes(x = time, y = med)) -> lowdm

#0.005
x %>%
  filter(mouseid%in%c("03-01", "03-02", "03-03")) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,  group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  geom_point(data=filter(x, type=="data" & mouseid%in%c("03-01", "03-02", "03-03")))+
  labs(x="day",y="")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  scale_color_manual(name="",labels=c("data","model"),
                     values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","dark orange"), guide=FALSE)+
  facet_grid(variable~mouseid, scales="free_y")+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    plot.background=element_rect(color="#fe9929"),
    strip.text = element_blank(),
    axis.title= element_blank(),
    legend.title=element_blank(),
    legend.position="none"
  )+
  geom_blank(data = blank_data, aes(x = time, y = med)) -> ultralowdm

x %>%
  filter(mouseid%in%c("04-01", "04-02", "04-03")) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,  group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  geom_point(data=filter(x, type=="data" & mouseid%in%c("04-01", "04-02", "04-03")))+
  labs(x="day",y="")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  scale_color_manual(name="",labels=c("data","model"),
                     values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","dark orange"), guide=FALSE)+
  facet_grid(variable~mouseid, scales="free_y",labeller=labeller(variable=data_labels))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    plot.background=element_rect(color="#fed98e"),
    strip.text.x = element_blank(),
    strip.text.y= element_text(family="Helvetica", angle=0),
    strip.background = element_blank(),    axis.title= element_blank(),
    axis.text.y= element_blank(),
    legend.title=element_blank(),
    legend.position="none"
  )+
  geom_blank(data = blank_data, aes(x = time, y = med)) -> nonedm

data.frame(
  variable = c("R", "R", "rbc", "rbc"), time = 0,
  med = c(1e5,1e7,1e6,1e7),
  lo=NA, hi=NA, color =NA, type="data"
) -> blank_data_uninf

data_labels_uninf <- c(
  "rbc"  = "RBCs",
  "R" = "R"
)

x %>%
  filter(mouseid%in%c("05-01", "05-02", "05-03") & !variable=="P") %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,  group=type,color=type,fill=type))+
  geom_ribbon(alpha=0.3, color="NA")+
  geom_line()+
  geom_point(
    data=filter(
      x,
      type=="data" & mouseid%in%c("05-01", "05-02", "05-03") & !variable=="P"))+
  labs(x="day",y="")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), 
                labels=trans_format('log10',math_format(10^.x)))+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  scale_color_manual(name="",labels=c("data","model"),
                     values=c("black","orange"))+
  scale_fill_manual(name="",values=c("black","dark orange"), guide=FALSE)+
  facet_grid(
    variable~mouseid, scales="free_y",
    labeller=labeller(variable=data_labels_uninf))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    plot.background=element_rect(linetype="dashed", color="black"),
    strip.text.x = element_blank(),
    strip.text.y= element_text(family="Helvetica", angle=0),
    strip.background = element_blank(),     axis.title.y= element_blank(),
    legend.title=element_blank(),
    legend.position="none"
  )+
  geom_blank(data = blank_data_uninf, aes(x = time, y = med)) -> uninfdm

ggplot()+
  geom_blank()+
  theme(panel.border = element_blank()) -> g

plot_grid(
  highdm, ultralowdm, 
  ncol = 1, 
  nrow = 2,
  align="v") -> parasitized_left

plot_grid(
  lowdm, nonedm,
  ncol = 1, 
  nrow = 2,
  align="v") -> parasitized_right

plot_grid(parasitized_left, parasitized_right) -> test

plot_grid(
  uninfdm,g,
  nrow=1,
  rel_widths =c(1,0.73)) -> bottom_row


plot_grid(
  test,bottom_row,
  nrow=2,
  rel_heights=c(1,0.4)) -> all_dm

all_dm
```

-------------------------------

### Figure S2

```{r figS2,purl=FALSE, fig.width=6, fig.height=5}
data_labels <- c(
  "01" =  "0.05%",
  "02"  = "0.005%",
  "03" = "0.0005%",
  "04" = "0%")

jfc %>%
  filter(variable %in% c("W","R","N","K", "rbc") & box!="05" & time<22) %>%
  mutate(variable=factor(variable, levels=c("K", "rbc", "R","N", "W")),
         lo=trnc(lo,c(1,1e8)),
         med=trnc(med,c(1,1e8)),
         hi=trnc(hi,c(1,1e8))) %>%
  ggplot(aes(x=time,y=med,ymin=lo,ymax=hi,group=variable,color=variable,fill=variable))+
  geom_line(aes(linetype=variable))+
  geom_ribbon(alpha=0.2,color=NA)+
  labs(x="day",y="density/ \u03BCl blood")+
  scale_linetype_manual(
    values=c("dashed", rep("solid",4)),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes",
             "Indiscriminate killing", "Targeted killing"))+
  scale_color_manual(
    values=c("grey","grey","pink","royal blue","darkorchid4"),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing", "Targeted killing"))+
  scale_fill_manual(
    values=c("grey","grey","pink","royal blue","darkorchid4"),
    labels=c("Parasitized RBCs","Total RBCs","Reticulocytes", 
             "Indiscriminate killing", "Targeted killing"), 
    guide="none")+  
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x), 
    labels=trans_format('log10',math_format(10^.x)),
    limits=c(1e4,1e8),oob=trnc)+
  scale_x_continuous(limits=c(0,21), breaks=seq(0,20,5))+
  theme(
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.text.x = element_blank(),
    strip.text.y = element_text(size=8, family="Helvetica", angle=0),
    strip.background = element_blank(),
    legend.text=element_text(size=8, family="Helvetica"),
    legend.background = element_rect(fill=alpha('white', 0)),
    axis.text = element_text(size=8, family="Helvetica", angle=0),
    legend.justification = c(1, 0),
    legend.title=element_blank(),
    axis.title=element_text(size=10, family="Helvetica", angle=0),
    panel.border = element_rect(size=0.25)
  )+
  facet_grid(box~mouse, labeller=labeller(box=data_labels)) -> figS2

figS2
```

-------------------------------

### Figure S3 

```{r figS3,purl=FALSE, fig.width=9, fig.height=9.1}
jfc %>%
  filter(time<22) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  group_by(mouseid,box) %>%
  summarize(
    time_peakK=time[which.max(K)]
  ) -> peak

jfc %>%
  filter(
    variable %in% c("lambda_w", "lambda_n","lambda_r","lambda_e", "lambda_u"),
    box%in%c("01","02")
  ) %>%
  left_join(peak, by=c("box","mouseid")) %>%
  mutate(
    variable=factor(
      variable, 
      levels=c("lambda_u","lambda_w","lambda_n","lambda_r","lambda_e"))
  ) %>%
  ggplot(aes(x=time,y=med,group=variable,fill=variable))+
  geom_area(color=NA)+
  geom_hline(aes(yintercept=1))+
  geom_vline(aes(xintercept=time_peakK), color="white")+
  labs(x="day",y= "offspring/ parasite")+
  scale_fill_manual(
    name="",
    values=c("chocolate4", "midnight blue","steelblue1","mistyrose", "light pink"), 
    labels= c("unrealized reproduction", "destroyed by targeted killing", 
              "destroyed by indiscriminate killing","emerged from reticulocytes",
              "emerged from erythrocytes"))+
  facet_grid(~mouseid)+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.text.y=element_text(size=6),
    axis.text.x=element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=8),
    legend.justification=c(0,1),
    legend.text = element_text(size=),
    legend.title = element_blank(),
    legend.position = "none"
  )+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(limits=c(0,6.6),breaks=seq(0,6,by=2)) -> p_12

jfc %>%
  filter(box%in%c("01","02")) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  mutate(leadR=lead(R),
         calc=lead(R-loss),
         turnover=pmin(leadR,loss),
         surplus=pmax(leadR-turnover,0),
         deficit=pmax(loss-turnover,0)) %>%
  select(time,mouseid,turnover,surplus,deficit) %>%
  gather(variable,value,-time,-mouseid) %>%
  mutate(variable=factor(variable,levels=c("surplus","deficit","turnover"))) %>%
  ggplot(aes(x=time,y=value,fill=variable,color=variable))+
  geom_bar(stat="identity")+
  scale_fill_manual(values=c(surplus="black",deficit="red",turnover="white"))+
  scale_color_manual(values=c(surplus=NA,deficit=NA,turnover="beige"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(limits=c(0,6e6),labels=aakmisc::scinot)+
  facet_grid(~mouseid)+
  labs(fill="",color="",x="",y="RBC/ \u03BCl blood")+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    axis.text.x=element_blank(),
    axis.text.y=element_text(size=6),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=8),
    legend.text=element_text(size=8),
    legend.justification=c(0,1),
    legend.title = element_blank(),
    legend.position = "none"
  ) -> rep_12

jfc %>%
  filter(
    variable %in% c("Qpn", "Qpw", "Qun", "Qps", "Qus"),
    box%in%c("01","02")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>% 
  gather(variable,value,-c(mouseid, time, box, mouse)) %>%
  mutate(variable=factor(variable,levels=c("Qun","Qps","Qpn", "Qpw","Qus"))) %>%
  left_join(peak, by=c("box","mouseid")) %>%
  ggplot()+
  geom_area(aes(x=time, y=value, fill=variable, group=variable))+
  geom_vline(aes(xintercept=time_peakK), color="white")+
  labs(x="day",y="fraction of RBCs")+ 
  scale_fill_manual(
    name="",
    values=c("seagreen","royal blue", "steelblue1",  "midnight blue", "blanchedalmond"), 
    labels=c(
      "unparasitized & destroyed by indiscriminate killing", 
      "parasitized & destroyed by emerging parasites", 
      "parasitized & destroyed by indiscriminate killing",
      "parasitized & destroyed by targeted killing",
      "survived"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  facet_grid(~mouseid)+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size=6),
    axis.title.x = element_blank(),
    axis.title = element_text(size=8),
    legend.text = element_text(size=6),
    legend.justification=c(0,1),
    legend.title = element_blank(),
    legend.position = "none"
  )+
  guides(fill = guide_legend(reverse=F, title=NULL)) -> d_12

plot_grid(
  p_12,rep_12,d_12,
  ncol = 1, 
  align="v", 
  axis="lr", 
  rel_heights=c(1,0.7,1)
) -> top_row

jfc %>%
  filter(
    variable %in% c("lambda_w", "lambda_n","lambda_r","lambda_e", "lambda_u"),
    box%in%c("03","04")
  ) %>%
  left_join(peak, by=c("box","mouseid")) %>%
  mutate(
    variable=factor(
      variable, 
      levels=c("lambda_u","lambda_w","lambda_n","lambda_r","lambda_e"))
  ) %>%
  ggplot(aes(x=time,y=med,group=variable,fill=variable))+
  geom_area(color=NA)+
  geom_hline(aes(yintercept=1))+
  geom_vline(aes(xintercept=time_peakK), color="white")+
  labs(x="day",y= "offspring/ parasite")+
  scale_fill_manual(
    name="",
    values=c("chocolate4", "midnight blue","steelblue1","mistyrose", "light pink"), 
    labels= c("unrealized reproduction", "destroyed by targeted killing", 
              "destroyed by indiscriminate killing",
              "emerged from reticulocytes", "emerged from erythrocytes"))+
  facet_grid(~mouseid)+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.text.x=element_blank(),
    axis.text.y=element_text(size=6),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=8),
    legend.justification=c(0,1),
    legend.text = element_text(size=),
    legend.title = element_blank(),
    legend.position = "none"
  )+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(limits=c(0,6.6),breaks=seq(0,6,by=2)) -> p_34

jfc %>%
  filter(box%in%c("03","04")) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  mutate(leadR=lead(R),
         calc=lead(R-loss),
         turnover=pmin(leadR,loss),
         surplus=pmax(leadR-turnover,0),
         deficit=pmax(loss-turnover,0)) %>%
  select(time,mouseid,turnover,surplus,deficit) %>%
  gather(variable,value,-time,-mouseid) %>%
  mutate(variable=factor(variable,levels=c("surplus","deficit","turnover"))) %>%
  ggplot(aes(x=time,y=value,fill=variable,color=variable))+
  geom_bar(stat="identity")+
  scale_fill_manual(values=c(surplus="black",deficit="red",turnover="white"))+
  scale_color_manual(values=c(surplus=NA,deficit=NA,turnover="beige"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(limits=c(0,6e6),labels=aakmisc::scinot)+
  facet_grid(~mouseid)+
  labs(fill="",color="",x="",y="RBC/ \u03BCl blood")+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    axis.text.x=element_blank(),
    axis.text.y=element_text(size=6),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=8),
    legend.text=element_text(size=6),
    legend.justification=c(0,1),
    legend.title = element_blank(),
    legend.position = "none"
  ) -> rep_34

jfc %>%
  filter(
    variable %in% c("Qpn", "Qpw", "Qun", "Qps", "Qus"),
    box%in%c("03","04")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>% 
  gather(variable,value,-c(mouseid, time, box, mouse)) %>%
  mutate(
    variable=factor(variable,levels=c("Qun","Qps","Qpn", "Qpw","Qus"))
  ) %>%
  left_join(peak, by=c("box","mouseid")) %>%
  ggplot()+
  geom_area(aes(x=time, y=value, fill=variable, group=variable))+
  geom_vline(aes(xintercept=time_peakK), color="white")+
  labs(x="day",y="fraction of RBCs")+ 
  scale_fill_manual(
    name="",
    values=c("seagreen","royal blue", "steelblue1",  "midnight blue", "blanchedalmond"), 
    labels=c("unparasitized & destroyed by indiscriminate killing", 
             "parasitized & destroyed by emerging parasites", 
             "parasitized & destroyed by indiscriminate killing",
             "parasitized & destroyed by targeted killing",
             "survived"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  facet_grid(~mouseid)+
  theme(
    strip.text=element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.text = element_text(size=6),
    axis.title.x = element_blank(),
    axis.title = element_text(size=8),
    legend.text = element_text(size=6),
    legend.justification=c(0,1),
    legend.title = element_blank(),
    legend.position = "none"
  )+
  guides(fill = guide_legend(reverse=F, title=NULL)) -> d_34

plot_grid(
  p_34,rep_34,d_34,
  ncol = 1, 
  align="v", 
  axis="lr", 
  rel_heights=c(1,0.7,1))+
  theme(
    plot.background = element_rect(color = "black"),
    axis.title=element_text(size=8)
  ) -> bottom_row

jfc %>%
  filter(
    variable %in% c("lambda_w", "lambda_n","lambda_r","lambda_e", "lambda_u"),
    mouseid%in%c("01-01")
  ) %>%
  left_join(peak, by=c("box","mouseid")) %>% 
  mutate(
    variable=factor(
      variable,
      levels=c("lambda_u","lambda_w","lambda_n","lambda_r","lambda_e"))
  ) %>%
  ggplot(aes(x=time,y=med,group=variable,fill=variable))+
  geom_area(color=NA)+
  geom_rect(fill = "white",xmin = -Inf,xmax = Inf, ymin = -Inf,ymax = Inf)+
  labs(x="day",y= "offspring/ parasite")+
  scale_fill_manual(
    name="",
    values=c("chocolate4", "midnight blue","steelblue1","mistyrose", "light pink"), 
    labels= c("unrealized reproduction", "destroyed by targeted killing", 
              "destroyed by indiscriminate killing",
              "emerged from reticulocytes", "emerged from erythrocytes"))+
  facet_grid(~mouseid)+
  theme(
    panel.border=element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.background = element_blank(), 
    strip.text = element_blank(), 
    axis.text= element_blank(), 
    axis.title = element_blank(), 
    axis.ticks=element_blank(),
    legend.position= c(0,0),
    legend.justification=c(0,0),
    legend.background = element_rect(fill=alpha('white', 0)),
    legend.key.size = unit(1, "lines"),
    legend.text=element_text(size=6, family="Helvetica"),
    legend.title=element_blank()
  ) + 
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(limits=c(0,6.6),breaks=seq(0,6,by=2)) -> p_dummy

jfc %>%
  filter(mouseid=="02-01") %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  mutate(leadR=lead(R),
         calc=lead(R-loss),
         turnover=pmin(leadR,loss),
         surplus=pmax(leadR-turnover,0),
         deficit=pmax(loss-turnover,0)) %>%
  select(time,mouseid,turnover,surplus,deficit) %>%
  gather(variable,value,-time,-mouseid) %>%
  mutate(variable=factor(variable,levels=c("surplus","deficit","turnover"))) %>%
  ggplot(aes(x=time,y=value,fill=variable,color=variable))+
  geom_area()+
  geom_rect(fill = "white", color = "white", xmin = -Inf,xmax = Inf, ymin = -Inf,ymax = Inf)+
  scale_fill_manual(
    values=c("black","red","white"),
    labels=c("surplus","deficit", "turnover")) +
  scale_color_manual(
    values=c("black","red","beige"),
    labels=c("surplus","deficit", "turnover"))+
  theme(panel.border=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.background = element_blank(), 
        strip.text = element_blank(), 
        axis.text= element_blank(), 
        axis.title = element_blank(), 
        axis.ticks=element_blank(),
        legend.position= c(0,0),
        legend.justification=c(0,0),
        legend.background = element_rect(fill=alpha('white', 0)),
        legend.key.size = unit(1, "lines"),
        legend.text=element_text(size=6, family="Helvetica"),
        legend.title=element_blank())+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  scale_y_continuous(limits=c(0,6e6),labels=aakmisc::scinot) -> rep_dummy

jfc %>%
  filter(
    variable %in% c("Qpn", "Qpw", "Qun", "Qps", "Qus"),
    box%in%c("03","04")
  ) %>%
  select(-c(lo,hi)) %>%
  spread(variable,med) %>% 
  gather(variable,value,-c(mouseid, time, box, mouse)) %>%
  mutate(variable=factor(variable,levels=c("Qun","Qps","Qpn", "Qpw","Qus"))) %>%
  left_join(peak, by=c("box","mouseid")) %>%
  ggplot()+
  geom_area(aes(x=time, y=value, fill=variable, group=variable))+
  geom_vline(aes(xintercept=time_peakK), color="white")+
  geom_rect(fill = "white",xmin = -Inf,xmax = Inf, ymin = -Inf,ymax = Inf)+
  labs(x="day",y="fraction of RBCs")+ 
  scale_fill_manual(
    name="",
    values=c("seagreen","royal blue", "steelblue1",  "midnight blue", "blanchedalmond"),
    labels=c("unparasitized & destroyed by indiscriminate killing", 
             "parasitized & destroyed by emerging parasites", 
             "parasitized & destroyed by indiscriminate killing",
             "parasitized & destroyed by targeted killing",
             "survived"))+
  scale_x_continuous(limits=c(0,21),breaks=seq(0,20,by=5))+
  facet_grid(~mouseid)+
  theme(panel.border=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.background = element_blank(), 
        strip.text = element_blank(), 
        axis.text= element_blank(), 
        axis.title = element_blank(), 
        axis.ticks=element_blank(),
        legend.position= c(0.0,0.0),
        legend.justification=c(0,0),
        legend.background = element_rect(fill=alpha('white', 0)),
        legend.key.size = unit(1, "lines"),
        legend.text=element_text(size=6, family="Helvetica"),
        legend.title=element_blank()) + 
  scale_y_continuous(limits=c(0,11.5),breaks=seq(0,10,by=2))+
  guides(fill = guide_legend(reverse=F, title=NULL)) -> d_dummy

plot_grid(
  p_dummy,rep_dummy,d_dummy,
  ncol = 1, 
  align="v", 
  axis="lr", 
  rel_heights=c(1,0.7,1)) -> legends

plot_grid(top_row, NULL,
          nrow = 1, 
          align= "vh",
          rel_widths = c(1,0.33)) -> top_row_blank

plot_grid(bottom_row, legends,
          nrow = 1, 
          align="v", 
          axis="lr", 
          rel_widths=c(1,0.33)) -> bottom_row_l


plot_grid(top_row_blank, bottom_row_l,
          ncol = 1, 
          axis="lr", 
          rel_widths=c(1,1)) -> complete

## code to do the above using cowplot below, adapted from
## https://stackoverflow.com/questions/33114380/centered-x-axis-label-for-muliplot-using-cowplot-package
ggdraw(add_sub(complete, "day",x=0.4, size=10)) -> figS3
figS3
```


-------------------------------

### Figure S4


```{r figS4,purl=FALSE, fig.width=9, fig.height=9.1}
jfc %>%   
  filter(time<22 & box!="05") %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  ddply(.(mouseid), mutate, 
        landmark_a=time[which.max(K)],
        landmark_b=time[which.max(N)]) %>%
  select("mouseid","time","K","R","N","W",  "box", "mouse", "landmark_a","landmark_b") %>%
  mutate(
    landmark_a=case_when(
      mouseid=="01-02" ~ 9L,
      mouseid=="02-02" ~ 7L,
      mouseid=="04-03" ~ 7L,
      TRUE ~ landmark_a),
    landmark_b=case_when(
      ## because i dont want it to appear on the plot but NA gives an error
      mouseid=="01-01" ~ 40L,
      TRUE ~ landmark_b)) -> s4


s4 %>%
  ggplot(aes(x=R,y=W,group=mouseid,color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(s4, time==0), aes(fill=mouseid), shape=21, color="black")+
  geom_point(data=filter(s4, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(s4, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(s4,time==landmark_a), label="a", color="black")+
  geom_text(data=filter(s4,time==landmark_b), label="b", color="black", nudge_x=0.05)+
  geom_text(data = filter(s4, mouseid == "01-01"),
            aes(x= Inf, y=Inf, label="†"), hjust = 1.5, vjust=1.5, color="black")+
  geom_text(data = filter(s4, mouseid %in% c("01-02","02-03")),
            aes(x= Inf, y=Inf, label="V"), hjust = 1.5, vjust=1.5,color="black")+
  labs(y="Targeted killing (W)",x="Reticulocyte supply (R)")+
  theme(
    axis.title.y=element_text(size=10, family="Helvetica"), 
    axis.title.x=element_blank(),
    axis.text= element_text(size=8, family="Helvetica"),
    legend.position="none",
    strip.text=element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm")
  )+
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)),
    limits=c(10^5.3,10^7.9))+
  scale_x_log10(
    limits=c(10^5.1,10^6.75),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)))+
  facet_grid(box~mouse)+
  scale_color_manual(values=rep("goldenrod4",12))+
  scale_fill_manual(values=rep("goldenrod4",12)) -> s_rw

data_labels <- c(
  "01" =  "0.05%",
  "02"  = "0.005%",
  "03" = "0.0005%",
  "04" = "0%"
)

s4 %>%
  ggplot(aes(x=K,y=W,group=mouseid,color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(s4, time==0), aes(fill=mouseid), shape=21, color="black")+
  geom_point(data=filter(s4, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(s4, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(s4,time==landmark_a), label="a", color="black", nudge_x = 0.2)+
  geom_text(data=filter(s4,time==landmark_b), label="b", color="black")+
  geom_text(data = filter(s4, mouseid %in% "01-01"),
            aes(x= Inf, y=Inf, label="†"), hjust = 1.5, vjust=1.5, color="black")+
  geom_text(data = filter(s4, mouseid %in% c("01-02","02-03")),
            aes(x= Inf, y=Inf, label="V"), hjust = 1.5, vjust=1.5,color="black")+
  labs(x="Parasite density (K)",y="Targeted killing (W)", nudge_y=-0.1)+
  theme(
    axis.title.x=element_text(size=10, family="Helvetica"),
    axis.text.x= element_text(size=8, family="Helvetica"),
    axis.title.y =element_blank(),
    axis.text.y =element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.text.x = element_blank(),
    strip.text.y = element_text(size=8, family="Helvetica", angle=0),
    strip.background = element_blank(),
    legend.position="none",
    plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm")
  )+
  scale_y_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)),
    limits=c(10^5.3,10^7.9))+
  scale_x_log10(
    limits=c(1e2,1e7),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)),
    oob=trnc)+
  facet_grid(box~mouse, labeller=labeller(box=data_labels))+
  scale_color_manual(values=rep("goldenrod4",12))+
  scale_fill_manual(values=rep("goldenrod4",12)) -> s_kw

data_labels <- c(
  "01" =  "0.05%",
  "02"  = "0.005%",
  "03" = "0.0005%",
  "04" = "0%")

s4 %>%
  ggplot(aes(x=R,y=N,group=mouseid,color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30))+
  geom_point(data=filter(s4, time==0), aes(fill=mouseid), shape=21, color="black")+
  geom_point(data=filter(s4, time%in%c(1:4,6:9,11:14,16:19)), shape=20, size=1)+
  geom_point(data=filter(s4, time%in%c(5,10,15)), shape=20, size=2)+
  geom_text(data=filter(s4,time==landmark_a), label="a", color="black", nudge_x=-0.05)+
  geom_text(data=filter(s4,time==landmark_b), label="b", color="black",nudge_y=0.1)+
  geom_text(data = filter(s4, mouseid %in% "01-01"),
            aes(x= Inf, y=Inf, label="†"), hjust = 1.5, vjust=1.5, color="black")+
  geom_text(data = filter(s4, mouseid %in% c("01-02","02-03")),
            aes(x= Inf, y=Inf, label="V"), hjust = 1.5, vjust=1.5,color="black")+
  labs(y="Indiscriminate killing (N)",x="Reticulocyte supply (R)")+
  theme(
    axis.title=element_text(size=10, family="Helvetica"),
    axis.text= element_text(size=8, family="Helvetica"),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.text = element_blank(),
    strip.background = element_blank(),
    legend.position="none",
    plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm")
  )+
  scale_y_log10(limits=c(10^5.2,1e7),breaks = trans_breaks("log10", function(x) 10^x),
                labels=trans_format('log10',math_format(10^.x)),oob=trnc)+
  scale_x_log10(limits=c(10^5.1,10^6.75),breaks = trans_breaks("log10", function(x) 10^x),
                labels=trans_format('log10',math_format(10^.x)))+
  facet_grid(box~mouse)+
  scale_color_manual(values=rep("goldenrod4",12))+
  scale_fill_manual(values=rep("goldenrod4",12)) -> s_rn

s4 %>%
  ggplot(aes(x=R,y=N,group=mouseid,color=mouseid))+
  geom_path(arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30),
            colour="goldenrod4")+
  geom_text(data=filter(s4, time%in%c(5,10,15)), aes(label=time), alpha=0.65)+
  geom_point(data=filter(s4, time==0), fill="goldenrod4", colour="goldenrod4", shape=21)+
  geom_text(data=filter(s4,time==landmark_a), label="a", color="black")+
  geom_text(data=filter(s4,time==landmark_b), label="b", color="black")+
  labs(y="Indiscriminate killing (W)",x="Reticulocyte supply (R)")+
  theme(
    axis.title=element_text(size=10, family="Helvetica"),
    axis.text= element_text(size=8, family="Helvetica"),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.text.x = element_blank(),
    strip.text.y = element_text(size=8, family="Helvetica", angle=0),
    strip.background = element_blank(),
    legend.position="none",
    plot.margin = unit(c(0.1, 0.2, 0.1, 0), "cm")
  )+
  scale_y_log10(
    limits=c(10^5.5,1e7),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)),oob=trnc)+
  scale_x_log10(
    limits=c(1e5,1e7),
    breaks = trans_breaks("log10", function(x) 10^x),
    labels=trans_format('log10',math_format(10^.x)))+
  facet_grid(box~mouse, labeller=labeller(box=data_labels)) -> s_rn_labels

s4 %>%
  ggplot(aes(x=R,y=N,group=mouseid, color=mouseid))+
  geom_path(
    arrow=arrow(ends="last", type="closed", length=unit(5,"pt"), angle=30),
    color=NA)+
  theme(
    panel.border=element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.background = element_blank(), 
    strip.text.x = element_blank(), 
    strip.text.y = element_text(size=8, family="Helvetica", angle=0),
    axis.text.y= element_blank(), 
    axis.text.x = element_text(size=8, family="Helvetica", color="white"),
    axis.title.y=element_blank(),
    axis.title.x = element_text(size=10, family="Helvetica",color="white"),
    axis.ticks = element_blank(),
    legend.position="none"
  )+
  scale_color_manual(values=rep("goldenrod4",12))+
  scale_fill_manual(values=rep("goldenrod4",12))+
  scale_y_log10(limits=c(10^5.5,1e7),breaks = trans_breaks("log10", function(x) 10^x),
                labels=trans_format('log10',math_format(10^.x)),oob=trnc)+
  scale_x_log10(limits=c(1e5,1e7),breaks = trans_breaks("log10", function(x) 10^x),
                labels=trans_format('log10',math_format(10^.x)))+
  facet_grid(box~., labeller=labeller(box=data_labels)) -> s_bottomlabels

plot_grid(
  s_rw,  s_kw,
  nrow = 1,
  align ="h",
  axis="tb",
  rel_widths=c(1,1),
  labels=c("A","B")) -> swirl_top

plot_grid(
  s_rn,NULL,s_bottomlabels, NULL,
  nrow = 1,
  axis= "tb",
  rel_widths=c(1,0.03,0.08,0.89),
  labels=c("C",NA,NA,NA)) -> swirl_bottom

plot_grid(
  swirl_top,swirl_bottom,
  ncol=1,
  axis="l",
  rel_heights=c(1,1)) -> whole

whole
```

-------------------------------

## Analyses used in paper

Here, we provide the evidence supporting several statements made in the text of the paper.
In each case, the supporting computations are followed by the quote from the text.

### RBC Limitation

#### Do infected and uninfected mice differ in reticulocyte supply over the first eight days?

```{r purl=FALSE}
jfc %>%
  filter(variable == "R") %>%
  select(-c(lo,hi)) %>%
  spread(variable,med) %>% 
  select(mouseid, time, box, R) %>%
  mutate(
    infection_status=case_when(
      box%in%c("01","02","03","04") ~ "infected",
      box=="05" ~ "uninfected"),
    inoculum=case_when(
      mouseid %in% c("01-02","02-03") ~ "underinoculated",
      !mouseid%in%c("01-02","02-03") ~ "intended inoculum")
  ) -> resupply

resupply %>%
  ggplot(aes(x=time, y=R, color=infection_status, linetype=inoculum, group=mouseid))+
  geom_line()+
  scale_x_continuous(limits=c(0,15), breaks=seq(0,15, by=1))
```

```{r lm,purl=FALSE}
resupply %>%
  filter(time<9) %>%
  group_by(mouseid,infection_status) %>%
  summarize(totretic_8=sum(R)) -> R_data

lm(totretic_8~infection_status, data=R_data) -> model8
summary(model8)
```

> **From the text:** 
> Strikingly, despite incurring huge losses in RBCs, infected mice do not increase the supply of reticulocytes any more than uninfected mice in the first 8 days of infection (*total reticulocyte count, days* 0--8, $F_{`r summary(model8)$fstatistic[2]`, `r summary(model8)$fstatistic[3]`} = `r signif(summary(model8)$fstatistic[1],3)`$, $p = `r summary(model8)$coefficients[2,4] %>% signif(1)`$).


-------------------------------

#### What is the maximum impact of resource limitation?

```{r purl=FALSE}
jfc %>%
  filter(!mouseid%in%c("05-01","05-02","05-03")) %>%
  filter(variable %in% c("lambda_u", "Beta")) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>% 
  select(mouseid, time, Beta, lambda_u) %>%
  group_by(mouseid,Beta) %>%
  summarize(
    peak_lambda_u=max(lambda_u,na.rm=TRUE),
    prop_lambda_u=peak_lambda_u/unique(Beta)
  ) -> df

df %>% pull(prop_lambda_u) %>% mean() %>% signif(2)
df %>% pull(prop_lambda_u) %>% quantile(prob=c(0.1,0.9)) %>% signif(2)
```

> **From the text:** 
> At its maximum, RBC limitation accounts for, on average, a `r 100*df$prop_lambda_u %>% mean() %>% signif(2)`% (interdecile range, `r 100*df$prop_lambda_u %>% quantile(0.1) %>% signif(2)`--`r 100*df$prop_lambda_u %>% quantile(0.9) %>% signif(2)`%) reduction in parasite reproductive potential.

-------------------------------

### The Indiscriminate Response 

#### At the beginning of an infection, what fraction of cells are killed by the indiscriminate response?

```{r purl=FALSE}
jfc %>%
  filter(
    variable %in% c("Qpn", "Qpw", "Qun", "Qps", "Qus"),
    !mouseid%in%c("05-01","05-02","05-03"),
    time<2
  ) %>%
  select(-c(lo,hi)) %>%
  spread(variable,med) %>%
  group_by(mouseid,time) %>%
  mutate(
    totloss=1-Qus, 
    UN_pcent=(Qun/totloss)*100,
    P_pcent=(Qps/totloss)*100
  ) %>%
  ungroup() -> x

x %>% filter(time==1) %>% pull(UN_pcent) %>%  mean() %>% signif(3)

x %>% filter(time==1) %>% pull(UN_pcent) %>% quantile(prob=c(0.1,0.5,0.9)) %>% signif(3)
```

> **From the text:** 
> At the beginning of the infection, almost all of the losses of RBCs can be attributed to the removal of uninfected cells by the indiscriminate response (mean, `r x %>% filter(time==1) %>% pull(UN_pcent) %>% mean() %>% signif(3)`%, interdecile range (`r x %>% filter(time==1) %>% pull(UN_pcent) %>% quantile(prob=0.1) %>% signif(2)`--`r x %>% filter(time==1) %>% pull(UN_pcent) %>% quantile(prob=0.9) %>% signif(2)`%)....

-------------------------------

#### In the pre-peak period, what is the greatest amount of killing that parasites are directly responsible for?

```{r purl=FALSE}
jfc %>%
  filter(!mouseid%in%c("05-01","05-02","05-03")) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  group_by(mouseid) %>%
  mutate(
    time_peakK=time[which.max(K)]
  ) %>%
  filter(time<=time_peakK) %>%
  mutate(
    totloss=1-Qus, 
    UN_pcent=(Qun/totloss)*100,
    totN_pcent=((Qun+Qpn)/totloss)*100,
    P_pcent=(Qps/totloss)*100) -> prepeak

prepeak %>%
  filter(P_pcent==max(P_pcent)) %>%
  ungroup() -> p

p %>%
  select(mouseid, time, P_pcent, UN_pcent, totN_pcent)

p %>% pull(P_pcent) %>% mean() %>% signif(2)
p %>% pull(P_pcent) %>% quantile(prob=c(0.1,0.5,0.9)) %>% signif(2)
```

> **From the text:** 
Even when the contribution of parasite emergence to RBC destruction is at its height, it is never responsible for more than half of the RBC destruction (maximum (across time) proportion of losses attributable to parasite emergence, mean (among mice) `r p %>% pull(P_pcent) %>% mean() %>% signif(2)`%, interdecile range `r p %>% pull(P_pcent) %>% quantile(0.1) %>% signif(2)`%--`r p %>% pull(P_pcent)  %>% quantile(0.9) %>% signif(2)`%)

-------------------------------

#### On the day when parasites are killing the largest fraction of RBCs, does this exceed the destruction of uninfected RBCs by the host?

```{r purl=FALSE}
prepeak %>%
  group_by(mouseid) %>%
  filter(P_pcent==max(P_pcent)) %>%
  mutate(answer=if_else(P_pcent>UN_pcent,"yes","no")) %>%
  select(mouseid,time,UN_pcent,P_pcent,answer) %>%
  print(n=40)
```

-------------------------------

#### On which days is the indiscriminate response not the greatest source of RBC loss?

```{r purl=FALSE}
jfc %>%
  filter(
    ! mouseid %in% c("05-01","05-02","05-03"),
    variable %in% c("Qpn","Qun","Qpw","Qps")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  filter(
    Qpn+Qun < pmin(Qps,Qpw)
  ) %>%
  select(-box,-mouse) %>%
  print(n=40)
```

-------------------------------

#### When does the parasite concentration peak?

```{r purl=FALSE}
jfc %>%
  filter(
    ! mouseid %in% c("05-01","05-02","05-03"),
    variable %in% c("K")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  group_by(mouseid) %>%
  filter(K==max(K)) %>%
  select(mouseid,peakday=time) -> peakdays

peakdays %>% print(n=40)
```

-------------------------------

#### Does it ever happen that the host removes more parasitized cells than unparasitized cells?  If so, by how much?

```{r purl=FALSE}
jfc %>%
  filter(
    ! mouseid %in% c("05-01","05-02","05-03"),
    variable %in% c("Qun","Qpn","Qpw")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  filter(
    Qpw + Qpn > Qun
  ) %>%
  mutate(
    ratio = (Qpw+Qpn)/Qun
  ) %>%
  select(-box,-mouse) %>%
  left_join(peakdays,by="mouseid") %>%
  mutate(time.to.peak=time-peakday) -> x

x %>% print(n=40)
x %>% pull(ratio) %>% mean() %>% signif(2)
x %>% pull(ratio) %>% quantile(prob=c(0.1,0.5,0.9)) %>% signif(3)

```

> **From the text:** 
...indeed there are only a few days, at or around the time of peak parasite density, when the host kills more parasitized cells than uninfected cells.
On these days, `r x %>% pull(ratio) %>% mean() %>% signif(2)` parasitized cells are removed on average for every 1 unparasitized cell (interdecile range `r x %>% pull(ratio) %>% quantile(0.1) %>% signif(2)`--`r x %>% pull(ratio) %>% quantile(0.9) %>% signif(2)`)....


-------------------------------

#### On days when the targeted response kills more parasitized cells than the indiscriminate response, by how much does it do so?

```{r purl=FALSE}
jfc %>%
  filter(variable %in% c("Qpn", "Qpw", "Qun", "Qps","K") & !mouseid%in%c("05-01","05-02","05-03")) %>%
  select(-c(lo,hi)) %>%
  spread(variable,med) %>%
  ddply(.(mouseid,time), mutate, 
        itvsun=Qpw/Qun) %>%
  ddply(.(mouseid), mutate,
        maxitvsun=max(itvsun),
        day_maxitvsun=time[which.max(itvsun)]) %>%
  filter(time==day_maxitvsun) %>%
  arrange(maxitvsun) -> it_vs_un

it_vs_un %>% pull(maxitvsun) %>% quantile(prob=c(0.1,0.5,0.9)) %>% signif(3)
```

-------------------------------

#### In the post-peak phase of infection, on what day is the indiscriminate response responsible for the highest proportion of killed parasitized cells? What is that proportion?

```{r purl=FALSE}
jfc %>%
  filter(
    variable %in% c("Qpn","Qpw"),
    ! mouseid %in% c("01-01","05-01","05-02","05-03")
  ) %>%
  select(-lo,-hi) %>%
  spread(variable,med) %>%
  left_join(peakdays,by="mouseid") %>%
  filter(
    time<22 & time>peakday
  ) %>%
  group_by(mouseid,time) %>%
  summarize(
    prop_pmort_n=100*Qpn/(Qpw+Qpn)
  ) %>%
  group_by(mouseid) %>%
  summarize(
    peak_prop_pmort=max(prop_pmort_n, na.rm=TRUE),
    maxday=time[which.max(prop_pmort_n)]
  ) %>%
  arrange(peak_prop_pmort) -> x

x %>% print(n=40)
x %>% pull(peak_prop_pmort) %>% mean() %>% signif(3)
x %>% pull(peak_prop_pmort) %>% quantile(prob=c(0.1,0.5,0.9)) %>% signif(3)
```

> **From the text:**
> At its peak, this response component is responsible for 
`r x %>% pull(peak_prop_pmort) %>% mean() %>% signif(2)`% 
(interdecile range 
`r x %>% pull(peak_prop_pmort) %>% quantile(0.1) %>% signif(2)`--`r x %>% pull(peak_prop_pmort) %>% quantile(0.9) %>% signif(2)`%) 
of parasite destruction on average.

-------------------------------

## Session information

```{r session-info,cache=FALSE}
bake(file="session_info.rds",{
  sessionInfo()
}) %>%
  print()
```

-----------------

## References
