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”.

Experimental data

The following codes import and pre-process the data. Experimental design and protocols are described in the paper. 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

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.

## 
## Call:
## lm(formula = log(Pd) ~ day + mouseid - 1, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63666 -0.28417 -0.02179  0.33097  0.76375 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## day           1.69193    0.05463  30.973  < 2e-16 ***
## mouseid01-01  7.13040    0.24045  29.655  < 2e-16 ***
## mouseid01-02  1.78189    0.28120   6.337 4.11e-07 ***
## mouseid01-03  6.49483    0.24045  27.011  < 2e-16 ***
## mouseid02-01  6.69657    0.24045  27.850  < 2e-16 ***
## mouseid02-02  5.84940    0.24045  24.327  < 2e-16 ***
## mouseid02-03  0.70050    0.33894   2.067   0.0469 *  
## mouseid03-01  6.52969    0.24045  27.156  < 2e-16 ***
## mouseid03-02  6.09318    0.24045  25.341  < 2e-16 ***
## mouseid03-03  5.90790    0.24045  24.570  < 2e-16 ***
## mouseid04-01  6.58500    0.24045  27.386  < 2e-16 ***
## mouseid04-02  5.86158    0.24045  24.378  < 2e-16 ***
## mouseid04-03  6.32530    0.24045  26.306  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3958 on 32 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.999,  Adjusted R-squared:  0.9986 
## F-statistic:  2387 on 13 and 32 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = log(Pd) ~ box:day + mouseid - 1, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.54582 -0.21023 -0.01221  0.12066  0.49245 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## mouseid01-01  6.67620    0.29204  22.860  < 2e-16 ***
## mouseid01-02  1.23685    0.34621   3.572  0.00126 ** 
## mouseid01-03  6.04063    0.29204  20.684  < 2e-16 ***
## mouseid02-01  6.40958    0.30581  20.959  < 2e-16 ***
## mouseid02-02  5.56240    0.30581  18.189  < 2e-16 ***
## mouseid02-03  0.29871    0.42943   0.696  0.49222    
## mouseid03-01  6.38173    0.27159  23.497  < 2e-16 ***
## mouseid03-02  5.94521    0.27159  21.890  < 2e-16 ***
## mouseid03-03  5.75993    0.27159  21.208  < 2e-16 ***
## mouseid04-01  7.29722    0.27159  26.868  < 2e-16 ***
## mouseid04-02  6.57380    0.27159  24.205  < 2e-16 ***
## mouseid04-03  7.03752    0.27159  25.912  < 2e-16 ***
## box01:day     1.87361    0.09602  19.512  < 2e-16 ***
## box02:day     1.80673    0.10265  17.600  < 2e-16 ***
## box03:day     1.75112    0.08589  20.389  < 2e-16 ***
## box04:day     1.40705    0.08589  16.383 3.38e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3326 on 29 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.9993, Adjusted R-squared:  0.999 
## F-statistic:  2747 on 16 and 29 DF,  p-value: < 2.2e-16

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

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 Pd, RBC, and Retic via a lognormal measurement model.

The following codes implement the model using the pomp package (1, 2). Documentation on this package is available at the pomp website.

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

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 (3) 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.

We display the results of this initial search.

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 (4), as implemented in the panelPomp package, 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.

The above required 1.4 hr on 250 cpus.

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.

The above required 0.68 hr on 250 cpus.

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

The above required 0.65 hr on 250 cpus.

Profile likelihood

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

The above required 6.3 hr on 250 cpus.

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

The above required 0.89 hr on 250 cpus.

Maximum likelihood estimate

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

Smoothing

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

The above required 0.16 hr on 250 cpus.

In the Supporting Information, 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.

Figures

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

Figure 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

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

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

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

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

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

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?

## 
## Call:
## lm(formula = totretic_8 ~ infection_status, data = R_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1179563  -583968  -280658   591709  1487696 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 4686079     242206   19.35 5.77e-11 ***
## infection_statusuninfected  -969367     541589   -1.79   0.0968 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 839000 on 13 degrees of freedom
## Multiple R-squared:  0.1977, Adjusted R-squared:  0.136 
## F-statistic: 3.204 on 1 and 13 DF,  p-value: 0.09678

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_{1, 13} = 3.2\), \(p = 0.1\)).

The Indiscriminate Response

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

## [1] 99.5
##   10%   50%   90% 
##  99.0  99.4 100.0

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, 99.5%, interdecile range (99–100%)….

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

## [1] 33
## 10% 50% 90% 
##  18  36  43

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) 33%, interdecile range 18%–43%)

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

## # A tibble: 12 x 5
## # Groups:   mouseid [12]
##    mouseid  time UN_pcent P_pcent answer
##    <chr>   <int>    <dbl>   <dbl> <chr> 
##  1 01-01       4     42.9    40.7 no    
##  2 01-02       7     61.6    17.9 no    
##  3 01-03       4     46.1    36.5 no    
##  4 02-01       5     28.5    28.7 yes   
##  5 02-02       5     27.4    42.8 yes   
##  6 02-03       8     68.7    15.3 no    
##  7 03-01       4     41.1    42.7 yes   
##  8 03-02       5     22.0    42.3 yes   
##  9 03-03       5     59.2    27.4 no    
## 10 04-01       5     34.2    41.9 yes   
## 11 04-02       5     68.4    21.4 no    
## 12 04-03       5     46.2    35.9 no

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

## # A tibble: 6 x 6
##   mouseid  time     Qpn    Qps    Qpw    Qun
##   <chr>   <int>   <dbl>  <dbl>  <dbl>  <dbl>
## 1 01-01       5 0.0157  0.0752 0.293  0.0533
## 2 01-03       5 0.00708 0.0813 0.121  0.0431
## 3 01-03       6 0.0141  0.0554 0.404  0.0403
## 4 02-02       6 0.0109  0.0685 0.298  0.0396
## 5 03-01       5 0.00440 0.0736 0.109  0.0306
## 6 03-02       5 0.00441 0.0742 0.0580 0.0386

When does the parasite concentration peak?

## # A tibble: 12 x 2
## # Groups:   mouseid [12]
##    mouseid peakday
##    <chr>     <int>
##  1 01-01         5
##  2 01-02         8
##  3 01-03         6
##  4 02-01         6
##  5 02-02         6
##  6 02-03         9
##  7 03-01         6
##  8 03-02         6
##  9 03-03         6
## 10 04-01         6
## 11 04-02         6
## 12 04-03         6

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

## # A tibble: 35 x 8
##    mouseid  time     Qpn    Qpw    Qun ratio peakday time.to.peak
##    <chr>   <int>   <dbl>  <dbl>  <dbl> <dbl>   <int>        <int>
##  1 01-01       5 0.0157  0.293  0.0533  5.79       5            0
##  2 01-01       6 0.0226  0.456  0.0540  8.87       5            1
##  3 01-01       7 0.0409  0.309  0.0896  3.91       5            2
##  4 01-01       8 0.0804  0.413  0.101   4.86       5            3
##  5 01-02       8 0.0111  0.176  0.0860  2.17       8            0
##  6 01-02       9 0.0121  0.229  0.0976  2.47       8            1
##  7 01-03       5 0.00708 0.121  0.0431  2.96       6           -1
##  8 01-03       6 0.0141  0.404  0.0403 10.4        6            0
##  9 01-03       7 0.0225  0.398  0.0696  6.05       6            1
## 10 01-03       8 0.0239  0.296  0.137   2.33       6            2
## 11 02-01       5 0.0171  0.117  0.0894  1.50       6           -1
## 12 02-01       6 0.0325  0.398  0.0719  5.99       6            0
## 13 02-01       7 0.0486  0.498  0.101   5.39       6            1
## 14 02-02       5 0.00464 0.0442 0.0451  1.08       6           -1
## 15 02-02       6 0.0109  0.298  0.0396  7.81       6            0
## 16 02-02       7 0.0131  0.442  0.0504  9.04       6            1
## 17 02-02       8 0.00835 0.164  0.122   1.41       6            2
## 18 02-03       9 0.0162  0.181  0.142   1.39       9            0
## 19 02-03      10 0.0112  0.164  0.134   1.31       9            1
## 20 03-01       5 0.00440 0.109  0.0306  3.71       6           -1
## 21 03-01       6 0.00650 0.367  0.0292 12.8        6            0
## 22 03-01       7 0.00797 0.205  0.0645  3.30       6            1
## 23 03-01       8 0.00555 0.176  0.113   1.60       6            2
## 24 03-02       5 0.00441 0.0580 0.0386  1.62       6           -1
## 25 03-02       6 0.00695 0.362  0.0363 10.2        6            0
## 26 03-02       7 0.00433 0.162  0.0780  2.13       6            1
## 27 03-03       6 0.0623  0.314  0.160   2.35       6            0
## 28 03-03       7 0.0805  0.355  0.166   2.62       6            1
## 29 03-03       8 0.0776  0.224  0.285   1.06       6            2
## 30 04-01       6 0.00862 0.249  0.0566  4.55       6            0
## 31 04-01       7 0.00673 0.135  0.0835  1.70       6            1
## 32 04-02       6 0.0154  0.121  0.125   1.09       6            0
## 33 04-02       7 0.00772 0.197  0.0855  2.39       6            1
## 34 04-03       6 0.0104  0.199  0.0758  2.76       6            0
## 35 04-03       7 0.00641 0.157  0.107   1.54       6            1
## [1] 4
##  10%  50%  90% 
## 1.34 2.62 8.97

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, 4 parasitized cells are removed on average for every 1 unparasitized cell (interdecile range 1.3–9)….

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?

## # A tibble: 11 x 3
##    mouseid peak_prop_pmort maxday
##    <chr>             <dbl>  <int>
##  1 01-02              24.0     12
##  2 04-03              38.0     11
##  3 03-01              42.1     11
##  4 02-03              42.6     15
##  5 04-01              44.4     11
##  6 04-02              56.0     11
##  7 02-02              56.4     11
##  8 03-02              56.6     11
##  9 01-03              67.4     12
## 10 03-03              69.7     11
## 11 02-01              82.3     12
## [1] 52.7
##  10%  50%  90% 
## 38.0 56.0 69.7

From the text: At its peak, this response component is responsible for 53% (interdecile range 38–70%) of parasite destruction on average.

Session information

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8    LC_NUMERIC=C           
##  [3] LC_TIME=en_US.UTF-8     LC_COLLATE=C           
##  [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
##  [7] LC_PAPER=en_US.UTF-8    LC_NAME=C              
##  [9] LC_ADDRESS=C            LC_TELEPHONE=C         
## [11] LC_MEASUREMENT=C        LC_IDENTIFICATION=C    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] doMPI_0.2.2       Rmpi_0.6-9        doParallel_1.0.15
##  [4] aakmisc_0.27-1    doRNG_1.7.1       rngtools_1.4     
##  [7] pkgmaker_0.27     registry_0.5-1    iterators_1.0.12 
## [10] foreach_1.4.7     panelPomp_0.9.1.3 pomp_2.3.0.5     
## [13] stringi_1.4.3     magrittr_1.5      forcats_0.4.0    
## [16] stringr_1.4.0     dplyr_0.8.3       purrr_0.3.2      
## [19] readr_1.3.1       tidyr_1.0.0       tibble_2.1.3     
## [22] ggplot2_3.2.1     tidyverse_1.2.1   plyr_1.8.4       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2        lubridate_1.7.4   mvtnorm_1.0-11   
##  [4] lattice_0.20-38   utf8_1.1.4        assertthat_0.2.1 
##  [7] zeallot_0.1.0     digest_0.6.21     R6_2.4.0         
## [10] cellranger_1.1.0  backports_1.1.5   coda_0.19-3      
## [13] httr_1.4.1        pillar_1.4.2      rlang_0.4.0      
## [16] lazyeval_0.2.2    RPostgreSQL_0.6-2 curl_4.2         
## [19] readxl_1.3.1      rstudioapi_0.10   munsell_0.5.0    
## [22] broom_0.5.2       compiler_3.6.1    modelr_0.1.5     
## [25] pkgconfig_2.0.3   tidyselect_0.2.5  codetools_0.2-16 
## [28] fansi_0.4.0       crayon_1.3.4      withr_2.1.2      
## [31] grid_3.6.1        nlme_3.1-141      jsonlite_1.6     
## [34] xtable_1.8-4      gtable_0.3.0      lifecycle_0.1.0  
## [37] DBI_1.0.0         scales_1.0.0      bibtex_0.4.2     
## [40] cli_1.1.0         reshape2_1.4.3    xml2_1.2.2       
## [43] ellipsis_0.3.0    generics_0.0.2    vctrs_0.2.0      
## [46] deSolve_1.24      tools_3.6.1       glue_1.3.1       
## [49] hms_0.5.1         colorspace_1.4-1  rvest_0.3.4      
## [52] haven_2.1.1

References

1. King AA, et al. (2019) pomp: Statistical inference for partially observed Markov processes Available at: https://kingaa.github.io/pomp/.

2. King AA, Nguyen D, Ionides EL (2016) Statistical inference for partially observed Markov processes via the R package pomp. Journal of Statistical Software 69(12):1–43.

3. Ionides EL, Nguyen D, Atchadé Y, Stoev S, King AA (2015) Inference for dynamic and latent variable models via iterated, perturbed Bayes maps. Proceedings of the National Academy of Sciences of the United States of America 112(3):719–724.

4. Bretó C, Ionides EL, King AA (2019) Panel data analysis via mechanistic models. Journal of the American Statistical Association in press.