This repository contains code used in the following paper.

Hélène Cecilia, Roosmarie Vriens, Jeroen Kortekaas, Paul J. Wichgers Schreur, Mariken de Wit, Raphaëlle Métras, Pauline Ezanno, Quirine A. ten Bosch (2021)
Heterogeneity of Rift Valley fever virus transmission potential across livestock hosts, quantified through a model-based analysis of host viral load and vector infection.
bioRxiv (https://doi.org/10.1101/2021.10.22.465395)

See README for more information on how to explore and use this code base.

Data/speciesData/lamb_dead_edited.csv

-"xvals","outcomeRNA","outcomeTCID","host_type","outcome_survival","id","study","inoculum","time_of_death"
-0,1.7,1.55,"lamb","die",2462,"main",52.6,FALSE
-1,4.34,1.55,"lamb","die",2462,"main",52.6,FALSE
-2,9.3,6.45,"lamb","die",2462,"main",52.6,FALSE
-3,9.91,6.21,"lamb","die",2462,"main",52.6,FALSE
-4,9.85,5.95,"lamb","die",2462,"main",52.6,TRUE
-0,1.7,1.55,"lamb","die",2464,"main",52.6,FALSE
-1,4.23,1.55,"lamb","die",2464,"main",52.6,FALSE
-2,8.78,3.88,"lamb","die",2464,"main",52.6,FALSE
-3,9.71,6.45,"lamb","die",2464,"main",52.6,FALSE
-4,8.98,4.82,"lamb","die",2464,"main",52.6,FALSE
-5,7.83,4.58,"lamb","die",2464,"main",52.6,FALSE
-6,7.21,1.79,"lamb","die",2464,"main",52.6,FALSE
-7,6.46,1.55,"lamb","die",2464,"main",52.6,TRUE
-0,1.7,1.55,"lamb","die",2465,"main",52.6,FALSE
-1,4.33,1.55,"lamb","die",2465,"main",52.6,FALSE
-2,8.76,5.75,"lamb","die",2465,"main",52.6,FALSE
-3,9.48,6.21,"lamb","die",2465,"main",52.6,TRUE
-0,1.7,1.55,"lamb","die",2467,"main",52.6,FALSE
-1,4.41,1.55,"lamb","die",2467,"main",52.6,FALSE
-2,8.95,5.98,"lamb","die",2467,"main",52.6,FALSE
-3,9.8,7.15,"lamb","die",2467,"main",52.6,FALSE
-4,9.57,5.51,"lamb","die",2467,"main",52.6,FALSE
-5,8.4,4.82,"lamb","die",2467,"main",52.6,FALSE
-6,8.79,4.12,"lamb","die",2467,"main",52.6,FALSE
-7,8.38,3.18,"lamb","die",2467,"main",52.6,TRUE
-0,1.7,1.55,"lamb","die",2468,"main",52.6,FALSE
-1,4.34,1.55,"lamb","die",2468,"main",52.6,FALSE
-2,8.6,6.68,"lamb","die",2468,"main",52.6,FALSE
-3,9.56,6.45,"lamb","die",2468,"main",52.6,FALSE
-4,9.27,5.98,"lamb","die",2468,"main",52.6,FALSE
-5,8.43,5.28,"lamb","die",2468,"main",52.6,FALSE
-6,8.12,3.42,"lamb","die",2468,"main",52.6,TRUE
-0,1.7,1.55,"lamb","die",2469,"main",52.6,FALSE
-1,4.95,1.55,"lamb","die",2469,"main",52.6,FALSE
-2,9.44,6.21,"lamb","die",2469,"main",52.6,FALSE
-3,9.77,7.15,"lamb","die",2469,"main",52.6,FALSE
-4,9.06,5.25,"lamb","die",2469,"main",52.6,FALSE
-5,8.47,5.28,"lamb","die",2469,"main",52.6,TRUE
-0,1.7,1.55,"lamb","die",1593,"supp",52.6,FALSE
-1,4.43,1.55,"lamb","die",1593,"supp",52.6,FALSE
-2,8.76,5.21,"lamb","die",1593,"supp",52.6,FALSE
-3,9.54,6.41,"lamb","die",1593,"supp",52.6,FALSE
-4,8.71,4.21,"lamb","die",1593,"supp",52.6,FALSE
-5,6.69,2.81,"lamb","die",1593,"supp",52.6,FALSE
-6,6.04,4.21,"lamb","die",1593,"supp",52.6,FALSE
-7,4.91,1.55,"lamb","die",1593,"supp",52.6,TRUE
-0,1.7,1.55,"lamb","die",1594,"supp",52.6,FALSE
-1,4.73,1.55,"lamb","die",1594,"supp",52.6,FALSE
-2,9.52,5.82,"lamb","die",1594,"supp",52.6,FALSE
-3,9.63,5.62,"lamb","die",1594,"supp",52.6,TRUE
-0,1.7,1.55,"lamb","die",1603,"supp",52.6,FALSE
-1,5.02,1.55,"lamb","die",1603,"supp",52.6,FALSE
-2,8.94,5.82,"lamb","die",1603,"supp",52.6,FALSE
-3,9.35,6.02,"lamb","die",1603,"supp",52.6,FALSE
-4,8.3,4.21,"lamb","die",1603,"supp",52.6,FALSE
-5,6.73,4.01,"lamb","die",1603,"supp",52.6,FALSE
-6,7.81,1.55,"lamb","die",1603,"supp",52.6,TRUE
-0,1.7,1.55,"lamb","die",1614,"supp",52.6,FALSE
-1,1.7,1.55,"lamb","die",1614,"supp",52.6,FALSE
-2,7.42,5.21,"lamb","die",1614,"supp",52.6,FALSE
-3,9.82,6.22,"lamb","die",1614,"supp",52.6,FALSE
-4,8.9,5.42,"lamb","die",1614,"supp",52.6,FALSE
-5,8.58,4.81,"lamb","die",1614,"supp",52.6,FALSE
-6,7.99,3.41,"lamb","die",1614,"supp",52.6,TRUE

manuscript_draft.Rmd

----
-title: "RVFV Within-host model"
-author: "Hélène Cecilia"
-date: "31 août 2020"
-output: html_document
-bibliography: RVF_within_host_collab.bib
----
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE)
-```
-### Introduction
-
-***
-**_NOTE:_** If we end up thinking that our results have a broad interest outside of RVF (e.g for other VBDs or other segmented genomes), include a more general first paragraph and present RVFV as a good case study?
-
-***
-
-1. Present Rift Valley fever
-  + impact worldwide (population-level aspects)
-  + type of virus, segmented genome, known and unknown aspects of immune response
-  + the need for a vaccine, limits of existing ones : is it directly relevant to our question? Else, discussion.
-2. The use of within-host models
-  + very few (only one? [@tuncer_2016]) within-host model for RVFV
-  + what has been shown for other viruses
-  + how it can complement between-host models
-3. The need for an RVF within-host model
-  + understand the processes driving RVFV viral load dynamics
-  + compare different experimental setups and see if "one model fits all"
-  
-***
-**_NOTE:_** Due to the wide knowledge gap in terms of RVFV within-host models, there are many research questions we could try to answer with the data available. Moving forward, we need to make sure of what is the focus of our study, to narrow the list of models we test.
-
-* Test if there's a difference between intraveinous infection and mosquito bites?
-* Estimate the inoculum from mosquito bites? [@li_handel_2014]
-* Assess the role of co-infections in shaping the infectious viral load? [@fonville_2015, @koelle_2019] Keeping in mind the opposing effects of complementary particles [@jacobs_2019] and defective interfering [@frensing_2013, @mapder_2019]
-* Explore different mode of action of innate/adaptive immune response? [@handel_2010, @cao_2015, @canini_carrat_2011, @hernandez-vargas_2014, @ben-shachar_koelle_2015]
-* Assess the role of temperature (fever) in viral decay? [@handel_2013]
-* Evaluate the added value of "dual-measurement" experiments for within-host model? [@petrie_2013]
-
-***
-
-<!-- Remark : In Li Y, Handel A. 2014 Modeling inoculum dose dependent patterns of acute virus infections. J. Theor. Biol. 347, 63–73. (doi:10.1016/j.jtbi.2014.01.008), they show a clear impact of the inoculum on the subsequent dynamics of the viral load. Should we estimate the viral load for experiments starting with a mosquito bite? with the number of positive bites, we could end up with an average inoculum per mosquito? See how they fit the viral loads (with phenomenological equation) and then see which type of immune component can recreate the correlation pattern with inoculum. Least square fitting -->
-<!-- Holder and Beauchemin 2011. Test of different models (ODE exponential time spent in state, Dirac (fixed delay), normal and log-normal delays) : show that normal and log-normal capture best the features, but Dirac can be a good enough approximation. Also shows that a simple four parameters equation (the one used in Li and Handel 2014) can fit well the data. This indicates that model selection based on viral titers only is not a good solution, because you end up with differences estimates of the key infection parameters. -->
-<!-- Least square fitting -->
-<!-- Parameter identifiability : cf Bellu et al. 2007 and Miao et al. 2011 -->
-<!-- Influence of temperature : Handel 2013 : influenza, relevant question because persistence in the environment + within the host, very different temperature ranges. Maybe not relevant for RVF, unless direct transmission (or temperature within mosquitoes) involves different temperatures. They include it as a dependance in the decay rate of the virus : c(T) = alpha.exp(gamma.T) -->
-
-### Material & Methods
-```{r, echo = FALSE}
-eq1 <- noquote(paste(expression("\\frac{dU}{dt} = -\\beta.U.V_{inf} \\\\ \\frac{dI}{dt} = \\beta.U.V_{inf} - \\delta.I \\\\ \\frac{dV_{inf}}{dt} = \\omega.p.I - \\kappa.V_{inf} - \\beta.U.V_{inf} \\\\ \\frac{dV_{non.inf}}{dt} = \\omega.(1-p).I - \\kappa.V_{non.inf}")))
-eq1_alt <- noquote(paste(expression("\\frac{dU}{dt} = -\\beta.U.V_{inf} \\\\ \\frac{dI}{dt} = \\beta.U.V_{inf} - \\delta.I \\\\ \\frac{dV_{inf}}{dt} = \\omega.p.I - \\kappa.V_{inf} - \\beta.U.V_{inf} \\\\ \\frac{dV_{non.inf}}{dt} = \\omega.(1-p).I - \\kappa.V_{non.inf} \\\\ \\frac{dV_{non.inf.M}}{dt} = \\epsilon.V_{non.inf}")))
-```
-
-1. Target-cell limited model
-  + 5 parameters
-  + infection resolves because there is no cell to infect anymore, and infected cells die and stop producing new virus.
-  + only $V_{inf}$ can infect new cells
-
-$$
-\begin{cases} `r eq1` \end{cases}
-$$
-* The PCR dataset used comes from the amplification of the M segment, and we want to fit $V_{inf}$ to TCID$_{50}$ and $V_{inf} + V_{non.inf}$ to rt-qPCR (RNA copies). Therefore the model does not represent the true level of non infectious particles, only those containing the M segment, to be consistent with what is measurable.
-* Use of the segment data : $\omega =$ number of progeny virion containing the M segment produced per infected cell per day. As the progeny dataset is the viral production 8h post-infection (when the latent phase just ended) and the stock dataset is several days post infection (with no knowledge of the number of cells contributing to the production), we cannot infer a value for $\omega$. However, we can set $p$ as the average fraction $\frac{SML}{virionContainsM}$ over the 25 cells present in the progeny dataset, which is 12% (median 9.1%). 
-
-**Update** We'll differentiate between the whole stock of non-infectious particles and the observable ones by including a new state:
-$$
-\begin{cases} `r eq1_alt` \end{cases}
-$$
-where $\epsilon$ will be the mean proportion of non infectious particles containing the M segment in the progeny dataset (fixed). RNA measures will be fitted to $V_{inf} + V_{non.inf.M}$. We need to see if non infectious particles should include the empty ones (with empty particles, $\epsilon = 14\%$, without $34\%$). It depends on the structure of the self-limited infection model, and whether these empty particles are assumed to play the same kind of role as the incomplete ones.
-
-* @handel_2013 introduce temperature-dependence (for influenza) by setting the decay rate $\kappa = \alpha.e^{\gamma.T}$. But he had data on virus persistence in the environment for different temperatures which allowed him to tranpose it within the host.
-* @petrie_2013, @schulze-horsel_2009, @iwami_2012 and @pinilla_2012 have interesting models with infectious and non-infectious virions, fitted to RNA and TCID$_{50}$, where they consider the loss of infectivity turning infectious particles into non-infectious ones. I should have a closer look at these. 
-
-**Update** The differenciation between $V_{non.inf}$ types was not done in @iwami_2012 but it might be because they were not dealing with a segmented genome. They had data to estimate the loss of infecitvity rate witout inference (0.039 day$^{-1}$). @petrie_2013, @petrie_2015 do rescale total loads to PCR. They mention a time-dependency of the ratio between the two measures, as well as the overall virion production. In @simon_2016, total $V_{RNA}$ and infectious $V_{TCID}$ are both produced by infected cells at independent rates, the fact that one is a proportion of the other is not explicit. Same for @pinilla_2012, who also had data to estimate the loss of infecitvity rate witout inference  (0.13h$^{-1}$).
-
-2. Introducing an eclipse phase
-
-It has been evidenced *in vitro* that the production of progeny virions starts around 8h after a cell has been infected, which could justify the introduction of a latently infected stage $E$. We can classically model the time spent in this state with an exponential distribution ($g = 1/8$ hour$^{-1}$). In this case, the mean time spent in the $E$ state is indeed 8 hours but the median is actually $8 \times ln2 = 5.5$:
-
-```{r, echo = FALSE}
-eq2 <- noquote(paste(expression("\\frac{dU}{dt} = -\\beta.U.V_{inf} \\\\ \\frac{dE}{dt} = \\beta.U.V_{inf} - g.E \\\\ \\frac{dI}{dt} = g.E - \\delta.I")))
-```
-
-$$
-\begin{cases} `r eq2` \end{cases}
-$$
-
-@holder_beauchemin_2011 advise a normal or log-normal distribution in order to fully capture the feature of the dynamic at the early stages. This can be approached by an Erlang distribution, which consists of several $E$ states, as shown in @krylova_earn_2013 (I can plot a few probability density function so we can choose a sensible $n$ and therefore not add a parameter to estimate):
-
-```{r, echo = FALSE}
-eq3 <- noquote(paste(expression("\\frac{dU}{dt} = -\\beta.U.V_{inf} \\\\ \\frac{dE_1}{dt} = \\beta.U.V_{inf} - n.g.E_1 \\\\ \\frac{dE_2}{dt} = n.g.E_1 - n.g.E_2 \\\\
-                                ... \\\\ \\frac{dE_n}{dt} = n.g.E_{n-1} - n.g.E_n \\\\ \\frac{dI}{dt} = g.E_n - \\delta.I")))
-```
-
-$$
-\begin{cases} `r eq3` \end{cases}
-$$
-
-3. Add a conversion factor from $V_{inf}$ to TCID$_{50}$ or PFU (plaque forming units)
-
-* In @handel_2010 (fitting EID titers), in the free virus compartment : $-\gamma.\beta.U.V$
-
-"The factor $\gamma$ in the equation for virions is needed for the conversion from infectious virions, the units of V in our model, to EID or PFU as measured by the experimental data. A value of $\gamma$ = 1 assumes that one virion corresponds to one EID or PFU, $\gamma$ < 1 allows for the possibility that more than one virion is required to produce one EID or PFU." $\gamma$ is fitted by the model, it ranges from 1 to 7.8.$10^{-8}$ depending on the model.
-
-* In @canini_2016 (fitting TCID$_{50}$ titers), a parameter $\sigma$ is used (in the same place as $\gamma$), with a slightly different definition:
-
-"A conversion factor $\sigma$ is needed for the conversion from infectious virus, V in our model, to TCID$_{50}$/mL as measured by the experimental data. Here we set $\sigma$ = 0.69 TCID$_{50}$/mL.cell$^{-1}$ which is consistent with 1 mL virus stock having half the number of plaque forming units (PFUs) as the TCID$_{50}$."
-
-4. What type of target cells are we modelling?
-
-We first intended to set our parameter values for the liver/hepatocytes, saying that it was the preferential target organ and that we could consider that the virus quantites measured in the serum were the result of what was happening in the liver, but it seemed a bit odd. Macrophages and monocytes (blood cells) are actually preferential target cells [@gommet_2011], so maybe we should reconsider? Target cell parameters in @clapham_2014 were calibrated on macrophages and monocytes as well.
-
-<!-- *** -->
-<!-- **_Modelling hypothesis (yes/no):_** -->
-
-<!-- * Neglect target cell renewal : sensible for acute infections -->
-<!-- * Neglect loss of free virus entering host cells : might not be true especially at the beginning of infection -->
-<!-- * Account for a conversion TCID to PFU [@canini_2016, @handel_2010] -->
-<!-- * Account for missing particles in rt-qPCR measures -->
-<!-- * Replace exponential by Dirac delta distribution for latent state [@holder_beauchemin_2011] -->
-
-<!-- *** -->
-
-***
-**_Fitting method:_**
-
-* Common errors when doing inference in within-host viral models [@nguyen_2016]
-* Issue when estimating key infection parameters with only virus data [@holder_beauchemin_2011]
-* Many papers used least-square methods to fit their models. I'm tempted to have a go at the **nloptr** package. Could we agree on one first model that we both test (Quirine with MCMC and Helene with least-square) ?
-
-***
-
-<!-- Use segment data to fix $p$ : directly from the data (cite the paper) or from our mini-model (put in Suppl.)? -->
-<!-- Do we define $p$ as $\frac{complete}{total}$ or $\frac{contains M \& complete}{contains M}$? Same goes for $\omega$. -->
-
-<!-- Should we (and if yes, how) account for viral clearance involving the degradation of infectious particles into non infectious ones? -->
-
-<!-- Changement probable : macrophages (origine sanguine) sont les cellules cibles privilégiées [@gommet_2011], revoir les chiffres sur les target cell density, lifespan (c'était également ce qu'avait consideré @clapham_2014) -->
-<!-- Si on veut rester sur les hepathocytes, probablement prévoir un coefficient entre la production de virus dans l'organe et ce qui est retrouvé dans le sérum. Pb : macrophages font aussi partie de l'immunité ??? Confusion des genres -->
-
-| **Parameter** | **Definition** | **Value** | **Source** |
-|:----:|--------------|-------------|-------------------|
-|$S_0$| Target cell density|$3.10^7$ cells.mL$^{-1}$| Liver : $9.08.10^{10}$ cells. Blood : 3000 mL. Ref?|
-|$\gamma^{-1}$| Healthy target cell lifespan | 45 days? | Ref? |
-|$A$| Target cell production |$\gamma.S_0$ (cells.mL$^{-1}$.day$^{-1}$)| Constant cell population|
-|$\sigma^{-1}$| Eclipse phase duration | 8 hours | Ref? |
-|$\delta^{-1}$| Infected cell lifespan | 4 days? | Ref? |
-|$\alpha$| Rate of removal of infected cells by immune response | 0.001 (immunity compotent.mL$^{-1}$)$^{-1}$.day$^{-1}$| Arbitrary (Clapham et al. 2014) |
-|$\beta$| Infection rate | to be estimated (viral particle.mL$^{-1}$)$^{-1}$.day$^{-1}$ | | 
-|$\omega$| Free virus production factor | ? (correlated with $\beta$, so should be fixed) | |
-|$\kappa$| Clearance rate of viral particles | to be estimated (day$^{-1}$) | |
-|$p$| Proportion of infectious progeny virions | to be estimated | |
-|$\eta$| Proliferation rate of immune response | to be estimated | |
-|$V_0^+$| Viral inoculum | 1 TCID$_{50}$.mL$^{-1}$ (include in sensitivity analysis)| For intravenous infection, 10$^5$ TCID$_{50}$ is injected, for 3000 mL of blood, which corresponds to $\simeq$ 30 viral particles per mL. This is consistent with the measure at $t_0$ being below the limit of detection (10$^{1.55}$)|
-|$Z_0$| Initial immune response | ? (should be small, probably smaller than 1.mL$^{-1}$) | Clapham et al. estimated 0.37|
-
-```{r, include=FALSE}
-    # gamma <- ifelse("gamma" %in% names(parameters), parameters[["gamma"]],1/(45*24)) #------------- susceptible cells live 45 days
-    # A <- ifelse("A" %in% names(parameters), parameters[["A"]], gamma * 3e+7) #--------------------- S0 = 3e+7 liver cells/mL of blood
-    # sigma <- ifelse("sigma" %in% names(parameters), parameters[["sigma"]], 1/8) #------------------ latent state during 8 hours
-    # delta <- ifelse("delta" %in% names(parameters), parameters[["delta"]], 1/(4*24)) #------------- infected cells last 4 days
-    # alpha <- ifelse("alpha" %in% names(parameters), parameters[["alpha"]], 0.001/24) #------------- rate of removal of infected cells by immune response
-    # omega <- ifelse("omega" %in% names(parameters), parameters[["omega"]], 1.4e+4/24) #------------ production of virions by infected cells
-    # 
-    # beta <- ifelse("beta" %in% names(parameters), parameters[["beta"]], 2e-7) #-------------------- rate of virions entering cells
-    # kappa.good <- ifelse("kappa.good" %in% names(parameters), parameters[["kappa.good"]], 1/6.5) #- natural clearance of good particles
-    # kappa.bad <- ifelse("kappa.bad" %in% names(parameters), parameters[["kappa.bad"]], 1/6.5) #---- natural clearance of bad particles
-    # p.inf <- ifelse("p.inf" %in% names(parameters), parameters[["p.inf"]], 3e-4) #----------------- proportion of complete particles among progeny virions
-    # eta <- ifelse("eta" %in% names(parameters), parameters[["eta"]], 5e-7) #----------------------- proliferation of immune response
-```
-
-### Discussion
-<!-- @handel_rohani_2015 (citations, do not use as is): -->
-<!-- The processes of pathogen invasion of the host, its sub-sequent spread, interplay with host immunity and the consequent pathogenesis impacts are central to understanding population-level transmission and mitigating the morbidity and mortality of infected hosts. Example of neuraminidase inhibitors. If instead there is a general theory on scaling from individual-level measurements of virus load and symptom severity to between-host transmission fitness, we could use more readily available within-host data to make quantitative predictions about outbreak mitigation. Progress towards a predictive multi-scale framework will require a more precise,quantitativeunderstanding of how infection dynamics, pathogen load, target cell depletion, immunology, symptomatology and other clinical features combine to shape pathogen transmission fitness at the population level.  -->
-<!-- Figure 6. Experimental set-up to determine infectiousness as function of within-host infection dynamics. The infected host is repeatedly sampled to determine as many infection-related quantities as possible (e.g. pathogen load, various immune response components, symptoms). In addition, sets of susceptible hosts are exposed to the infected host at various intervals to determine transmission. This could allow one to obtain a quantitative mapping between quantities such as pathogen load and symptoms and transmission potential [85]. It might even be possible to set up the experiment in such a way that potential contact behaviour changes in the infected host or the contacts can be measured.   -->
-<!-- If, for some pathogen, increased transmission is mainly a function of host behaviour,a different strategy is called for compared with a situation where increased transmission is mainly associated with specific types of symptoms or high pathogen load. Beyond intervention strategies, linking the within-host and between-host scales will be important in obtaining a more complete and predictive understanding of host–pathogen ecology and evolution. -->
-<!-- We have focused only on the question of scaling from within-host to between-host levels from a single pathogen genotype point of view, without considering changes in genotype, i.e. we did not consider explicitly evolutionary processes. This ignores the possibility of competitive dynamics between genotypes, which can be especially important for pathogens with high mutation rates and those that lead to long-term infection. -->
-<!-- We have also not discussed how to include a distinct transmission stage in the process of scaling from individual infection to population-level transmission (vector borne diseases is not mentionned but I think it relates to this point) -->
-<!-- Idée: aller voir si virus peak est négativement correlé à durée infection -->
-
-### References
\ No newline at end of file