1. Ecology
  2. Epidemiology and Global Health

Accelerated viral dynamics in bat cell lines, with implications for zoonotic emergence

  1. Cara E Brook  Is a corresponding author
  2. Mike Boots
  3. Kartik Chandran
  4. Andrew P Dobson
  5. Christian Drosten
  6. Andrea L Graham
  7. Bryan T Grenfell
  8. Marcel A Müller
  9. Melinda Ng
  10. Lin-Fa Wang
  11. Anieke van Leeuwen
  1. Department of Integrative Biology, University of California, Berkeley, United States
  2. Department of Ecology and Evolutionary Biology, Princeton University, United States
  3. Department of Microbiology and Immunology, Albert Einstein College of Medicine, United States
  4. Institute of Virology, Charité-Universitätsmedizin Berlin, corporate member of Freie Universität Berlin, Humboldt-Universität zu Berlin, and Berlin Institute of Health, Germany
  5. Fogarty International Center, National Institutes of Health, United States
  6. Martsinovsky Institute of Medical Parasitology, Tropical and Vector Borne Diseases, Sechenov University, Russian Federation
  7. Emerging Infectious Diseases Program, Duke-National University of Singapore Medical School, Singapore
  8. Royal Netherlands Institute for Sea Research, Department of Coastal Systems, and Utrecht University, Netherlands
https://doi.org/10.7554/eLife.48401
Share this article
Cite this article
  1. Cara E Brook
  2. Mike Boots
  3. Kartik Chandran
  4. Andrew P Dobson
  5. Christian Drosten
  6. Andrea L Graham
  7. Bryan T Grenfell
  8. Marcel A Müller
  9. Melinda Ng
  10. Lin-Fa Wang
  11. Anieke van Leeuwen
(2020)
Accelerated viral dynamics in bat cell lines, with implications for zoonotic emergence
eLife 9:e48401.
https://doi.org/10.7554/eLife.48401
  2. Download BibTeX
  3. Download .RIS

Abstract

Bats host virulent zoonotic viruses without experiencing disease. A mechanistic understanding of the impact of bats’ virus hosting capacities, including uniquely constitutive immune pathways, on cellular-scale viral dynamics is needed to elucidate zoonotic emergence. We carried out virus infectivity assays on bat cell lines expressing induced and constitutive immune phenotypes, then developed a theoretical model of our in vitro system, which we fit to empirical data. Best fit models recapitulated expected immune phenotypes for representative cell lines, supporting robust antiviral defenses in bat cells that correlated with higher estimates for within-host viral propagation rates. In general, heightened immune responses limit pathogen-induced cellular morbidity, which can facilitate the establishment of rapidly-propagating persistent infections within-host. Rapidly-transmitting viruses that have evolved with bat immune systems will likely cause enhanced virulence following emergence into secondary hosts with immune systems that diverge from those unique to bats.

eLife digest

Bats can carry viruses that are deadly to other mammals without themselves showing serious symptoms. In fact, bats are natural reservoirs for viruses that have some of the highest fatality rates of any viruses that people acquire from wild animals – including rabies, Ebola and the SARS coronavirus.

Bats have a suite of antiviral defenses that keep the amount of virus in check. For example, some bats have an antiviral immune response called the interferon pathway perpetually switched on. In most other mammals, having such a hyper-vigilant immune response would cause harmful inflammation. Bats, however, have adapted anti-inflammatory traits that protect them from such harm, include the loss of certain genes that normally promote inflammation. However, no one has previously explored how these unique antiviral defenses of bats impact the viruses themselves.

Now, Brook et al. have studied this exact question using bat cells grown in the laboratory. The experiments made use of cells from one bat species – the black flying fox – in which the interferon pathway is always on, and another – the Egyptian fruit bat – in which this pathway is only activated during an infection. The bat cells were infected with three different viruses, and then Brook et al. observed how the interferon pathway helped keep the infections in check, before creating a computer model of this response.

The experiments and model helped reveal that the bats’ defenses may have a potential downside for other animals, including humans. In both bat species, the strongest antiviral responses were countered by the virus spreading more quickly from cell to cell. This suggests that bat immune defenses may drive the evolution of faster transmitting viruses, and while bats are well protected from the harmful effects of their own prolific viruses, other creatures like humans are not.

The findings may help to explain why bats are often the source for viruses that are deadly in humans. Learning more about bats' antiviral defenses and how they drive virus evolution may help scientists develop better ways to predict, prevent or limit the spread of viruses from bats to humans. More studies are needed in bats to help these efforts. In the meantime, the experiments highlight the importance of warning people to avoid direct contact with wild bats.

Introduction

Bats have received much attention in recent years for their role as reservoir hosts for emerging viral zoonoses, including rabies and related lyssaviruses, Hendra and Nipah henipaviruses, Ebola and Marburg filoviruses, and SARS coronavirus (Calisher et al., 2006; Wang and Anderson, 2019). In most non-Chiropteran mammals, henipaviruses, filoviruses, and coronaviruses induce substantial morbidity and mortality, display short durations of infection, and elicit robust, long-term immunity in hosts surviving infection (Nicholls et al., 2003; Hooper et al., 2001; Mahanty and Bray, 2004). Bats, by contrast, demonstrate no obvious disease symptoms upon infection with pathogens that are highly virulent in non-volant mammals (Schountz et al., 2017) but may, instead, support viruses as long-term persistent infections, rather than transient, immunizing pathologies (Plowright et al., 2016).

Recent research advances are beginning to shed light on the molecular mechanisms by which bats avoid pathology from these otherwise virulent pathogens (Brook and Dobson, 2015). Bats leverage a suite of species-specific mechanisms to limit viral load, which include host receptor sequence incompatibilities for some bat-virus combinations (Ng et al., 2015; Takadate et al., 2020) and constitutive expression of the antiviral cytokine, IFN-α, for others (Zhou et al., 2016). Typically, the presence of viral RNA or DNA in the cytoplasm of mammalian cells will induce secretion of type I interferon proteins (IFN-α and IFN-β), which promote expression and translation of interferon-stimulated genes (ISGs) in neighboring cells and render them effectively antiviral (Stetson and Medzhitov, 2006). In some bat cells, the transcriptomic blueprints for this IFN response are expressed constitutively, even in the absence of stimulation by viral RNA or DNA (Zhou et al., 2016). In non-flying mammals, constitutive IFN expression would likely elicit widespread inflammation and concomitant immunopathology upon viral infection, but bats support unique adaptations to combat inflammation (Zhang et al., 2013; Ahn et al., 2019; Xie et al., 2018; Pavlovich et al., 2018) that may have evolved to mitigate metabolic damage induced during flight (Kacprzyk et al., 2017). The extent to which constitutive IFN-α expression signifies constitutive antiviral defense in the form of functional IFN-α protein remains unresolved. In bat cells constitutively expressing IFN-α, some protein-stimulated, downstream ISGs appear to be also constitutively expressed, but additional ISG induction is nonetheless possible following viral challenge and stimulation of IFN-β (Zhou et al., 2016; Xie et al., 2018). Despite recent advances in molecular understanding of bat viral tolerance, the consequences of this unique bat immunity on within-host virus dynamics—and its implications for understanding zoonotic emergence—have yet to be elucidated.

The field of ‘virus dynamics’ was first developed to describe the mechanistic underpinnings of long-term patterns of steady-state viral load exhibited by patients in chronic phase infections with HIV, who appeared to produce and clear virus at equivalent rates (Nowak and May, 2000; Ho et al., 1995). Models of simple target cell depletion, in which viral load is dictated by a bottom-up resource supply of infection-susceptible host cells, were first developed for HIV (Perelson, 2002) but have since been applied to other chronic infections, including hepatitis-C virus (Neumann et al., 1998), hepatitis-B virus (Nowak et al., 1996) and cytomegalovirus (Emery et al., 1999). Recent work has adopted similar techniques to model the within-host dynamics of acute infections, such as influenza A and measles, inspiring debate over the extent to which explicit modeling of top-down immune control can improve inference beyond the basic resource limitation assumptions of the target cell model (Baccam et al., 2006; Pawelek et al., 2012; Saenz et al., 2010; Morris et al., 2018).

To investigate the impact of unique bat immune processes on in vitro viral kinetics, we first undertook a series of virus infection experiments on bat cell lines expressing divergent interferon phenotypes, then developed a theoretical model elucidating the dynamics of within-host viral spread. We evaluated our theoretical model analytically independent of the data, then fit the model to data recovered from in vitro experimental trials in order to estimate rates of within-host virus transmission and cellular progression to antiviral status under diverse assumptions of absent, induced, and constitutive immunity. Finally, we confirmed our findings in spatially-explicit stochastic simulations of fitted time series from our mean field model. We hypothesized that top-down immune processes would overrule classical resource-limitation in bat cell lines described as constitutively antiviral in the literature, offering a testable prediction for models fit to empirical data. We further predicted that the most robust antiviral responses would be associated with the most rapid within-host virus propagation rates but also protect cells against virus-induced mortality to support the longest enduring infections in tissue culture.

Results

Virus infection experiments in antiviral bat cell cultures yield reduced cell mortality and elongated epidemics

We first explored the influence of innate immune phenotype on within-host viral propagation in a series of infection experiments in cell culture. We conducted plaque assays on six-well plate monolayers of three immortalized mammalian kidney cell lines: [1] Vero (African green monkey) cells, which are IFN-defective and thus limited in antiviral capacity (Desmyter et al., 1968); [2] RoNi/7.1 (Rousettus aegyptiacus) cells which demonstrate idiosyncratic induced interferon responses upon viral challenge (Kuzmin et al., 2017; Arnold et al., 2018; Biesold et al., 2011; Pavlovich et al., 2018); and [3] PaKiT01 (Pteropus alecto) cells which constitutively express IFN-α (Zhou et al., 2016; Crameri et al., 2009). To intensify cell line-specific differences in constitutive immunity, we carried out infectivity assays with GFP-tagged, replication-competent vesicular stomatitis Indiana viruses: rVSV-G, rVSV-EBOV, and rVSV-MARV, which have been previously described (Miller et al., 2012; Wong et al., 2010). Two of these viruses, rVSV-EBOV and rVSV-MARV, are recombinants for which cell entry is mediated by the glycoprotein of the bat-evolved filoviruses, Ebola (EBOV) and Marburg (MARV), thus allowing us to modulate the extent of structural, as well as immunological, antiviral defense at play in each infection. Previous work in this lab has demonstrated incompatibilities in the NPC1 filovirus receptor which render PaKiT01 cells refractory to infection with rVSV-MARV (Ng and Chandrab, 2018, Unpublished results), making them structurally antiviral, over and above their constitutive expression of IFN-α. All three cell lines were challenged with all three viruses at two multiplicities of infection (MOI): 0.001 and 0.0001. Between 18 and 39 trials were run at each cell-virus-MOI combination, excepting rVSV-MARV infections on PaKiT01 cells at MOI = 0.001, for which only eight trials were run (see Materials and methods; Figure 1—figure supplements 1–3, Supplementary file 1).

Because plaque assays restrict viral transmission neighbor-to-neighbor in two-dimensional cellular space (Howat et al., 2006), we were able to track the spread of GFP-expressing virus-infected cells across tissue monolayers via inverted fluorescence microscopy. For each infection trial, we monitored and re-imaged plates for up to 200 hr of observations or until total monolayer destruction, processed resulting images, and generated a time series of the proportion of infectious-cell occupied plate space across the duration of each trial (see Materials and methods). We used generalized additive models to infer the time course of all cell culture replicates and construct the multi-trial dataset to which we eventually fit our mechanistic transmission model for each cell line-virus-specific combination (Figure 1; Figure 1—figure supplements 1–5).

Figure 1 with 7 supplements see all
Download asset Open asset
Fitted time series of infectious cell proportions from mean field model for rVSV-G, rVSV-EBOV, and rVSV-MARV infections (columns) on Vero, RoNi/7.1, and PaKiT01 cell lines (rows) at MOI = 0.001.

Results are shown for the best fit immune absent model on Vero cells, induced immunity model on RoNi/7.1 cells, and constitutive (for rVSV-VSVG and rVSV-EBOV) and induced (for rVSV-MARV) immunity models on PaKiT01 cells. Raw data across all trials are shown as open circles (statistical smoothers from each trial used for fitting are available in Figure 1—figure supplements 2–3). Model output is shown as a solid crimson line (95% confidence intervals by standard error = red shading). Panel background corresponds to empirical outcome of the average stochastic cell culture trial (persistent infection = white; virus-induced epidemic extinction = gray; immune-mediated epidemic extinction = black). Parameter values are listed in Table 1 and Supplementary file 4. Results for absent/induced/constitutive fitted models across all cell lines are shown in Figure 1—figure supplement 4 (MOI = 0.001) and Figure 1—figure supplement 5 (MOI = 0.0001).

All three recombinant vesicular stomatitis viruses (rVSV-G, rVSV-EBOV, and rVSV-MARV) infected Vero, RoNi/7.1, and PaKiT01 tissue cultures at both focal MOIs. Post-invasion, virus spread rapidly across most cell monolayers, resulting in virus-induced epidemic extinction. Epidemics were less severe in bat cell cultures, especially when infected with the recombinant filoviruses, rVSV-EBOV and rVSV-MARV. Monolayer destruction was avoided in the case of rVSV-EBOV and rVSV-MARV infections on PaKiT01 cells: in the former, persistent viral infection was maintained throughout the 200 hr duration of each experiment, while, in the latter, infection was eliminated early in the time series, preserving a large proportion of live, uninfectious cells across the duration of the experiment. We assumed this pattern to be the result of immune-mediated epidemic extinction (Figure 1). Patterns from MOI = 0.001 were largely recapitulated at MOI = 0.0001, though at somewhat reduced total proportions (Figure 1—figure supplement 5).

A theoretical model fit to in vitro data recapitulates expected immune phenotypes for bat cells

We next developed a within-host model to fit to these data to elucidate the effects of induced and constitutive immunity on the dynamics of viral spread in host tissue (Figure 1). The compartmental within-host system mimicked our two-dimensional cell culture monolayer, with cells occupying five distinct infection states: susceptible (S), antiviral (A), exposed (E), infectious (I), and dead (D). We modeled exposed cells as infected but not yet infectious, capturing the ‘eclipse phase’ of viral integration into a host cell which precedes viral replication. Antiviral cells were immune to viral infection, in accordance with the 'antiviral state' induced from interferon stimulation of ISGs in tissues adjacent to infection (Stetson and Medzhitov, 2006). Because we aimed to translate available data into modeled processes, we did not explicitly model interferon dynamics but instead scaled the rate of cell progression from susceptible to antiviral (ρ) by the proportion of exposed cells (globally) in the system. In systems permitting constitutive immunity, a second rate of cellular acquisition of antiviral status (ε) additionally scaled with the global proportion of susceptible cells in the model. Compared with virus, IFN particles are small and highly diffusive, justifying this global signaling assumption at the limited spatial extent of a six-well plate and maintaining consistency with previous modeling approximations of IFN signaling in plaque assay (Howat et al., 2006).

To best represent our empirical monolayer system, we expressed our state variables as proportions (PS, PA, PE, PI, and PD), under assumptions of frequency-dependent transmission in a well-mixed population (Keeling and Rohani, 2008), though note that the inclusion of PD (representing the proportion of dead space in the modeled tissue) had the functional effect of varying transmission with infectious cell density. This resulted in the following system of ordinary differential equations:

(1) dPSdt=bPD(PS+ PA)−βPSPI−μPS−ρPEPS− εPS+cPA
(2) dPAdt=ρPEPS+ εPS−cPA−μPA
(3) dPEdt=βPSPI-σPE-μPE
(4) dPIdt=σPE-αPI-μPI
(5) dPDdt=μ(PS+PE+ PI+ PA)+αPI−bPD(PS+ PA)

We defined 'induced immunity' as complete, modeling all cells as susceptible to viral invasion at disease-free equilibrium, with defenses induced subsequent to viral exposure through the term ρ. By contrast, we allowed the extent of constitutive immunity to vary across the parameter range of ε > 0, defining a 'constitutive' system as one containing any antiviral cells at disease-free equilibrium. In fitting this model to tissue culture data, we independently estimated both ρ and ε, as well as the cell-to-cell transmission rate, β, for each cell-virus combination. Since the extent to which constitutively-expressed IFN-α is constitutively translated into functional protein is not yet known for bat hosts (Zhou et al., 2016), this approach permitted our tissue culture data to drive modeling inference: even in PaKiT01 cell lines known to constitutively express IFN-α, the true constitutive extent of the system (i.e. the quantity of antiviral cells present at disease-free equilibrium) was allowed to vary through estimation of ε. For the purposes of model-fitting, we fixed the value of c, the return rate of antiviral cells to susceptible status, at 0. The small spatial scale and short time course (max 200 hours) of our experiments likely prohibited any return of antiviral cells to susceptible status in our empirical system; nonetheless, we retained the term c in analytical evaluations of our model because regression from antiviral to susceptible status is possible over long time periods in vitro and at the scale of a complete organism (Radke et al., 1974; Rasmussen and Farley, 1975; Samuel and Knutson, 1982).

Before fitting to empirical time series, we undertook bifurcation analysis of our theoretical model and generated testable hypotheses on the basis of model outcomes. From our within-host model system (Equation 1-5), we derived the following expression for R0, the pathogen basic reproduction number (Supplementary file 2):

(6) R0=βσ(b-μ)(c+μ)bσ+μα+μc+μ+ε

Pathogens can invade a host tissue culture when R0>1. Rapid rates of constitutive antiviral acquisition (ε) will drive R0<1: tissue cultures with highly constitutive antiviral immunity will be therefore resistant to virus invasion from the outset. Since, by definition, induced immunity is stimulated following initial virus invasion, the rate of induced antiviral acquisition (ρ) is not incorporated into the equation for R0; while induced immune processes can control virus after initial invasion, they cannot prevent it from occurring to begin with. In cases of fully induced or absent immunity (ε=0), the R0 equation thus reduces to a form typical of the classic SEIR model:

(7) R0=βσb-μbα+μσ+μ

At equilibrium, the theoretical, mean field model demonstrates one of three infection states: endemic equilibrium, stable limit cycles, or no infection (Figure 2). Respectively, these states approximate the persistent infection, virus-induced epidemic extinction, and immune-mediated epidemic extinction phenotypes previously witnessed in tissue culture experiments (Figure 1). Theoretically, endemic equilibrium is maintained when new infections are generated at the same rate at which infections are lost, while limit cycles represent parameter space under which infectious and susceptible populations are locked in predictable oscillations. Endemic equilibria resulting from cellular regeneration (i.e. births) have been described in vivo for HIV (Coffin, 1995) and in vitro for herpesvirus plaque assays (Howat et al., 2006), but, because they so closely approach zero, true limit cycles likely only occur theoretically, instead yielding stochastic extinctions in empirical time series.

Figure 2
Download asset Open asset
Two parameter bifurcations of the mean field model, showing variation in the transmission rate, β, against variation in the pathogen-induced mortality rate, α, under diverse immune assumptions.

Panel (A) depicts dynamics under variably constitutive immunity, ranging from absent (left: ε=0) to high (right: ε=.0025). In all panel (A) plots, the rate of induced immune antiviral acquisition (ρ) was fixed at 0.01. Panel (B) depicts dynamics under variably induced immunity, ranging from absent (left: ρ=0) to high (right: ρ=1). In all panel (B) plots, the rate of constitutive antiviral acquisition (ε)) was fixed at 0.0001 Branch point curves are represented as solid lines and Hopf curves as dashed lines. White space indicates endemic equilibrium (persistence), gray space indicates limit cycles, and black space indicates no infection (extinction). Other parameter values for equilibrium analysis were fixed at: b = .025, μ = .001, σ = 1/6, c = 0. Special points from bifurcations analyses are listed in Supplementary file 3.

Bifurcation analysis of our mean field model revealed that regions of no infection (pathogen extinction) were bounded at lower threshold (Branch point) values for β, below which the pathogen was unable to invade. We found no upper threshold to invasion for β under any circumstances (i.e. β high enough to drive pathogen-induced extinction), but high β values resulted in Hopf bifurcations, which delineate regions of parameter space characterized by limit cycles. Since limit cycles so closely approach zero, high βs recovered in this range would likely produce virus-induced epidemic extinctions under experimental conditions. Under more robust representations of immunity, with higher values for either or both induced (ρ) and constitutive (ε) rates of antiviral acquisition, Hopf bifurcations occurred at increasingly higher values for β, meaning that persistent infections could establish at higher viral transmission rates (Figure 2). Consistent with our derivation for R0, we found that the Branch point threshold for viral invasion was independent of changes to the induced immune parameter (ρ) but saturated at high values of ε that characterize highly constitutive immunity (Figure 3).

Figure 3
Download asset Open asset
Two parameter bifurcations of the mean field model, showing variation in the transmission rate, β, against variation in: (A) the induced immunity rate of antiviral acquisition (ρ) and (B) the constitutive immunity rate of antiviral acquisition (ε).

Panels show variation in the extent of immunity, from absent (left) to high (right). Branch point curves are represented as solid lines and Hopf curves as dashed lines. White space indicates endemic equilibrium (persistence), gray space indicates limit cycling, and black space indicates no infection (extinction). Other parameter values for equilibrium analysis were fixed at: b = .025, μ = .001, σ = 1/6, α = 1/6, c = 0. Special points from bifurcations analyses are listed in Supplementary file 3.

We next fit our theoretical model by least squares to each cell line-virus combination, under absent, induced, and constitutive assumptions of immunity. In general, best fit models recapitulated expected outcomes based on the immune phenotype of the cell line in question, as described in the general literature (Table 1; Supplementary file 4). The absent immune model offered the most accurate approximation of IFN-deficient Vero cell time series, the induced immune model best recovered the RoNi/7.1 cell trials, and, in most cases, the constitutive immune model most closely recaptured infection dynamics across constitutively IFN-α-expressing PaKiT01 cell lines (Figure 1, Figure 1—figure supplements 4–5, Supplementary file 4). Ironically, the induced immune model offered a slightly better fit than the constitutive to rVSV-MARV infections on the PaKiT01 cell line (the one cell line-virus combination for which we know a constitutively antiviral cell-receptor incompatibility to be at play). Because constitutive immune assumptions can prohibit pathogen invasion (R0<1), model fits to this time series under constitutive assumptions were handicapped by overestimations of ε, which prohibited pathogen invasion. Only by incorporating an exceedingly rapid rate of induced antiviral acquisition could the model guarantee that initial infection would be permitted and then rapidly controlled.

Table 1
Optimized parameters from best fit deterministic model and spatial approximation at MOI = 0.001
Cell lineVirusImmune assumptionAIC reduction 
from next-best model		Antiviral rateε
[lci – uci] *		ρ
[lci – uci] *		β
[lci – uci] *		Mean field
R0		Spatial
β
VerorVSV-GAbsent200 [0–0]0 [0–0]2.44
[1.52–3.36]		8.7324.418
 rVSV-EBOVAbsent200 [0–0]0 [0–0]1.5
[1.06–1.94]		5.4214.996
 rVSV-MARVAbsent200 [0–0]0 [0–0]0.975
[0.558–1.39]		3.459.752
RoNi/7.1rVSV-GInduced27.03 × 10−50 [0–0]0.089
[0–0.432]		2.47
[1.49–3.45]		10.9124.705
 rVSV-EBOVInduced2.012.87 × 10−50 [0–0]0.0363
[0–0.343]		0.685
[0.451–0.919]		3.046.849
 rVSV-MARVInduced21.40 × 10−50 [0–0]0.0177
[0–0.257]		1.23
[0.917–1.55]		5.4812.324
PaKiT01rVSV-GConstitutive29.9.002090.00602
[0–0.019]		8.26 × 10−8
[0–4.75 × 10−7]		3.45
[1.07–5.84]		6.2034.516
 rVSV-EBOVConstitutive27.9.004990.0478
[0–0.0958]		4.46 × 10−8
[0–4.37 × 10−7]		34.5
[28.7–40.2]		18.82344.821
 rVSV-MARVInduced2.006870 [0–0]13.1
[0–37.9]		3.25
[0–41.3]		8.8332.452

  1. Improvement in AIC from next best model for same cell line-virus-MOI combination. All δ-AIC are reported in Supplementary file 4.

    *lci = lower and uci = upper 95% confidence interval. No confidence interval is shown for spatial β which was fixed at 10 times the estimated mean for the mean field model fits when paired with equivalent values of ε and ρ.

  2. All other parameters were fixed at: b = 0.025 (mean field), 0.15 (spatial); α = 1/6; c = 0; μ = 1/121 (Vero), 1/191 (RoNi/7.1), and 1/84 (PaKiT01).

Robust immunity is linked to rapid within-host virus transmission rates in fitted models

In fitting our theoretical model to in vitro data, we estimated the within-host virus transmission rate (β) and the rate(s) of cellular acquisition to antiviral status (ρ or ρ + ε) (Table 1; Supplementary file 4). Under absent immune assumptions, ρ and ε were fixed at 0 while β was estimated; under induced immune assumptions, ε was fixed at 0 while ρ and β were estimated; and under constitutive immune assumptions, all three parameters (ρ, ε, and β) were simultaneously estimated for each cell-virus combination. Best fit parameter estimates for MOI=0.001 data are visualized in conjunction with β – ρ and β – ε bifurcations in Figure 4; all general patterns were recapitulated at lower values for β on MOI=0.0001 trials (Figure 4—figure supplement 1).

As anticipated, the immune absent model (a simple target cell model) offered the best fit to IFN-deficient Vero cell infections (Figure 4; Table 1; Supplementary file 4). Among Vero cell trials, infections with rVSV-G produced the highest β estimates, followed by infections with rVSV-EBOV and rVSV-MARV. Best fit parameter estimates on Vero cell lines localized in the region of parameter space corresponding to theoretical limit cycles, consistent with observed virus-induced epidemic extinctions in stochastic tissue cultures.

Figure 4 with 1 supplement see all
Download asset Open asset
Best fit parameter estimates for β and ρ or ε from mean-field model fits to MOI=0.001 time series data, atop (A,B) β – ρ and (C) β – ε bifurcation.

Fits and bifurcations are grouped by immune phenotype: (A) absent; (B) induced; (C) constitutive immunity, with cell lines differentiated by shape (Vero=circles; RoNi/7.1 = triangles; PaKiT01=squares) and viral infections by color (rVSV-G = green, rVSV-EBOV = magenta, rVSV-MARV = blue). Note that y-axis values are ten-fold higher in panel (C). Branch point curves (solid lines) and Hopf curves (dashed lines) are reproduced from Figure 3. White space indicates endemic equilibrium (pathogen persistence), gray space indicates limit cycling (virus-induced epidemic extinction), and black space indicates no infection (immune-mediated pathogen extinction). In panel (A) and (B), ε is fixed at 0; in panel (C), ρ is fixed at 5x10−8 for bifurcation curves and estimated at 4x10−8 and 8x10−8 for rVSV-EBOV and rVSV-G parameter points, respectively. Other parameter values were fixed at: b = .025, μ = 0.001, σ = 1/6, α = 1/6, and c = 0 across all panels. Raw fitted values and corresponding 95% confidence intervals for β, ρ, and ε, background parameter values, and AIC recovered from model fit, are reported in Supplementary file 4. Parameter fits at MOI=0.0001 are visualized in Figure 4—figure supplement 1.

In contrast to Vero cells, the induced immunity model offered the best fit to all RoNi/7.1 data, consistent with reported patterns in the literature and our own validation by qPCR (Table 1; Figure 1—figure supplement 6; Arnold et al., 2018; Kuzmin et al., 2017; Biesold et al., 2011; Pavlovich et al., 2018). As in Vero cell trials, we estimated highest β values for rVSV-G infections on RoNi/7.1 cell lines but here recovered higher β estimates for rVSV-MARV than for rVSV-EBOV. This reversal was balanced by a higher estimated rate of acquisition to antiviral status (ρ) for rVSV-EBOV versus rVSV-MARV. In general, we observed that more rapid rates of antiviral acquisition (either induced, ρ, constitutive, ε, or both) correlated with higher transmission rates (β). When offset by ρ, β values estimated for RoNi/7.1 infections maintained the same amplitude as those estimated for immune-absent Vero cell lines but caused gentler epidemics and reduced cellular mortality (Figure 1). RoNi/7.1 parameter estimates localized in the region corresponding to endemic equilibrium for the deterministic, theoretical model (Figure 4), yielding less acute epidemics which nonetheless went extinct in stochastic experiments.

Finally, rVSV-G and rVSV-EBOV trials on PaKiT01 cells were best fit by models assuming constitutive immunity, while rVSV-MARV infections on PaKiT01 were matched equivalently by models assuming either induced or constitutive immunity—with induced models favored over constitutive in AIC comparisons because one fewer parameter was estimated (Figure 1—figure supplements 4–5; Supplementary file 4). For all virus infections, PaKiT01 cell lines yielded β estimates a full order of magnitude higher than Vero or RoNi/7.1 cells, with each β balanced by an immune response (either ρ, or ρ combined with ε) also an order of magnitude higher than that recovered for the other cell lines (Figure 4; Table 1). As in RoNi/7.1 cells, PaKiT01 parameter fits localized in the region corresponding to endemic equilibrium for the deterministic theoretical model. Because constitutive immune processes can actually prohibit initial pathogen invasion, constitutive immune fits to rVSV-MARV infections on PaKiT01 cell lines consistently localized at or below the Branch point threshold for virus invasion (R0=1). During model fitting for optimization of ε, any parameter tests of ε values producing R0<1 resulted in no infection and, consequently, produced an exceedingly poor fit to infectious time series data. In all model fits assuming constitutive immunity, across all cell lines, parameter estimates for ρ and ε traded off, with one parameter optimized at values approximating zero, such that the immune response was modeled as almost entirely induced or entirely constitutive (Table 1; Supplementary file 4). For RoNi/7.1 cells, even when constitutive immunity was allowed, the immune response was estimated as almost entirely induced, while for rVSV-G and rVSV-EBOV fits on PaKiT01 cells, the immune response optimized as almost entirely constitutive. For rVSV-MARV on PaKiT01 cells, however, estimation of ρ was high under all assumptions, such that any additional antiviral contributions from ε prohibited virus from invading at all. The induced immune model thus produced a more parsimonious recapitulation of these data because virus invasion was always permitted, then rapidly controlled.

Antiviral cells safeguard live cells against rapid cell mortality to elongate epidemic duration in vitro

In order to compare the relative contributions of each cell line’s disparate immune processes to epidemic dynamics, we next used our mean field parameter estimates to calculate the initial ‘antiviral rate’—the initial accumulation rate of antiviral cells upon virus invasion for each cell-virus-MOI combination—based on the following equation:

(8) AntiviralRate=ρPEPs−ϵPs

where PE was calculated from the initial infectious dose (MOI) of each infection experiment and PS was estimated at disease-free equilibrium:

(9) PE= 1-e-MOI
(10) PS= (b-μ)(c+μ)b(c+μ+ε)

Because ρ and ε both contribute to this initial antiviral rate, induced and constitutive immune assumptions are capable of yielding equally rapid rates, depending on parameter fits. Indeed, under fully induced immune assumptions, the induced antiviral acquisition rate (ρ) estimated for rVSV-MARV infection on PaKiT01 cells was so high that the initial antiviral rate exceeded even that estimated under constitutive assumptions for this cell-virus combination (Supplementary file 4). In reality, we know that NPC1 receptor incompatibilities make PaKiT01 cell lines constitutively refractory to rVSV-MARV infection (Ng and Chandrab, 2018, Unpublished results) and that PaKiT01 cells also constitutively express the antiviral cytokine, IFN-α. Model fitting results suggest that this constitutive expression of IFN-α may act more as a rapidly inducible immune response following virus invasion than as a constitutive secretion of functional IFN-α protein. Nonetheless, as hypothesized, PaKiT01 cell lines were by far the most antiviral of any in our study—with initial antiviral rates estimated several orders of magnitude higher than any others in our study, under either induced or constitutive assumptions (Table 1; Supplementary file 4). RoNi/7.1 cells displayed the second-most-pronounced signature of immunity, followed by Vero cells, for which the initial antiviral rate was essentially zero even under forced assumptions of induced or constitutive immunity (Table 1; Supplementary file 4).

Using fitted parameters for β and ε, we additionally calculated R0, the basic reproduction number for the virus, for each cell line-virus-MOI combination (Table 1; Supplementary file 4). We found that R0 was essentially unchanged across differing immune assumptions for RoNi/7.1 and Vero cells, for which the initial antiviral rate was low. In the case of PaKiT01 cells, a high initial antiviral rate under either induced or constitutive immunity resulted in a correspondingly high estimation of β (and, consequently, R0) which still produced the same epidemic curve that resulted from the much lower estimates for β and R0 paired with absent immunity. These findings suggest that antiviral immune responses protect host tissues against virus-induced cell mortality and may facilitate the establishment of more rapid within-host transmission rates.

Total monolayer destruction occurred in all cell-virus combinations excepting rVSV-EBOV infections on RoNi/7.1 cells and rVSV-EBOV and rVSV-MARV infections on PaKiT01 cells. Monolayer destruction corresponded to susceptible cell depletion and epidemic turnover where R-effective (the product of R0 and the proportion susceptible) was reduced below one (Figure 5). For rVSV-EBOV infections on RoNi/7.1, induced antiviral cells safeguarded remnant live cells, which birthed new susceptible cells late in the time series. In rVSV-EBOV and rVSV-MARV infections on PaKiT01 cells, this antiviral protection halted the epidemic (Figure 5; R-effective <1) before susceptibles fully declined. In the case of rVSV-EBOV on PaKiT01, the birth of new susceptibles from remnant live cells protected by antiviral status maintained late-stage transmission to facilitate long-term epidemic persistence. Importantly, under fixed parameter values for the infection incubation rate (σ) and infection-induced mortality rate (α), models were unable to reproduce the longer-term infectious time series captured in data from rVSV-EBOV infections on PaKiT01 cell lines without incorporation of cell births, an assumption adopted in previous modeling representations of IFN-mediated viral dynamics in tissue culture (Howat et al., 2006). In our experiments, we observed that cellular reproduction took place as plaque assays achieved confluency.

Figure 5 with 3 supplements see all
Download asset Open asset
Fitted time series of susceptible (green shading) and antiviral (blue shading) cell proportions from the mean field model for rVSV-G, rVSV-EBOV, and rVSV-MARV infections (columns) on Vero, RoNi/7.1, and PaKiT01 cell lines (rows) at MOI = 0.001.

Results are shown for the best fit immune absent model on Vero cells, induced immunity model on RoNi/7.1 cells and constitutive (rVSV-G and rVSV-EBOV) and induced (rVSV-MARV) immune models on PaKiT01 cells. Combined live, uninfectious cell populations (S + A + E) are shown in tan shading, with raw live, uninfectious cell data from Hoechst stains visualized as open circles. The right-hand y-axis corresponds to R-effective (pink solid line) across each time series; R-effective = 1 is a pink dashed, horizontal line. Panel background corresponds to empirical outcome of the average stochastic cell culture trial (persistent infection = white; virus-induced epidemic extinction = gray; immune-mediated epidemic extinction = black). Parameter values are listed in Supplementary file 4 and results for absent/induced/constitutive fitted models across all cell lines in Figure 5—figure supplement 1 (MOI = 0.001) and Figure 5—figure supplement 2 (MOI = 0.0001).

Finally, because the protective effect of antiviral cells is more clearly observable spatially, we confirmed our results by simulating fitted time series in a spatially-explicit, stochastic reconstruction of our mean field model. In spatial simulations, rates of antiviral acquisition were fixed at fitted values for ρ and ε derived from mean field estimates, while transmission rates (β) were fixed at values ten times greater than those estimated under mean field conditions, accounting for the intensification of parameter thresholds permitting pathogen invasion in local spatial interactions (see Materials and methods; Videos 1–3; Figure 5—figure supplement 3; Supplementary file 5; Webb et al., 2007). In immune capable time series, spatial antiviral cells acted as ‘refugia’ which protected live cells from infection as each initial epidemic wave ‘washed’ across a cell monolayer. Eventual birth of new susceptibles from these living refugia allowed for sustained epidemic transmission in cases where some infectious cells persisted at later timepoints in simulation (Videos 1–3; Figure 5—figure supplement 3).

Video 1
Download asset
Two hundred hour time series of spatial stochastic model for rVSV-EBOV infection on 10,000 cell grid for PaKiT01, assuming conditions of absent immunity: (A) spatial spread of infection, (B) time series of state variables.

Parameter values are listed in Supplementary file 4.

Video 2
Download asset
Two hundred hour time series of spatial stochastic model for rVSV-EBOV infection on 10,000 cell grid for PaKiT01, assuming conditions of induced immunity: (A) spatial spread of infection, (B) time series of state variables.

Parameter values are listed in Supplementary file 4.

Video 3
Download asset
Two hundred hour time series of spatial stochastic model for rVSV-EBOV infection on 10,000 cell grid for PaKiT01, assuming conditions of constitutive immunity: (A) spatial spread of infection, (B) time series of state variables.

Parameter values are listed in Supplementary file 4.

Discussion

Bats are reservoirs for several important emerging zoonoses but appear not to experience disease from otherwise virulent viral pathogens. Though the molecular biological literature has made great progress in elucidating the mechanisms by which bats tolerate viral infections (Zhou et al., 2016; Ahn et al., 2019; Xie et al., 2018; Pavlovich et al., 2018; Zhang et al., 2013), the impact of unique bat immunity on virus dynamics within-host has not been well-elucidated. We used an innovative combination of in vitro experimentation and within-host modeling to explore the impact of unique bat immunity on virus dynamics. Critically, we found that bat cell lines demonstrated a signature of enhanced interferon-mediated immune response, of either constitutive or induced form, which allowed for establishment of rapid within-host, cell-to-cell virus transmission rates (β). These results were supported by both data-independent bifurcation analysis of our mean field theoretical model, as well as fitting of this model to viral infection time series established in bat cell culture. Additionally, we demonstrated that the antiviral state induced by the interferon pathway protects live cells from mortality in tissue culture, resulting in in vitro epidemics of extended duration that enhance the probability of establishing a long-term persistent infection. Our findings suggest that viruses evolved in bat reservoirs possessing enhanced IFN capabilities could achieve more rapid within-host transmission rates without causing pathology to their hosts. Such rapidly-reproducing viruses would likely generate extreme virulence upon spillover to hosts lacking similar immune capacities to bats.

To achieve these results, we first developed a novel, within-host, theoretical model elucidating the effects of unique bat immunity, then undertook bifurcation analysis of the model’s equilibrium properties under immune absent, induced, and constitutive assumptions. We considered a cell line to be constitutively immune if possessing any number of antiviral cells at disease-free equilibrium but allowed the extent of constitutive immunity to vary across the parameter range for ε, the constitutive rate of antiviral acquisition. In deriving the equation for ε, the basic reproduction number, which defines threshold conditions for virus invasion of a tissue (R0>1), we demonstrated how the invasion threshold is elevated at high values of constitutive antiviral acquisition, ε. Constitutive immune processes can thus prohibit pathogen invasion, while induced responses, by definition, can only control infections post-hoc. Once thresholds for pathogen invasion have been met, assumptions of constitutive immunity will limit the cellular mortality (virulence) incurred at high transmission rates. Regardless of mechanism (induced or constitutive), interferon-stimulated antiviral cells appear to play a key role in maintaining longer term or persistent infections by safeguarding susceptible cells from rapid infection and concomitant cell death.

Fitting of our model to in vitro data supported expected immune phenotypes for different bat cell lines as described in the literature. Simple target cell models that ignore the effects of immunity best recapitulated infectious time series derived from IFN-deficient Vero cells, while models assuming induced immune processes most accurately reproduced trials derived from RoNi/7.1 (Rousettus aegyptiacus) cells, which possess a standard virus-induced IFN-response. In most cases, models assuming constitutive immune processes best recreated virus epidemics produced on PaKiT01 (Pteropus alecto) cells, which are known to constitutively express the antiviral cytokine, IFN-α (Zhou et al., 2016). Model support for induced immune assumptions in fits to rVSV-MARV infections on PaKiT01 cells suggests that the constitutive IFN-α expression characteristic of P. alecto cells may represent more of a constitutive immune priming process than a perpetual, functional, antiviral defense. Results from mean field model fitting were additionally confirmed in spatially explicit stochastic simulations of each time series.

As previously demonstrated in within-host models for HIV (Coffin, 1995; Perelson et al., 1996; Nowak et al., 1995; Bonhoeffer et al., 1997; Ho et al., 1995), assumptions of simple target-cell depletion can often provide satisfactory approximations of viral dynamics, especially those reproduced in simple in vitro systems. Critically, our model fitting emphasizes the need for incorporation of top-down effects of immune control in order to accurately reproduce infectious time series derived from bat cell tissue cultures, especially those resulting from the robustly antiviral PaKiT01 P. alecto cell line. These findings indicate that enhanced IFN-mediated immune pathways in bat reservoirs may promote elevated within-host virus replication rates prior to cross-species emergence. We nonetheless acknowledge the limitations imposed by in vitro experiments in tissue culture, especially involving recombinant viruses and immortalized cell lines. Future work should extend these cell culture studies to include measurements of multiple state variables (i.e. antiviral cells) to enhance epidemiological inference.

The continued recurrence of Ebola epidemics across central Africa highlights the importance of understanding bats’ roles as reservoirs for virulent zoonotic disease. The past decade has born witness to emerging consensus regarding the unique pathways by which bats resist and tolerate highly virulent infections (Brook and Dobson, 2015; Xie et al., 2018; Zhang et al., 2013; Ahn et al., 2019; Zhou et al., 2016; Ng et al., 2015; Pavlovich et al., 2018). Nonetheless, an understanding of the mechanisms by which bats support endemic pathogens at the population level, or promote the evolution of virulent pathogens at the individual level, remains elusive. Endemic maintenance of infection is a defining characteristic of a pathogen reservoir (Haydon et al., 2002), and bats appear to merit such a title, supporting long-term persistence of highly transmissible viral infections in isolated island populations well below expected critical community sizes (Peel et al., 2012). Researchers debate the relative influence of population-level and within-host mechanisms which might explain these trends (Plowright et al., 2016), but increasingly, field data are difficult to reconcile without acknowledgement of a role for persistent infections (Peel et al., 2018; Brook et al., 2019). We present general methods to study cross-scale viral dynamics, which suggest that within-host persistence is supported by robust antiviral responses characteristic of bat immune processes. Viruses which evolve rapid replication rates under these robust antiviral defenses may pose the greatest hazard for cross-species pathogen emergence into spillover hosts with immune systems that differ from those unique to bats.

Materials and methods

Key resources table
Reagent type (species)
or resource		DesignationSource or referenceIdentifiersAdditional information
Cell line (Vero)Kidney (normal, epithelial, adult)ATCCCCL-81
Cell line (Rousettus aegyptiacus)Kidney (normal, epithelial, adult)(Biesold et al., 2011; Kühl et al., 2011)RoNi/7.1
Cell line (Pteropus alecto)Kidney (normal, epithelial, adult)(Crameri et al., 2009)PaKiT01
Virus strainReplication competent, recombinant vesicular stomatitis Indiana virus expressing eGFP(Miller et al., 2012; Wong et al., 2010)rVSV-G
Virus strainReplication competent, recombinant vesicular stomatitis Indiana virus expressing eGFP and EBOV GP in place of VSV G(Miller et al., 2012; Wong et al., 2010)rVSV-EBOV
Virus strainReplication competent, recombinant vesicular stomatitis Indiana virus expressing eGFP and MARV GP in place of VSV G(Miller et al., 2012; Wong et al., 2010)rVSV-MARV
ReagentHoechst 33342 Fluorescent StainThermoFishercat #: 62249
ReagentL-Glutamine SolutionThermoFishercat #: 25030081
ReagentGibco HEPESThermoFishercat #: 15630080
ReagentiTaq Universal SYBR Green SupermixBioRadcat #: 1725120
Commercial assay or kitQuick RNA Mini Prep KitZymocat #:
R1054
Commercial assay or kitInvitrogen Superscript III cDNA Synthesis KitThermoFishercat #: 18080051
SoftwareMatCont (version 2.2)(Dhooge et al., 2008)MatCont
RR version 3.6.0(R Development Core Team, 2019)R

  1. *Note that primers for R. aegyptiacus and P. alecto β-Actin, IFN-α, and IFN-β genes are listed in the Supplementary file 6.

Cell culture experiments

Cells

All experiments were carried out on three immortalized mammalian kidney cell lines: Vero (African green monkey), RoNi/7.1 (Rousettus aegyptiacus) (Kühl et al., 2011; Biesold et al., 2011) and PaKiT01 (Pteropus alecto) (Crameri et al., 2009). The species identifications of all bat cell lines was confirmed morphologically and genetically in the publications in which they were originally described (Kühl et al., 2011; Biesold et al., 2011; Crameri et al., 2009). Vero cells were obtained from ATCC.

Monolayers of each cell line were grown to 90% confluency (~9×105 cells) in 6-well plates. Cells were maintained in a humidified 37 °C, 5% CO2 incubator and cultured in Dulbecco’s modified Eagle medium (DMEM) (Life Technologies, Grand Island, NY), supplemented with 2% fetal bovine serum (FBS) (Gemini Bio Products, West Sacramento, CA), and 1% penicillin-streptomycin (Life Technologies). Cells were tested monthly for mycoplasma contamination while experiments were taking place; all cells assayed negative for contamination at every testing.

Previous work has demonstrated that all cell lines used are capable of mounting a type I IFN response upon viral challenge, with the exception of Vero cells, which possess an IFN-β deficiency (Desmyter et al., 1968; Rhim et al., 1969; Emeny and Morgan, 1979). RoNi/7.1 cells have been shown to mount idiosyncratic induced IFN defenses upon viral infection (Pavlovich et al., 2018; Kuzmin et al., 2017; Arnold et al., 2018; Kühl et al., 2011; Biesold et al., 2011), while PaKiT01 cells are known to constitutively express the antiviral cytokine, IFN-α (Zhou et al., 2016). This work is the first documentation of IFN signaling induced upon challenge with the particular recombinant VSVs outlined below. We verified known antiviral immune phenotypes via qPCR. Results were consistent with the literature, indicating a less pronounced role for interferon defense against viral infection in RoNi/7.1 versus PaKiT01 cells.

Viruses

Replication-capable recombinant vesicular stomatitis Indiana viruses, expressing filovirus glycoproteins in place of wild type G (rVSV-G, rVSV-EBOV, and rVSV-MARV) have been previously described (Wong et al., 2010; Miller et al., 2012). Viruses were selected to represent a broad range of anticipated antiviral responses from host cells, based on a range of past evolutionary histories between the virus glycoprotein mediating cell entry and the host cell’s entry receptor. These interactions ranged from the total absence of evolutionary history in the case of rVSV-G infections on all cell lines to a known receptor-level cell entry incompatibility in the case of rVSV-MARV infections on PaKiT01 cell lines.

To measure infectivities of rVSVs on each of the cell lines outlined above, so as to calculate the correct viral dose for each MOI, NH4Cl (20 mM) was added to infected cell cultures at 1–2 hr post-infection to block viral spread, and individual eGFP-positive cells were manually counted at 12–14 hr post-infection.

Innate immune phenotypes via qPCR of IFN genes

Request a detailed protocol

Previously published work indicates that immortalized kidney cell lines of Rousettus aegyptiacus (RoNi/7.1) and Pteropus alecto (PaKiT01) exhibit different innate antiviral immune phenotypes through, respectively, induced (Biesold et al., 2011; Pavlovich et al., 2018; Kühl et al., 2011; Arnold et al., 2018) and constitutive (Zhou et al., 2016) expression of type I interferon genes. We verified these published phenotypes on our own cell lines infected with rVSV-G, rVSV-EBOV, and rVSV-MARV via qPCR of IFN-α and IFN-β genes across a longitudinal time series of infection.

Specifically, we carried out multiple time series of infection of each cell line with each of the viruses described above, under mock infection conditions and at MOIs of 0.0001 and 0.001—with the exception of rVSV-MARV on PaKiT01 cell lines, for which infection was only performed at MOI = 0.0001 due to limited viral stocks and the extremely low infectivity of this virus on this cell line (thus requiring high viral loads for initial infection). All experiments were run in duplicate on 6-well plates, such that a typical plate for any of the three viruses had two control (mock) wells, two MOI = 0.0001 wells and two MOI = 0.001 wells, excepting PaKiT01 plates, which had two control and four MOI = 0.0001 wells at a given time. We justify this PaKiT01 exemption through the expectation that IFN-α expression is constitutive for these cells, and by the assumption that any expression exhibited at the lower MOI should also be present at the higher MOI.

For these gene expression time series, four 6-well plates for each cell line–virus combination were incubated with virus for one hour at 37 °C. Following incubation, virus was aspirated off, and cell monolayers were washed in PBS, then covered with an agar plaque assay overlay to mimic conditions under which infection trials were run. Plates were then harvested sequentially at timepoints of roughly 5, 10, 15, and 20 hr post-infection (exact timing varied as multiple trials were running simultaneously). Upon harvest of each plate, agar overlay was removed, and virus was lysed and RNA extracted from cells using the Zymo Quick RNA Mini Prep kit, according to the manufacturer’s instructions and including the step for cellular DNA digestion. Post-extraction, RNA quality was verified via nanodrop, and RNA was converted to cDNA using the Invitrogen Superscript III cDNA synthesis kit, according to the manufacturer’s instructions. cDNA was then stored at 4 °C and as a frozen stock at −20 °C to await qPCR.

We undertook qPCR of cDNA to assess expression of the type I interferon genes, IFN-α and IFN-β, and the housekeeping gene, β-Actin, using primers previously reported in the literature (Supplementary file 6). For qPCR, 2 µl of each cDNA sample was incubated with 7 µl of deionized water, 1 µl of 5 UM forward/reverse primer mix and 10 µl of iTaq Universal SYBR Green, then cycled on a QuantStudio3 Real-Time PCR machine under the following conditions: initial denaturation at 94 °C for 2 min followed by 40 cycles of: denaturation at 95 °C (5 s), annealing at 58 °C (15 s), and extension at 72 °C (10 s).

We report simple δ-Ct values for each run, with raw Ct of the target gene of interest (IFN-α or IFN-β) subtracted from raw Ct of the β-Actin housekeeping gene in Figure 1—figure supplement 6. Calculation of fold change upon viral infection in comparison to mock using the δ-δ-Ct method (Livak and Schmittgen, 2001) was inappropriate in this case, as we wished to demonstrate constitutive expression of IFN-α in PaKiT01 cells, whereby data from mock cells was identical to that produced from infected cells.

Plaque assays and time series imaging

Request a detailed protocol

After being grown to ~90% confluency, cells were incubated with pelleted rVSVs expressing eGFP (rVSV-G, rVSV-EBOV, rVSV-MARV). Cell lines were challenged with both a low (0.0001) and high (0.001) multiplicity of infection (MOI) for each virus. In a cell monolayer infected at a given MOI (m), the proportion of cells (P), infected by k viral particles can be described by the Poisson distribution: R0, such that the number of initially infected cells in an experiment equals: R0>1. We assumed that a ~90% confluent culture at each trial’s origin was comprised of ~9x105 cells and conducted all experiments at MOIs of 0.0001 and 0.001, meaning that we began each trial by introducing virus to, respectively, ~81 or 810 cells, representing the state variable ‘E’ in our theoretical model. Low MOIs were selected to best approximate the dynamics of mean field infection and limit artifacts of spatial structuring, such as premature epidemic extinction when growing plaques collide with plate walls in cell culture.

Six-well plates were prepared with each infection in duplicate or triplicate, such that a control well (no virus) and 2–3 wells each at MOI 0.001 and 0.0001 were incubated simultaneously on the same plate. In total, we ran between 18 and 39 trials at each cell-virus-MOI combination, excepting r-VSV-MARV infections on PaKiT01 cells at MOI = 0.001, for which we ran only eight trials due to the low infectivity of this virus on this cell line, which required high viral loads for initial infection. Cells were incubated with virus for one hour at 37 °C. Following incubation, virus was aspirated off, and cell monolayers were washed in PBS, then covered with a molten viscous overlay (50% 2X MEM/L-glutamine; 5% FBS; 3% HEPES; 42% agarose), cooled for 20 min, and re-incubated in their original humidified 37 °C, 5% CO2 environment.

After application of the overlay, plates were monitored periodically using an inverted fluorescence microscope until the first signs of GFP expression were witnessed (~6–9.5 hr post-infection, depending on the cell line and virus under investigation). From that time forward, a square subset of the center of each well (comprised of either 64- or 36-subframes and corresponding to roughly 60% and 40% of the entire well space) was imaged periodically, using a CellInsight CX5 High Content Screening (HCS) Platform with a 4X air objective (ThermoFisher, Inc, Waltham, MA). Microscope settings were held standard across all trials, with exposure time fixed at 0.0006 s for each image. One color channel was imaged, such that images produced show GFP-expressing cells in white and non-GFP-expressing cells in black (Figure 1—figure supplement 1).

Wells were photographed in rotation, as frequently as possible, from the onset of GFP expression until the time that the majority of cells in the well were surmised to be dead, GFP expression could no longer be detected, or early termination was desired to permit Hoechst staining.

In the case of PaKiT01 cells infected with rVSV-EBOV, where an apparently persistent infection established, the assay was terminated after 200+ hours (8+ days) of continuous observation. Upon termination of all trials, cells were fixed in formaldehyde (4% for 15 min), incubated with Hoechst stain (0.0005% for 15 min) (ThermoFisher, Inc, Waltham, MA), then imaged at 4X on the CellInsight CX5 High Content Screening (HCS) Platform. The machine was allowed to find optimal focus for each Hoechst stain image. One color channel was permitted such that images produced showed live nuclei in white and dead cells in black.

Hoechst staining

Request a detailed protocol

Hoechst stain colors cellular DNA, and viral infection is thought to interfere with the clarity of the stain (Dembowski and DeLuca, 2015). As such, infection termination, cell fixation, and Hoechst staining enables generation of a rough time series of uninfectious live cells (i.e. susceptible + antiviral cells) to complement the images which produced time series of proportions infectious. Due to uncertainty over the exact epidemic state of Hoechst-stained cells (i.e. exposed but not yet infectious cells may still stain), we elected to fit our models only to the infectious time series derived from GFP-expressing images and used Hoechst stain images as a post hoc visual check on our fit only (Figure 5; Figure 5—figure supplements 1–2).

Image processing

Request a detailed protocol

Images recovered from the time series above were processed into binary (‘infectious’ vs. ‘non-infectious’ or, for Hoechst-stained images, ‘live’ vs. ‘dead’) form using the EBImage package (Pau et al., 2010) in R version 3.6 for MacIntosh, after methods further detailed in Supplementary file 7. Binary images were then further processed into time series of infectious or, for Hoechst-stained images, live cells using a series of cell counting scripts. Because of logistical constraints (i.e. many plates of simultaneously running infection trials and only one available imaging microscope), the time course of imaging across the duration of each trial was quite variable. As such, we fitted a series of statistical models to our processed image data to reconstruct reliable values of the infectious proportion of each well per hour for each distinct trial in all cell line–virus-MOI combinations (Figure 1—figure supplements 2–3), as well as for declining live cell counts from control well data derived from the Hoestch time series (Supplementary file 1; Supplementary file 7; Figure 1—figure supplement 7). All original and processed images, image processing and counting code, and resulting time series data are freely available for download at the following FigShare repository: DOI: 10.6084/m9.figshare.8312807.

Mean field model

Theoretical model details

Request a detailed protocol

To derive the expression for R0, the basic pathogen reproductive number in vitro, we used Next Generation Matrix (NGM) techniques (Diekmann et al., 1990; Heffernan et al., 2005), employing Wolfram Mathematica (version 11.2) as an analytical tool. R0 describes the number of new infections generated by an existing infection in a completely susceptible host population; a pathogen will invade a population when R0>1 (Supplementary file 2). We then analyzed stability properties of the system, exploring dynamics across a range of parameter spaces, using MatCont (version 2.2) (Dhooge et al., 2008) for Matlab (version R2018a) (Supplementary file 3).

Theoretical model fitting

Request a detailed protocol

The birth rate, b, and natural mortality rate, μ, balance to yield a population-level growth rate, such that it is impossible to estimate both b and μ simultaneously from total population size data alone. As such, we fixed b at. 025 and estimated μ by fitting an infection-absent version of our mean field model to the susceptible time series derived via Hoechst staining of control wells for each of the three cell lines (Figure 1—figure supplement 7). This yielded a natural mortality rate, μ, corresponding to a lifespan of approximately 121, 191, and 84 hours, respectively, for Vero, RoNi/7.1, and PaKiT01 cell lines (Figure 1—figure supplement 7). We then fixed the virus incubation rate, σ, as the inverse of the shortest observed duration of time from initial infection to the observation of the first infectious cells via fluorescent microscope for all nine cell line – virus combinations (ranging 6 to 9.5 hours). We fixed α, the infection-induced mortality rate, at 1/6, an accepted standard for general viral kinetics (Howat et al., 2006), and held c, the rate of antiviral cell regression to susceptible status, at 0 for the timespan (<200 hours) of the experimental cell line infection trials.

We estimated cell line–virus-MOI-specific values for β, ρ, and ε by fitting the deterministic output of infectious proportions in our mean field model to the full suite of statistical outputs of all trials for each infected cell culture time series (Figure 1—figure supplements 2–3). Fitting was performed by minimizing the sum of squared differences between the deterministic model output and cell line–virus-MOI-specific infectious proportion of the data at each timestep. We optimized parameters for MOI = 0.001 and 0.0001 simultaneously to leverage statistical power across the two datasets, estimating a different transmission rate, β, for trials run at each infectious dose but, where applicable, estimating the same rates of ρ and ε across the two time series. We used the differential equation solver lsoda() in the R package deSolve (Soetaert et al., 2010) to obtain numerical solutions for the mean field model and carried out minimization using the ‘Nelder-Mead’ algorithm of the optim() function in base R. All model fits were conducted using consistent starting guesses for the parameters, β (β = 3), and where applicable, ρ (ρ = 0.001) and ε (ε = 0.001). In the case of failed fits or indefinite hessians, we generated a series of random guesses around the starting conditions and continued estimation until successful fits were achieved.

All eighteen cell line–virus-MOI combinations of data were fit by an immune absent (ε = ρ = 0) version of the theoretical model and, subsequently, an induced immunity (ε = 0; ρ >0) and constitutive immunity (ε >0; ρ >0) version of the model. Finally, we compared fits across each cell line–virus-MOI combination via AIC. In calculating AIC, the number of fitted parameters in each model (k) varied across the immune phenotypes, with one parameter (β) estimated for absent immune assumptions, two (β and ρ) for induced immune assumptions, and three (β, ρ, and ε) for constitutive immune assumptions. The sample size (n) corresponded to the number of discrete time steps across all empirical infectious trials to which the model was fitted for each cell-line virus combination. All fitting and model comparison scripts are freely available for download at the following FigShare repository: DOI: 10.6084/m9.figshare.8312807.

Spatial model simulations

Request a detailed protocol

Finally, we verified all mean field fits in a spatial context, in order to more thoroughly elucidate the role of antiviral cells in each time series. We constructed our spatial model in C++ implemented in R using the packages Rcpp and RcppArmadillo (Eddelbuettel and Francois, 2011; Eddelbuettel and Sanderson, 2017). Following Nagai and Honda (2001) and Howat et al. (2006), we modeled this system on a two-dimensional hexagonal lattice, using a ten-minute epidemic timestep for cell state transitions. At the initialization of each simulation, we randomly assigned a duration of natural lifespan, incubation period, infectivity period, and time from antiviral to susceptible status to all cells in a theoretical monolayer. Parameter durations were drawn from a normal distribution centered at the inverse of the respective fixed rates of μ, σ, α, and c, as reported with our mean field model. Transitions involving the induced (ρ) and constitutive (ε) rates of antiviral acquisition were governed probabilistically and adjusted dynamically at each timestep based on the global environment. As such, we fixed these parameters at the same values estimated in the mean field model, and multiplied both ρ and ε by the global proportion of, respectively, exposed and susceptible cells at a given timestep.

In contrast to antiviral acquisition rates, transitions involving the birth rate (b) and the transmission rate (β) occurred probabilistically based on each cell’s local environment. The birth rate, b, was multiplied by the proportion of susceptible cells within a six-neighbor circumference of a focal dead cell, while β was multiplied by the proportion of infectious cells within a thirty-six neighbor vicinity of a focal susceptible cell, thus allowing viral transmission to extend beyond the immediate nearest-neighbor boundaries of an infectious cell. To compensate for higher thresholds to cellular persistence and virus invasion which occur under local spatial conditions (Webb et al., 2007), we increased the birth rate, b, and the cell-to-cell transmission rate, β, respectively, to six and ten times the values used in the mean field model (Supplementary file 4). We derived these increases based on the assumption that births took place exclusively based on pairwise nearest-neighbor interactions (the six immediately adjacent cells to a focal dead cell), while viral transmission was locally concentrated but included a small (7.5%) global contribution, representing the thirty-six cell surrounding vicinity of a focal susceptible. We justify these increases and derive their origins further in Supplementary file 5.

We simulated ten stochastic spatial time series for all cell-virus combinations under all three immune assumptions at a population size of 10,000 cells and compared model output with data in Figure 5—figure supplement 3. Spatial model code is available for public access at the following FigShare repository: DOI: 10.6084/m9.figshare.8312807.

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files. All images and code used in this study have been made available for download at the following Figshare repository: https://doi.org/10.6084/m9.figshare.8312807.

The following data sets were generated
    1. Brook CE
    2. Ng M
    3. Boots M
    4. Dobson A
    5. Graham A
    6. Grenfell B
    7. Chandran KC
    8. van Leeuwen A
    (2019) figshare
    Data and Code from: Accelerated viral dynamics in bat cell lines, with implications for zoonotic emergence.
    https://doi.org/10.6084/m9.figshare.8312807

References

  25. Book
    1. Keeling MJ
    2. Rohani P
    (2008)
    Modeling Infectious Diseases in Humans and Animals
    Princeton University Press: Princeton, NJ.
  38. Book
    1. Nowak MA
    2. May RM
    (2000)
    Virus Dynamics: Mathematical Principles of Immunology and Virology
    Oxford: Oxford University Press.
  47. Software
    1. R Development Core Team
    (2019) R: A Language and Environment for Statistical Computing
    R Foundation for Statistical Computing, Vienna, Austria.
    http://www.r-project.org
    1. Samuel CE
    2. Knutson GS
    (1982)
    Mechanism of interferon action
    Journal of Biological Chemistry 257:11791–11795.

Decision letter

  1. Dan Haydon
    Reviewing Editor; University of Glasgow, United Kingdom
  2. Neil M Ferguson
    Senior Editor; Imperial College London, United Kingdom
  3. Dan Haydon
    Reviewer; University of Glasgow, United Kingdom

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This paper is an excellent example of how in vitro models of cell-virus interactions can be used to shape and formulate more general and larger-scale hypotheses about epidemiological dynamics. In this case, the choice of bat cell lines expressing induced and constitutive immune phenotypes enables estimates of different viral propagation rates. The results suggest that if bat cells do have greater constitutive immunity, this could lead to situations in which viruses that do propagate in bats will do so with much greater vigour (and possibly virulence and transmissability) should they 'spill-over' into non-bat hosts. The paper should be of wide general interest to those with interests in emerging disease dynamics and to quantitative biologists interested in the mathematical modelling of in vitro systems.

Decision letter after peer review:

Thank you for submitting your article "Within-host dynamics of virulent viruses in bat reservoirs for emerging zoonotic disease" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Dan Haydon as the Reviewing Editor and Reviewer #3, and the evaluation has been overseen by Neil Ferguson as the Senior Editor.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

There is a view that bats (both as individuals and populations) are able to maintain a large number of virus species (often with little or no pathology) that are highly pathogenic outside of this order, and this paper seeks mechanistic reasons for why this might be. Specifically the authors set out to explore whether this might be due to uniquely constitutive immune capabilities of bat cells. To explore this hypothesis, the authors study the dynamics of 3 strains of vesicular stomatitis virus that express different sorts of glycoprotein, in 3 different cell cultures: Vero (non-bat cells with limited antiviral capabilities); Egyptian Fruit Bat cells which demonstrate an idiosyncratic induced interferon response, and Black Fruit Bat cells which constitutively express interferon. Well-mixed and imperfectly-mixed frequency dependent models of the infection dynamics are fitted to the means of these replicated time-series, and the models constructed to reflect the possibly differential anti-viral capabilities of the cell lines. The parameter estimates from the model fitting process lead to the key conclusion that 'induced immune responses favor slower cell-to-cell transmission rates (a proxy for viral replication rates) […] while constitutive immunity amplifies cellular transmission rates, in conjunction with less rapid antiviral responses.' The authors further conclude that '[i]f hosts with constitutive immune defenses favor evolution of rapidly transmitting viruses, such pathogens are likely to cause extreme virulence in spillover hosts lacking similarly constitutive defenses'.

One reviewer indicates that if the paper has set out to answer the question as to whether the immune phenotype of bats alters the capacity for viruses to persist within them, the answers are not sufficiently clear to merit publication in eLife. This is an understandable point of view: it seems that none of the best-fitting well-mixed models are constitutive, and only one of the spatially explicit models is constitutive where we might expect it to be (and one where we wouldn't expect it to fit best). However, I have decided to go with the majority view and provide an opportunity for revision, but the authors should consider carefully whether the essential points below can be met (or convincingly rebutted).

Essential revisions:

1) Overall, the reviewers felt that this paper was quite a bit more complicated in the presentation of its results than it needs to be, and I encourage the authors to find ways of simplifying the Results section and linearizing the framing of the key messages. I wondered what the benefits of presenting both the well-mixed and spatial models were? Which is the most appropriate? Do they say importantly different things? Could the manuscript be simplified by focusing on just one?

2) Another key question the authors should consider is whether the data in Figure 4 really support (or even need) the cell turnover embedded within the formulation of the model? Could the main point of the paper be more clearly made by fitting simpler models more consistent with the short time frame of the experimental data, and concentrating on the initial infection spread and sometime declines, rather than a scenario requiring cell turnover? If the cell cultures were capable of maintaining infection over the longer term, why weren't the experiments run over longer time frames to demonstrate this?

Introduction:

3) The authors cannot claim to be studying the within-host dynamics of bat viruses. Rather, they experimentally examine and model the dynamics of recombinant viruses in vitro. The authors should explicitly recognize throughout their manuscript that their study is at this more limited scale and, as a result, that its impact is also limited.

Results/Materials and methods:

4) The three cell lines are presented as examples of no immunity, induced immunity and constitutive immunity. The questions set out at the end of the Introduction are not answered directly (for example at the beginning of the Discussion) with attention instead diverted to measures of rate of spread. Ensure a more identifiable match between the questions posed at the end of the Introduction, and the start of the Discussion.

5) Add more explanation of what these virus lines are and why they are chosen. The reader is not told that rVSV-MARV and rVSV-EBOV are recombinant viruses expressing Ebola and Marburg Virus (Results first paragraph). Moreover, it is not stated why the authors employ two viruses (they do not seem to be interested in contrasting them but they cannot be said to be true replicates) or these two in particular.

6) How many replicates were there for each cell-line/clone combination? The only detail mentioned in this regard is that there were 2-3 technical replicates per plate. Such details should be stated prominently in the Results, since large numbers of points are plotted, but the reader has no way to discern which of them are independent.

7) There is a strong case for requiring the model is fitted to the individual replicates and not the mean of the replicates. If this case is rejected, then careful justification must be provided.

8) The cell cultures are grown to 90% confluency, and are modeled as systems in which there is cell birth and natural death (birth rate about 4x the natural death rate). Readers are likely to be left with a number of questions that the authors should address more explicitly: Where does the birth rate parameter come from? Is this a realistic estimate reflective of the dynamics of the cell culture, or is it a deduction based on the requirement to sustain the infection dynamics (either endemic or oscillatory?). The Materials and methods suggest that it is what would be required to maintain sustainable live cell populations, but is this what actually happens in the culture during the infection? What sort of turnover does this result in, in the absence of infection? What is the justification for frequency over density dependent dynamics? At the heart of these questions is whether it is appropriate to think of these cell cultures as a sort of viral chemostat which could maintain steady state or oscillatory dynamics indefinitely, or whether the system would be better modeled as an infection process in a fixed non-reproducing host population. The latter view would substantially simplify the paper if such a perspective was justifiable.

9) The potential influence of the assumed value for the rate at which antiviral cells regress to susceptibility is unclear. At present, it is not clear where the parameter value used comes from and what impact it could have on the results. This is particularly pertinent since sustained epidemics are ascribed to the antiviral cells returning to susceptibility – can one confidently separate this from birth of new susceptibles?

10) Why is ε not estimated? It seems to be fixed at 0 or 1. In subsection “Fitting of theoretical model to cell culture data demonstrates higher within-host transmission rates under constitutive immune assumptions” this seems to be recognized but Supplementary file 4 doesn't provide any intermediate values. It wasn't clear to why ε is bounded at 1. It is a rate?

11) It is surely odd that as b (and u) is lowered towards zero, R0 also tends to zero. Surely one would expect a virus to be able to spread in a non-renewing monolayer?

Discussion:

12) The use of the word 'favor' in the first paragraph of the Discussion is a little misleading and should be re-worded. These findings are not based on any evolutionary process.

13) The Discussion should make clear that the question of whether bats in general have constitutively active intracellular immunity is far from resolved. As touched on in the Introduction, this partly depends on definition – does constitutive expression of e.g. interferon-α actually provide constitutive protection? What about other immune genes which seem to be defective (e.g. Xie et al., 2018)? More generally, the immune systems of only a very small fraction of bat diversity (and of non-lab animals in general) has been studied, so it remains unclear how well we can generalise to all bats, or whether bats are truly unusual relative to other wild animals. Where possible the authors should list examples of studies showing constitutive immunity specifically in the species (or close relatives of these) thought to be reservoirs for the highly virulent viruses discussed.

https://doi.org/10.7554/eLife.48401.sa1

Author response

Essential revisions:

1) Overall, the reviewers felt that this paper was quite a bit more complicated in the presentation of its results than it needs to be, and I encourage the authors to find ways of simplifying the Results section and linearizing the framing of the key messages. I wondered what the benefits of presenting both the well-mixed and spatial models were? Which is the most appropriate? Do they say importantly different things? Could the manuscript be simplified by focusing on just one?

Thank you for taking the time to decipher our paper, and we apologize for any unnecessary complexity and confusing organization and/or writing. As highlighted above, our new version of the manuscript presents a completely reworked version of the model and, correspondingly, a complete rewrite of the text. The main messages of our paper are threefold: (1) we find support for immune phenotypes suggested from the literature for each of the three cell lines, in particular showing a necessary effect of top-down immune response in describing data for the bat cell lines (especially the constitutively antiviral PaKiT01 cell line); (2) we find that more pronounced immune responses (either induced or constitutive) are correlated with higher cell-to-cell viral transmission rates (β), suggesting that bats’ robust antiviral immune phenotype could be driving the evolution of higher within-host viral replication rates likely to cause virulence upon emergence into hosts lacking bat immune capacities; and (3) IFN-induced antiviral cells promote the maintenance of persistent within-host infections by protecting live cell refugia against rapid viral infection and mortality in the initial epidemic takeoff such that live cells are maintained that can sustain transmission at later timepoints.

Key differences in the new version of the model which depart from the original model include (a) no regression of antiviral cells back to susceptible status and (b) estimation of the rate of constitutive cellular acquisition of antiviral status (represented as ε), rather than fixation as a 0-1 scalar as in the first version of submission. We elaborate on the effects of these major changes to the model in answer to questions #2, 9, and 10 below.

Finally, in response to your questions above, we have largely relegated the spatial model to the supplementary materials and are here using it only to validate and visualize results from our mean field model fitting. The mean field model demonstrates all of our key takeaway messages from this manuscript, with the added benefit of being generalizable to the whole organism when considered theoretically in our bifurcation analysis. However, since our data (in plaque assay) are spatial, we felt that it was important to simulate the dynamics in a spatial context as well. We find that, in doing this, our key takeaways are reaffirmed and validated. We retain the spatial stochastic videos in our supplementary text because they clearly illustrate how the epidemic wave ‘washes over’ antiviral refugia, which can become infected at later timepoints to sustain transmission. This pattern is seen in the mean field model as well (Figure 5, main text, and Figure 5—figure supplements 1 and 2) but much easier to interpret when visualized spatially.

2) Another key question the authors should consider is whether the data in Figure 4 really support (or even need) the cell turnover embedded within the formulation of the model? Could the main point of the paper be more clearly made by fitting simpler models more consistent with the short time frame of the experimental data, and concentrating on the initial infection spread and sometime declines, rather than a scenario requiring cell turnover? If the cell cultures were capable of maintaining infection over the longer term, why weren't the experiments run over longer time frames to demonstrate this?

In the newest version of the manuscript, the original Figure 4 is now Figure 1, showing the fitted infectious time series to these data. In response to this comment, we attempted to fit several versions of the model that incorporated only the natural mortality rate (fitted for each cell line) and the epidemic dynamics of each time series, with no birth rate included in the model. However, these versions of the model were unable to re-capture the extended time series witnessed in the rVSV-EBOV infections on PaKiT01 cell lines, using fitted mortality (μ) and incubation (σ) rates and a fixed, standard infection-induced mortality rate (α= 1/6). Thus, cell turnover was needed in order to recapture the data.

Additionally, we feel confident that cell turnover did take place in the experiments because we grew cells to 90% confluency prior to initiating the plaque assay, then watched them subsequently fill the plate when they could not be split at later timepoints. In including cell turnover, our analysis additionally maintains consistency with previous in vitromodeling of IFN dynamics in cell culture in the general literature (Howat et al., 2006).

Introduction:

3) The authors cannot claim to be studying the within-host dynamics of bat viruses. Rather, they experimentally examine and model the dynamics of recombinant viruses in vitro. The authors should explicitly recognize throughout their manuscript that their study is at this more limited scale and, as a result, that its impact is also limited.

Fair point, and thank you for bringing up this issue. We do believe that our mean field model can be generalized to the level of the whole host organism, as a description of viral spread throughout a host’s tissue, which is why we have undertaken the broad bifurcation analysis, in addition to fitting our model to the in vitrodata. Figure 2 and 3, as well as the derivation for R0 (as a threshold for tissue invasion) presented in our text are, thus, truly representative of within-host dynamics. Nonetheless, we appreciate your concerns and have scaled back our claims to highlight the terms ‘in vitro’ or ‘bat cell culture’ where appropriate, including in the title.

Results/Materials and methods:

4) The three cell lines are presented as examples of no immunity, induced immunity and constitutive immunity. The questions set out at the end of the Introduction are not answered directly (for example at the beginning of the Discussion) with attention instead diverted to measures of rate of spread. Ensure a more identifiable match between the questions posed at the end of the Introduction, and the start of the Discussion.

Thanks for this! Good point, and apologies for any digressions in our writing. As mentioned above, our newest version of the manuscript is a complete rewrite of the original, and we have addressed this discrepancy. The end of our Introduction now highlights our main hypotheses:

“(a) We hypothesized that model fitting would demonstrate the most significant role for top-down immune processes in bat cell lines described as constitutively antiviral in the literature. We further predicted that the most robust antiviral responses would (b) be correlated with the most rapid within-host virus propagation rates while (c) protecting cells against virus-induced mortality to support the longest enduring infections in tissue culture.”

And the beginning of our Discussion addresses these hypotheses:

“Critically, we found that – as hypothesized – (a) bat cell lines demonstrated a signature of enhanced interferon-mediated immune response, of either constitutive or induced form, which (b) allowed for establishment of rapid within-host, cell-to-cell virus transmission rates (β). These results were supported by both data-independent bifurcation analysis of our mean field theoretical model, as well as fitting of this model to viral infection time series established in bat cell culture. Additionally, we (c) demonstrated how the antiviral state induced by the interferon pathway protects live cells from mortality in tissue culture, resulting in in vitro epidemics of extended duration and an enhanced probability of establishing a long-term persistent infection.”

5) Add more explanation of what these virus lines are and why they are chosen. The reader is not told that rVSV-MARV and rVSV-EBOV are recombinant viruses expressing Ebola and Marburg Virus (Results first paragraph). Moreover, it is not stated why the authors employ two viruses (they do not seem to be interested in contrasting them but they cannot be said to be true replicates) or these two in particular.

Thank you for bringing this up. We’ve added to our Viruses section in the Materials and methods portion of the manuscript to further elaborate on and explain why the viruses were selected, as well as extended our Introduction to these time series upon first presentation in the Results section. We included both rVSV-EBOV and rVSV-MARV because we wanted to include a known constitutively antiviral cell line-virus combination (rVSV-MARV on PaKiT01 cells), over and above the constitutive expression of IFN-α in these cells.

6) How many replicates were there for each cell-line/clone combination? The only detail mentioned in this regard is that there were 2-3 technical replicates per plate. Such details should be stated prominently in the Results, since large numbers of points are plotted, but the reader has no way to discern which of them are independent.

We now highlight the number of replicates (18-39) in the Results section and the Materials and methods. Additionally, we include dashed statistical fits to each trial (see response to #7 below) in two figure supplements (Figure 1—figure supplement 2 and Figure 1—figure supplement 3) to further differentiate each time series.

7) There is a strong case for requiring the model is fitted to the individual replicates and not the mean of the replicates. If this case is rejected, then careful justification must be provided.

Thank you, and yes, we agree. In our newest version of the manuscript, we have refitted the mean field model to all individual replicates for each cell line-virus-MOI combination, rather than to the mean of replicates. We find our fits to be much improved and now consistent with our hypotheses.

8) The cell cultures are grown to 90% confluency, and are modeled as systems in which there is cell birth and natural death (birth rate about 4x the natural death rate). Readers are likely to be left with a number of questions that the authors should address more explicitly: Where does the birth rate parameter come from? Is this a realistic estimate reflective of the dynamics of the cell culture, or is it a deduction based on the requirement to sustain the infection dynamics (either endemic or oscillatory?). The Materials and methods suggest that it is what would be required to maintain sustainable live cell populations, but is this what actually happens in the culture during the infection? What sort of turnover does this result in, in the absence of infection? What is the justification for frequency over density dependent dynamics? At the heart of these questions is whether it is appropriate to think of these cell cultures as a sort of viral chemostat which could maintain steady state or oscillatory dynamics indefinitely, or whether the system would be better modeled as an infection process in a fixed non-reproducing host population. The latter view would substantially simplify the paper if such a perspective was justifiable.

Thanks for these valid and important insights and questions. We touched a bit on these questions in our responses to point #2 above, but we will elaborate further here. We did attempt to fit several versions of the model that incorporated only the natural mortality rate (fitted for each cell line) and the epidemic dynamics of each time series, with no birth rate included in the model. These versions of the model were unable to re-capture the extended time series witnessed in the rVSV-EBOV infections on PaKiT01 cell lines, using fitted mortality (μ) and incubation (σ) rates and a fixed, standard infection-induced mortality rate (α= 1/6). Thus, cell turnover was needed in order to recapture the data, or we would need to re-evaluate the infection-induced mortality parameter. In order to recapture our data, however, the duration of the infectious period this parameter would necessarily need to be extended to allow for cells to produce infectious virus for over a week without dying in order to recover the sustained 200 hour infectious time series witnessed in our data. This is very unlikely.

In our most recent version of the manuscript, in order to address the criticisms highlighted above, we chose to fit a unique natural mortality rate (μ) for each cell line, rather than fixing all μ at a constant and consistent value as done previously. In order to do this (as highlighted in the Figure 1—figure supplement 7) and discussed in the Materials and methods), we fixed the birth rate, b, at.025, then fit a μ for each cell line to recover the statistical mean of the declining susceptible count of all control well trials for that cell line. The net effect of these birth-death fits is still cell loss in tissue culture over time but with some cell turnover taking place. We feel confident that cell turnover did take place because we grew cells to 90% confluency prior to initiating the plaque assay, then watched them subsequently fill the plate when they could not be split at later timepoints. In including cell turnover, our analysis additionally maintains consistency with previous in vitromodeling of IFN dynamics in cell culture in the general literature (Howat et al., 2006).

The reviewers additionally questioned our decision to model these dynamics in a frequency vs. density-dependent context. We feel this assumption was well supported in the current context because plaque assays restrict transmission neighbor-to-neighbor in two dimensional space, thus limiting the number of possible infectious contacts for a single cell at all timesteps, in keeping with frequency-dependent dynamics. Intriguingly, because dead cells occupy space in our model, however, this has the effect of reducing the transmission rate in a density-dependent fashion: declines in susceptible and infectious cell density are realized dynamically as declines in their proportions where dead cells fill their place and do not contribute to epidemic dynamics. This effect makes the mean field model more generalizable to the level of the entire organism.

9) The potential influence of the assumed value for the rate at which antiviral cells regress to susceptibility is unclear. At present, it is not clear where the parameter value used comes from and what impact it could have on the results. This is particularly pertinent since sustained epidemics are ascribed to the antiviral cells returning to susceptibility – can one confidently separate this from birth of new susceptibles?

Agreed, and thank you for highlighting this potential point of controversy. Certainly, antiviral cells do have the capacity to regress to susceptible status (see Samuel and Knutson, 1982; Rasmussen and Farley, 1975; Radke et al., 1974), but after reconsideration of our assumptions and more intensive review of the literature, we agree with this reviewer critique that this is unlikely to occur on the 200 hour timescale of our tissue cultures, especially in such a limited spatial environment with live virus still present (and therefore IFN-inducing) in the tissue culture. As a result of this, we have fixed the value of c at 0 in all model fits (and bifurcation analyses) in the newest version of the manuscript, such that the regeneration of susceptibles to maintain late stage epidemic transmission can now only occur as a result of susceptible cell births. This greatly simplifies our model and also produces more robust model fits to the data. We have maintained the term c in our model equations because this antiviral regression phenomenon is relevant at the scale of the whole organism, and we wished to keep our derivations for R0 generalizable.

10) Why is ε not estimated? It seems to be fixed at 0 or 1. In subsection “Fitting of theoretical model to cell culture data demonstrates higher within-host transmission rates under constitutive immune assumptions” this seems to be recognized but Supplementary file 4 doesn't provide any intermediate values. It wasn't clear to why ε is bounded at 1. It is a rate?

Previously, ε was modeled as a 0-1 scaler that modified the rate of antiviral acquisition (ρ), which is why it was fixed. However, your comments inspired us to re-structure the model to estimate the induced immune rate of antiviral acquisition (ρ), which is multiplied by the proportion of exposed cells in tissue culture, separately from the constitutive rate of antiviral acquisition (ε), such that – to answer your questions – (a) ε is now estimated under models assuming some degree of constitutive immunity (it is fixed at 0 for induced-only and absent immune models) and (b) yes, ε is a rate (though it was not previously). We find this a much simpler way to present the model, allowing for greater flexibility in parameter estimation and more accurate fits to the data.

11) It is surely odd that as b (and u) is lowered towards zero, R0 also tends to zero. Surely one would expect a virus to be able to spread in a non-renewing monolayer?

In a simple, density-dependent model with no recovery from infection, R0 will be written as βN/(α + μ), indicating that the pathogen is unable to spread at low population densities. Though our transmission rates are cell-to-cell frequency-dependent, our model is really density -dependent as a result of the inclusion of dead cells as a state variable and, therefore, has cell density thresholds for pathogen invasion. In the supplementary material of this manuscript, we derive the equation for Ps at disease free equilibrium, which is equal to (b-μ)(c+μ)b(c+μ+ε). This configuration of Ps is incorporated into our equation for R0 because a pathogen will not be able to invade a tissue with too high an initial proportion of dead cells (PD) relative to live cells (PS). Under conditions of no constitutive immunity (ε=0), this equation for PSreduces to (b-μ)b which then surfaces in the equation for R0. Because our system of differential equations was constructed including birth rates, the birth and death rates are incorporated in the formulation of the susceptible proportion at DFE. Thus, a monolayer with a higher death rate than birth rate will be composed of completely dead cells at DFE and, therefore, not be invadable by a pathogen.

We hope that this explanation provides some clarity. If easier to understand, we could rewrite R0 with the term PSinstead of birth and death rates. Please let us know if you would prefer it this way.

Discussion:

12) The use of the word 'favor' in the first paragraph of the Discussion is a little misleading and should be re-worded. These findings are not based on any evolutionary process.

Very true. Thanks for this. We’ve dropped the word wherever it was used in the text.

13) The Discussion should make clear that the question of whether bats in general have constitutively active intracellular immunity is far from resolved. As touched on in the Introduction, this partly depends on definition – does constitutive expression of e.g. interferon-α actually provide constitutive protection? What about other immune genes which seem to be defective (e.g. Xie et al., 2018)? More generally, the immune systems of only a very small fraction of bat diversity (and of non-lab animals in general) has been studied, so it remains unclear how well we can generalise to all bats, or whether bats are truly unusual relative to other wild animals. Where possible the authors should list examples of studies showing constitutive immunity specifically in the species (or close relatives of these) thought to be reservoirs for the highly virulent viruses discussed.

Great point – and thank you! We agree that this issue is far from resolved, and, in fact, we believe that our results largely support the opposite conclusion: that constitutive IFN expression functions more like a primed, inducible defense than a constitutive, functional antiviral protein. We have emphasized this debate in the Introduction as well as discussed it in the Results and in the Discussion.

https://doi.org/10.7554/eLife.48401.sa2

Article and author information

Author details

  1. Cara E Brook

    1. Department of Integrative Biology, University of California, Berkeley, Berkeley, United States
    2. Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Visualization, Methodology
    For correspondence
    cbrook@berkeley.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4276-073X

  2. Mike Boots

    Department of Integrative Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Conceptualization, Supervision
    Competing interests
    No competing interests declared

  3. Kartik Chandran

    Department of Microbiology and Immunology, Albert Einstein College of Medicine, Bronx, United States
    Contribution
    Resources, Supervision, Funding acquisition
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0232-7077

  4. Andrew P Dobson

    Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Supervision
    Competing interests
    No competing interests declared

  5. Christian Drosten

    Institute of Virology, Charité-Universitätsmedizin Berlin, corporate member of Freie Universität Berlin, Humboldt-Universität zu Berlin, and Berlin Institute of Health, Berlin, Germany
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared

  6. Andrea L Graham

    Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    Contribution
    Conceptualization, Supervision, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6580-2755

  7. Bryan T Grenfell

    1. Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    2. Fogarty International Center, National Institutes of Health, Bethesda, United States
    Contribution
    Conceptualization, Supervision, Investigation
    Competing interests
    No competing interests declared

  8. Marcel A Müller

    1. Institute of Virology, Charité-Universitätsmedizin Berlin, corporate member of Freie Universität Berlin, Humboldt-Universität zu Berlin, and Berlin Institute of Health, Berlin, Germany
    2. Martsinovsky Institute of Medical Parasitology, Tropical and Vector Borne Diseases, Sechenov University, Moscow, Russian Federation
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2242-5117

  9. Melinda Ng

    Department of Microbiology and Immunology, Albert Einstein College of Medicine, Bronx, United States
    Contribution
    Conceptualization, Data curation, Validation, Investigation, Methodology
    Competing interests
    No competing interests declared

  10. Lin-Fa Wang

    Emerging Infectious Diseases Program, Duke-National University of Singapore Medical School, Singapore, Singapore
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2752-0535

  11. Anieke van Leeuwen

    1. Department of Ecology and Evolutionary Biology, Princeton University, Princeton, United States
    2. Royal Netherlands Institute for Sea Research, Department of Coastal Systems, and Utrecht University, Den Burg, Netherlands
    Contribution
    Conceptualization, Formal analysis, Supervision, Investigation, Visualization
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1987-1458

Funding

National Science Foundation (Graduate Research Fellowship)

  • Cara E Brook

Adolph C. and Mary Sprague Miller Institute for Basic Research in Science, University of California Berkeley (Postdoctoral Fellowship)

  • Cara E Brook

National Institutes of Health (R01-AI134824)

  • Kartik Chandran

Singapore National Research Foundation (NRF2012NRF-CRP001-056)

  • Lin-Fa Wang

Singapore National Research Foundation (NRF2016NRF-NSFC002-013)

  • Lin-Fa Wang

Deutsche Forschungsgemeinschaft (DR 772/10-2)

  • Christian Drosten

Bundesministerium für Bildung und Forschung (RAPID #01KI1723A)

  • Christian Drosten

Horizon 2020 (#653316)

  • Christian Drosten

DARPA (PREEMPT Program - Cooperative Agreement no. D18AC00031)

  • Cara E Brook

National Institutes of Health (1R01AI129822-01)

  • Cara E Brook

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

CEB was supported by a National Science Foundation Graduate Research Fellowship at Princeton University, a Miller Institute for Basic Research Fellowship at UC Berkeley, a DARPA PREEMPT program Cooperative Agreement grant D18AC00031, and an NIH grant 1R01AI129822-01. Tissue culture experiments were funded by an NIH grant R01 AI134824 to KC. Work in LFW’s lab was funded by the Singapore National Research Foundation grants (NRF2012NRF-CRP001-056 and NRF2016NRF-NSFC002-013). CD was supported by the German Research Council (DFG) grant DFG SPP1596 (DR 772/10–2), the Federal Ministry of Education and Research (BMBF) grant RAPID (#01KI1723A) and the EU Horizon 2020 grant EVAg (#653316). The authors thank the Chandran lab at Albert Einstein College of Medicine – in particular, Cecilia Harold, Megan Slough, Rohit Jangra, and Tanwee Alkutkar – for technical support during tissue culture experiments. The authors further thank Jessica Metcalf and the Graham lab at Princeton for conceptual guidance throughout the project’s development.

Senior Editor

  1. Neil M Ferguson, Imperial College London, United Kingdom

Reviewing Editor

  1. Dan Haydon, University of Glasgow, United Kingdom

Reviewer

  1. Dan Haydon, University of Glasgow, United Kingdom

Publication history

  1. Received: May 13, 2019
  2. Accepted: February 2, 2020
  3. Accepted Manuscript published: February 3, 2020 (version 1)
  4. Version of Record published: March 10, 2020 (version 2)
  5. Version of Record updated: March 12, 2020 (version 3)

Copyright

© 2020, Brook et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 34,821
    Page views
  • 3,669
    Downloads
  • 70
    Citations

Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Cara E Brook
  2. Mike Boots
  3. Kartik Chandran
  4. Andrew P Dobson
  5. Christian Drosten
  6. Andrea L Graham
  7. Bryan T Grenfell
  8. Marcel A Müller
  9. Melinda Ng
  10. Lin-Fa Wang
  11. Anieke van Leeuwen
(2020)
Accelerated viral dynamics in bat cell lines, with implications for zoonotic emergence
eLife 9:e48401.
https://doi.org/10.7554/eLife.48401

Categories and tags

Further reading

  1. Hear Cara Brook discuss why bats carry so many viruses deadly to other species

    Podcast

    From Ebola to COVID-19, why are bats so often the hosts of dreadful diseases?

    1. Ecology
    2. Evolutionary Biology

    Thermal phenotypic plasticity of pre- and post-copulatory male harm buffers sexual conflict in wild Drosophila melanogaster

    Claudia Londoño-Nieto, Roberto García-Roa ... Pau Carazo
    Research Article Updated

    Strong sexual selection frequently leads to sexual conflict and ensuing male harm, whereby males increase their reproductive success at the expense of harming females. Male harm is a widespread evolutionary phenomenon with a strong bearing on population viability. Thus, understanding how it unfolds in the wild is a current priority. Here, we sampled a wild Drosophila melanogaster population and studied male harm across the normal range of temperatures under which it reproduces optimally in nature by comparing female lifetime reproductive success and underlying male harm mechanisms under monogamy (i.e. low male competition/harm) vs. polyandry (i.e. high male competition/harm). While females had equal lifetime reproductive success across temperatures under monogamy, polyandry resulted in a maximum decrease of female fitness at 24°C (35%), reducing its impact at both 20°C (22%), and 28°C (10%). Furthermore, female fitness components and pre- (i.e. harassment) and post-copulatory (i.e. ejaculate toxicity) mechanisms of male harm were asymmetrically affected by temperature. At 20°C, male harassment of females was reduced, and polyandry accelerated female actuarial aging. In contrast, the effect of mating on female receptivity (a component of ejaculate toxicity) was affected at 28°C, where the mating costs for females decreased and polyandry mostly resulted in accelerated reproductive aging. We thus show that, across a natural thermal range, sexual conflict processes and their effects on female fitness components are plastic and complex. As a result, the net effect of male harm on overall population viability is likely to be lower than previously surmised. We discuss how such plasticity may affect selection, adaptation and, ultimately, evolutionary rescue under a warming climate.

    1. Ecology
    2. Evolutionary Biology

    Sensory conflict disrupts circadian rhythms in the sea anemone Nematostella vectensis

    Cory A Berger, Ann M Tarrant
    Research Article Updated

    Circadian clocks infer time of day by integrating information from cyclic environmental factors called zeitgebers, including light and temperature. Single zeitgebers entrain circadian rhythms, but few studies have addressed how multiple, simultaneous zeitgeber cycles interact to affect clock behavior. Misalignment between zeitgebers (‘sensory conflict’) can disrupt circadian rhythms, or alternatively clocks may privilege information from one zeitgeber over another. Here, we show that temperature cycles modulate circadian locomotor rhythms in Nematostella vectensis, a model system for cnidarian circadian biology. We conduct behavioral experiments across a comprehensive range of light and temperature cycles and find that Nematostella’s circadian behavior is disrupted by chronic misalignment between light and temperature, which involves disruption of the endogenous clock itself rather than a simple masking effect. Sensory conflict also disrupts the rhythmic transcriptome, with numerous genes losing rhythmic expression. However, many metabolic genes remained rhythmic and in-phase with temperature, and other genes even gained rhythmicity, implying that some rhythmic metabolic processes persist even when behavior is disrupted. Our results show that a cnidarian clock relies on information from light and temperature, rather than prioritizing one signal over the other. Although we identify limits to the clock’s ability to integrate conflicting sensory information, there is also a surprising robustness of behavioral and transcriptional rhythmicity.