The role of spatial structure in the evolution of viral innate immunity evasion: A diffusion-reaction cellular automaton model

Abstract

Most viruses have evolved strategies for preventing interferon (IFN) secretion and evading innate immunity. Recent work has shown that viral shutdown of IFN secretion can be viewed as a social trait, since the ability of a given virus to evade IFN-mediated immunity depends on the phenotype of neighbor viruses. Following this idea, we investigate the role of spatial structure in the evolution of innate immunity evasion. For this, we model IFN signaling and viral spread using a spatially explicit approximation that combines a diffusion-reaction model and cellular automaton. Our results indicate that the benefits of preventing IFN secretion for a virus are strongly determined by spatial structure through paracrine IFN signaling. Therefore, innate immunity evasion can evolve as a cooperative or even altruistic trait based on indirect fitness effects that IFN shutdown exerts on other members of the viral population. We identify key factors determining whether evasion from IFN-mediated immunity should evolve, such as population bottlenecks occurring during viral transmission, the relative speed of cellular infection and IFN secretion, and the diffusion properties of the medium.

About the Authors:

Ernesto Segredo-Otero

Roles Conceptualization, Investigation, Methodology, Software

Affiliation: Institute for Integrative Systems Biology (I2SysBio), Consejo Superior de Investigaciones Científicas-Universitat de València, València, Spain

Rafael Sanjuán

Roles Conceptualization, Funding acquisition, Investigation, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Institute for Integrative Systems Biology (I2SysBio), Consejo Superior de Investigaciones Científicas-Universitat de València, València, Spain

ORCID logo http://orcid.org/0000-0002-1844-545X

Introduction

Innate immune signaling during early infection constitutes an important line of defense against viruses. Infected cells recognize pathogen-associated molecular patterns and secrete type-I interferons (IFNs), which act in autocrine and paracrine manners to halt the spread of the infection locally [1, 2]. Autocrine signaling triggers antiviral responses in the secretor cell including gene expression arrest and apoptosis, whereas paracrine signaling induces a virus-resistant state in neighbor cells. Innate immunity imposes a strong selective pressure on viruses, which have consequently evolved a variety of evasion mechanisms at the level of viral sensing, signal transduction, and/or gene expression [2–4]. If successful, these mechanisms lead to avoidance or shutdown of IFN secretion.

We have recently shown that a mutant of vesicular stomatitis virus (VSV) that fails to block IFN secretion reduces the fitness of neighbor, IFN-blocking VSV by triggering local antiviral responses [5]. This interference defines IFN shutdown as a social trait, since the fitness of a given virus infecting a given cell (focal virus) depends on other members of the viral population (neighborhood). A neighborhood constituted by viruses that block or do not stimulate IFN provides a fitness advantage relative to being located near IFN-stimulating viruses. This fitness effect is indirect because it is determined by cells not infected by the focal virus, is exerted through paracrine signaling, and depends on the spatial distribution of infected cells and immunized cells.

Blocking IFN should also have a direct effect on the focal virus determined by autocrine signaling. In principle, since autocrine signaling is aimed at reducing viral progeny production, IFN blockade should directly benefit the focal virus. However, if the viral infection cycle proceeds faster than autocrine antiviral responses, this direct benefit may not be realized. In addition, directing viral proteins to prevent IFN secretion (by avoiding recognition or blocking production) might have a direct cost for the focal virus. For instance, in VSV, the matrix protein M prevents IFN secretion by inhibiting host gene expression [6], but M is also an essential structural component of the virion. Hence, using M proteins for blocking gene expression in the nucleus might be costly for virion morphogenesis. Another potential source of direct costs to IFN blockade is that anti-IFN viral proteins could downregulate cellular gene expression or trigger premature apoptosis, reducing the availability of cellular resources for viral progeny production. The existence of such costs is suggested by our previous results with VSV [5].

Therefore, IFN shutdown should have direct and indirect effects on viral fitness. Because indirect effects depend on other members of the viral and cellular populations, tackling the evolution of innate immunity evasion requires considering the spatial structure and the dynamics of viral spread and IFN-mediated responses. In general, the term spatial structure refers to a non-random arrangement of individuals in space, and can originate from physical barriers, limited dispersal, a tendency to aggregate, demographical history, and other processes related to population dynamics [7]. For instance, in viral infections, host-to-host transmission and intra-host dissemination is limited by anatomical barriers. Yet, even in the absence of such barriers, most viruses infecting solid tissues exhibit spatial structure, as evidenced by the formation of infection foci. Such foci can result from a simple diffusion-reaction process involving adsorption of viral particles to cells, production of viral progeny, and diffusion of progeny particles in the medium before reaching a new cell [8–11]. Analogously, the spatial structure of immune responses mediated by IFN and other cytokines can be determined by the diffusion of signaling molecules and the effect of these signals in receptor cells [12]. In addition to infection and immunity, diffusion-reaction models have been used for addressing a variety of eco-evolutionary questions such as, for instance, how bacterial CRISPR limits viral spread [13] and how spatial obstacles perturb population expansion fronts and promote random genetic drift [14].

Another useful approach to investigating viral spread is provided by cell automaton models, in which individual cells are simulated as elements of a grid, making spatial structure explicit [15–21]. This approach has been previously used for studying how the interaction between viruses and the immune system is influenced by space. For instance, it has been shown that spatial clustering of infected cells (which is a consequence of viral spread in foci) tends to reduce the ability of cytotoxic T-lymphocytes (CTLs) to clear the infection, but that clustering actually increases clearing efficacy if chemotaxis allows CTLs to perform a non-random search for infected cells [22]. Another work explored how spatial structure influences the evolution of CTL escape mutants in HIV-1 [23]. CTL escape mutant should be able to avoid CTL-mediated lysis regardless of the presence of non-mutant viruses in the same host because CTL-mediated lysis is epitope-specific. Nevertheless, CTLs can also secrete cytokines capable of repressing infection in an epitope-independent manner, for instance by reducing the overall susceptibility of cells to viral infection. CTL escape mutants might be sensitive to this mode of action and, consequently, their fitness might be dependent on the presence of CTL-stimulating viruses in the neighborhood.

Here, we investigate how natural selection determines the evolution of viral innate immunity evasion in the presence of spatial structure. For this, we use a cell automaton in which infection and IFN-mediated immunity occur as diffusion-reaction processes. From these simulations, we infer direct and indirect fitness effects associated to IFN shutdown, as well as quantitative descriptors of spatial structure. Within this framework, we analyze how viral demography, the timing of infection and immunity, and the physical properties of the virions and the medium should determine the evolution of IFN evasion.

Model

We considered NS cells susceptible to infection by a virus (uninfected), NE infected cells that are not yet virion-producers (eclipse phase), NP virus-producer cells, and ND cells killed by the virus (Fig 1A). Killed cells were not replaced by new cells and hence became equivalent to empty positions in the grid. We considered two virus variants, W and D. For simplicity, we ignored cells co-infected with both variants. The W variant blocked or did not stimulate IFN production (actively preventing IFN production or avoiding recognition by the innate immune system were equivalent in the context of this model). In contrast, cells infected with the D variant detected the virus and became IFN producers. The model allowed the relative speed of infection and IFN production to vary, such that eclipse-phase cells could be primed for IFN production (NEε) or be IFN producers (NEπ), and the same applied to D-producing cells (NPε and NPπ, respectively).

[Figure omitted. See PDF.]

Fig 1. Model.

A. Scheme of the infection and immunization processes. NS: susceptible (non-infected) cells. NE: infected cells in eclipse phase. NP: virion-producing cells. ND: cells killed by the virus. NR: immunized cells (blue). W-infected cells: red. D-infected cells: green. An eclipse phase between viral sensing (Nε) and IFN secretion (Nπ) is also considered. The relative speed of infection and IFN production can vary. Hence, D-infected cells can be in four possible stages (NEε, NPε, NEπ, NPπ). kv: virion infectivity (infection rate). kR: IFN immunization rate. τEP: viral eclipse half time. τεπ: half time between infection and IFN secretion. τPD: half time between virion production and cell death. Hence, τPD +τEP is the total duration of the infection cycle. rv(W) and rv(D): virion production rates for W and D, respectively. K = rvτPD is thus the number of virions produced per infected cell. δV: virion degradation/outflow rate. rβ: IFN production rate of immunized cells. δβ: IFN degradation/outflow rate. Values for each parameter are provided in Table 1. In this model, post-infection IFN effects are not allowed. S2 Fig shows a model with IFN-induced apoptosis. B. Fitness partition according to direct and indirect (neighborhood) components. For the W virus in a W neighborhood, fW|W = 0. The term c is the direct fitness effect of being an IFN-suppressor (c > 0: IFN suppression is costly; c < 0: IFN suppression is beneficial). For a D virus in an W neighborhood, fD|W = c. The term b is the indirect benefit of being in an IFN-suppressive (W) neighborhood, independent of direct effects. The derivation of the overall relative fitness of each variant (fW and fD) is shown in the text. rW quantifies how strongly W viruses are influenced by neighbors of the same type, and analogously for rD. The term r measures the spatial segregation between W and D in terms of infection and immunity. Three possible scenarios are shown (red: W-infected cells, green: D-infected, light blue: IFN-immunized cells). Left: full segregation (r = 1). Center: partial segregation (0 < r < 1). Right: complete mixing (r = 0).

https://doi.org/10.1371/journal.pcbi.1007656.g001

Based on known dynamics [24–26], though, we initially assumed that virion release preceded IFN secretion (NEπ ≈ 0). We also incorporated the observation that uninfected cells respond to IFN in a dose-dependent manner [5, 27] and become immunized (NR), but do not produce IFN themselves [28]. Infection (NS→NE) and immunization (NS→NR) were simulated as Poisson stochastic processes occurring with probability(1)for each cell and simulation time unit Δt (0.5 min). For infection, λ = kVVNS, where kV is the infection rate (infectivity), V the local virion concentration and NS equals 1 or 0. For immunization, λ = kRβNS, where kR is the immunization rate, β the local IFN concentration, and NS equals 1 or 0. All other cellular state transitions were modeled as random processes occurring with cumulative probability P = 0 for t < τ, P = 0.5 for t = τ, and P = 1 for t > τ, where τ is the half time of the corresponding transition (e.g. τEP for NE → NP; Fig 1A).

Viral dynamics obeyed a diffusion-reaction process in an orthogonally divided space, the partial derivative of the virion concentration with time in every grid position being(2)The first two terms correspond to the reaction part, where rv is the virion production rate of infected cells. For the D virus, the NP term was replaced with NPε + NPπ (that is, all D-producing cells regardless of IFN status). Wherein needed, we used different virion production rates for W and D, rv(W) and rv(D), respectively. The δV parameter is the virion degradation/outflow rate. We ignored the loss of virions due to adsorption to cells. The last term describes diffusion, where D is the Stokes-Einstein diffusion coefficient, Δx is the grid size and ΔV the virion concentration difference between grid subunits. To calculate it, we considered only adjacent cells for each grid position (x, y), such that:(3)In the limits of the grid, we used a continuous-system approach, meaning that the neighbors of the rightmost cells were the leftmost cells, and the neighbors of the uppermost cells were the bottommost cells. For reasons of computational efficiency simulations were performed using a time unit Δt = 0.5 min, but a finer time resolution was needed to calculate the diffusion process. We thus used shorter time units for the diffusion part only and assumed interim quasi-steady states for all other variables in the system.

Analogously, for IFN:(4)

Viral spread occurred exclusively through the above diffusion-reaction process (for illustration, S1 Fig shows the diffusion process of virions and IFN alone, without the reaction term). We did not consider other types of spread, such as cell-to-cell spread, active transport of viruses, or cell mobility. Cellular proliferation was also ignored, since the time scale of the infection was shorter than the typical cell division time of normal (non-tumoral) cells, with some exceptions (e.g activated lymphocytes). For simplicity, we also ignored details of the intracellular infection and immunization processes, and we assumed that, in each cell, after the eclipse phase both IFN and virion production was linear with time. In nature, viruses can follow different replication mechanisms [29] and, hence, virion release is not necessarily linear with time. Also, IFN production is regulated by positive feedback loops in the infected cell [30], which we did not consider.

Initially, we assumed that IFN had no effect on already infected cells but we later relaxed this assumption to include autocrine effects as well as post-infection paracrine effects. For this, we allowed infected cells to respond to IFN by undergoing apoptosis (S2 Fig; NA cells). We modeled IFN-triggered cell death (NE→NA, NP→NA) as a Poisson stochastic process as above, with λ = kAβNP for virion-producing cells and λ = kAβNE for cells in eclipse phase, where kA is the IFN-induced apoptosis rate of infected cells.

The cellular automaton diffusion-reaction simulations were performed using MATLAB R2018b scripts (S1 File).

We calculated the growth rate of each virus variant as(5)where NI(0) is the initial number of infected cells and NI the number of cells infected at time t including both producer and eclipse-phase cells (NI = NE + NP for the W variant; NI = NEε + NEπ + NPε + NPπ for the D variant). We used log10 instead of the natural logarithm to facilitate visualization of the results.

To obtain a time-integrated metric of viral fitness, we calculated , where N* is the cumulative number of cells infected throughout the simulation and tf the final time point (36 h).

Following our previous work [5], we defined the fitness (f) of a focal virus conditional to the presence of other viruses in the same population (neighborhood; Fig 1B). For this, we used the pure W infection (i.e. W virus in a W neighborhood) as reference and we expressed fitness relative to it, such that:(6)where RW|W is the growth rate of a pure W infection, or(7)for time-integrated fitness. Therefore, by definition pure W infections had fitness f = 0. Let us then consider the extreme situation in which, due to the presence of D neighbors, the spread of the W virus would be hampered by as many IFN-immunized cells as in a pure D infection. For this scenario, we defined the fitness of W as fW|D = −b (W virus in D neighborhood). Hence, −b is the cost of being in a D neighborhood and, reversely, b is the benefit of being in a W neighborhood. Although b depends on paracrine signaling (that is, the benefit of avoiding IFN secretion), it may also be indirectly influenced by other features of the neighborhood such as, for instance, local virion abundance, which determines competition for cellular resources. Next, let us consider a situation in which the spread of a D virus would not be hampered by the immunization of neighbor cells (D virus in W neighborhood). We defined the fitness of D infections accordingly as fD|W = c. Here, c can be interpreted as the direct effect of blocking IFN for the focal virus, independent of neighborhood. If c > 0, blocking IFN has a direct cost, whereas if c < 0 blocking IFN has a direct benefit. Finally, we defined the fitness of a D virus in a D neighborhood (pure D infection) as fD|D = c−b, which is simply the sum of the independent direct and indirect terms.

Having considered the four extreme cases in which the neighborhood was purely D or purely W, we then allowed for intermediate scenarios (Fig 1B). For this, we defined the overall fitness of W as(8)where rW is a parameter that quantifies how strongly W viruses are influenced by W neighbors. It follows that fW = −b(1−rW). Analogously, we defined fD = rDfD|D+(1−rD)fD|W, where rD quantifies how strongly D viruses are influenced by D neighbors. Thus, fD = c−rDb. Natural selection will favor the W variant and, therefore, will promote the evolution of IFN evasion if fW−fD>0. Consequently:(9)The quantity rW+rD−1 describes the difference between the neighborhoods of W and D. More precisely, rW+rD−1 is the difference between the indirect IFN-mediated effects experienced by W and D viruses. By denoting r = rW+rD−1, the condition for the evolution of IFN shutdown is:(10)This expression corresponds to Hamilton´s rule [31]. Hence, the condition for natural selection to favor viruses that block or avoid IFN signaling can be written in terms of classical social evolution theory. In the social evolution field, r has been often defined as the genetic relatedness between interacting partners for the relevant trait. In our context, such relatedness depends strictly on the spatial structure of infection and immunity. If the two variants are fully insulated, the benefits/costs of preventing/stimulating paracrine signaling will be felt only by viruses of the same kind as the focal virus and, hence, r = 1 (Fig 1B). In contrast, if the innate immune response triggered by the IFN-stimulating variant affects both variants equally, r = 0. If IFN-blocking and IFN-stimulating viruses are spatially segregated to some extent (r > 0), selection may indirectly favor IFN-blocking viruses because they take a greater share of the benefits of preventing paracrine signaling. Notice that all parameters considered (b, c, rW, rD, and r) are time-dependent because fitness and spatial structure vary as the infection and the immune response progress. Also, these parameters are defined at the population level, and are not inferred from an analysis of each individual cell in the grid but, instead, from total cell counts in a given population. An Excel spreadsheet for the calculation of these parameters is available upon request.

Results

Selection for IFN shutdown is determined by spatial structure

We investigated the dynamics of viral spread and innate immune responses in a square grid containing 47,961 cells. Our model considered virion and IFN diffusion as well as cell infection and immunization (Fig 1A). For this, we used the empirically-determined size and degradation rate of VSV particles and the measured diffusion coefficient of chicken IFN. Generally, parameter values corresponded to a rapidly growing lytic virus and adherent IFN-producing cells. Yet, the model is generally applicable to different types of viruses. Details of parameter values are provided in Table 1.

[Figure omitted. See PDF.]

Table 1. Default parameters used in the simulations.

https://doi.org/10.1371/journal.pcbi.1007656.t001

The simulated infections progressed as foci as a result of the diffusion-reaction process, reproducing the typical spread mode of many viruses (Fig 2A). In simulations containing only one type of virus variant, IFN-blocking (W) and IFN-stimulating (D) viruses initially spread at similar rates. However, D infections subsequently became halted by innate immunity, whereas W infections progressed and invaded the entire cell population. In mixed infections initiated with an equal input of W and D, the growth of both variants was halted as the immune response was deployed (Fig 2A). Hence, in mixed infections the fitness of W and D relative to a pure W infection (fW and fD, respectively) decreased with time.

[Figure omitted. See PDF.]

Fig 2. Dynamics of the infection-immunization process and inference of fitness components.

A. Left: spatial structure of the infection and immune response at different time points. Right: cell counts. Susceptible cells are shown in dark blue (cell counts not shown), cells infected with the W virus in red, cells infected with the D virus in green, immunized cells in light blue, and dead cells in white (counts now shown). Infections were initiated with 50 infected cells of pure W, pure D, or a 1:1 mix of W and D. The last row shows an unstructured infection (r = 0) in which the concentration of virions was homogenized at every time unit of the simulation. Parameter values are given in Table 1. B. Inferred fitness components. Top graph: fitness values of W and D in mixed infections relative to pure W infections (fW and fD). Middle: inferred direct (c) and indirect (b) fitness components. For pure D infections fD = fD|D = c−b whereas, in mixed infections lacking spatial structure, fW−fD = -c. This allows c and b to be inferred at each time point. In mixed infections exhibiting spatial structure fW−fD = rb−c, which gives the rb term. Bottom: parameters describing spatial structure. These were calculated as rW = 1+fW/b, rD = (c−fD)/b, and r = rW+rD−1. The plots of all parameters are averages obtained from 10 simulations. The time range shown encompasses from the start of IFN release (9 h) to the approximate time point in which the pure D infection died out (36 h, endpoint). The spatial structure parameters could not be measured reliably before 13 h because b was too close to 0. In these simulations, IFNs suppression had no direct fitness effects (c = 0).

https://doi.org/10.1371/journal.pcbi.1007656.g002

The fitness difference between W and D in mixed infections can be expressed as fW−fD = rb−c (see Model; Fig 1B) where b is the indirect benefit of being in a W, IFN suppressive neighborhood, c is the direct fitness effect of suppressing IFN on the focal virus, and r is the difference between the neighborhoods of W and D determined by spatial structure. In the absence of spatial structure, r = 0 and hence fW−fD = -c. Hence, to infer c we performed simulations of mixed infections in which spatial structure was disrupted by equalizing virion concentration in the grid at every time unit Δt (0.5 min). Knowing c, we obtained b by comparing the fitness of a pure D infection relative to a pure W infection, fD|D = c−b (Fig 2B). In these simulations, IFN had no effect on already infected cells and, hence, autocrine signaling was not allowed. Consequently, we obtained c = 0. In contrast, b increased with time as the cell population became immunized.

In mixed infections, initially r was close to 1 because each virus was only weakly influenced by neighbors from other infection foci. Yet, as the infection progressed and IFN-mediated immunity expanded throughout the cell population, r-values decreased (Fig 2B). We calculated the parameters describing spatial structure associated to each variant as rW = 1+fW/b, rD = (c−fD)/b, and r = rW+rD−1, where rW describes to what extent the fitness of W viruses is determined by neighbor viruses of the same type W and, analogously, rD describes to what extent the fitness of D viruses is determined by neighbor viruses of the same type D (see Model for details; Fig 1B). We found that rW decayed with time faster than rD, indicating that the neighborhood of both variants gradually became more similar to a pure D infection than to a pure W infection. This shows that IFN functions as a harmful diffusible molecule for the virus, akin to a pollutant or poison. Thus, we define IFN as a “public bad”, by contraposition to well-known diffusible public goods (e.g. microbial secreted enzymes).

We ended simulations at 36 h, a time point at which pure D infections died out, and calculated a time-integrated fitness value using the cumulative number of infected cells (N*), from which we obtained the corresponding time-integrated values of each parameter in our model (Table 2). Despite the inhibitory effects of IFN, W remained fitter than D. Because c = 0, this was strictly due to spatial structure which, in turn, emerged from the diffusion-reaction infection and immunity processes and allowed W to be on average less adversely affected by IFN-mediated immunity than D. These results suggest that spatial structure is a key determinant of the evolution of IFN suppression in viruses.

[Figure omitted. See PDF.]

Table 2. Fitness and spatial structure descriptors inferred from 10 replicate simulations using parameter values provided in Table 1.

https://doi.org/10.1371/journal.pcbi.1007656.t002

Factors promoting spatial insulation of IFN-suppressing and IFN-stimulating variants

As shown above, the fitness of IFN-suppressing virus variants depends on their spatial insulation from IFN-stimulating variants. A straightforward factor determining such segregation is the initial fraction of infected cells. To illustrate this, we performed simulations in which we varied the initial number of infected cells from NI(0) = 2 to NI(0) = 1000, keeping constant the initial frequency of each variant at 50% as well as the initial number of susceptible cells. We found that for NI(0) = 2, the W variant was largely unaffected by the presence of the D virus, whereas for NI(0) = 1000 the interference was such that both variants showed similarly low fitness (Fig 3). The relative fitness of the D variant also decayed as NI(0) increased due to self-interference (except for NI(0) = 1000, due to saturation of W spread). The indirect benefit of suppressing IFN secretion (b) tended to increase with NI(0) because larger NI(0) values produced more immunized cells in the absence of IFN suppression. However, spatial structure dropped from r ≈ 1 for NI(0) = 2 to r ≈ 0 for NI(0) = 1000. This effect was driven by a gradual populating of the W neighborhood with IFN-immunized cells, as shown by the decreasing rW values. Spatial structure determined the fitness advantage of W over D through the rb term, such that selection for IFN suppression was stronger for lower NI(0) values. These results suggest that IFN evasion should more easily evolve in viruses that experience severe population bottlenecks during transmission and hence produce few, isolated initial infection foci.

[Figure omitted. See PDF.]

Fig 3. Effect of the founder population size on the relative fitness of IFN supressing and IFN-stimulating virus variants.

Top: spatial structure of the infection at endpoint (36 h) for different numbers of initially infected cells, ranging from 2 to 1000. Pure W infections and mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, cells infected with the D virus in green, and immunized cells in light blue. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as a function of the initial number of infected cells (mean ± SEM values from ten replicate simulations are shown). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.g003

Another factor that should determine selection for IFN suppression is medium viscosity, since viscosity determines the diffusion coefficient of IFN and virions. We found that the spatial structure parameter r increased with medium viscosity. The more local nature of infection and immunization at higher viscosities allowed the W variant to be less adversely affected by D viruses, increasing the fitness difference between W and D and thus making selection for IFN suppression stronger (S3 Fig). We also set out to test whether selection for IFN suppression depended on virion size. For this, we performed simulations in which we varied virion size from 50 nm to 500 nm. We found that the spatial structure parameter r tended to increase with larger virion sizes, but that this effect was weak (S4 Fig) because IFN diffusion was unaffected.

IFN suppression as an altruistic trait

The above results indicate that IFN shutdown can be favored by natural selection acting on indirect fitness effects, as long as infection and immunity are spatially structured. This suggests that IFN shutdown could evolve even if it imposes a direct fitness cost to the actor (c > 0). From a social evolution perspective, this implies that innate immunity evasion could evolve as an altruistic trait. To illustrate this, we made IFN shutdown costly by assigning a direct advantage to D in terms of virion production rate (higher rv and thus, higher number of virions produced per cell). Initially, D expanded faster than W throughout the cell population. However, as cells became immunized, the spread of D was halted whereas W continued to expand (S5 Fig). Using the cumulative number of infected cells to calculate fitness as above, we found that W remained fitter than D for a virion production cost of up to twofold, i.e. rv(D) = 2rv(W) (Fig 4). As expected, c-values increased as the advantage of D in terms of virion productivity became larger. Less intuitively, b also increased despite the fact that parameters controlling paracrine innate immunity were not changed. This occurred because D neighborhoods contained more virions than W neighborhoods and hence experienced stronger local competition for cells. Hence, b captured fitness effects associated to spatial structure including but not limited to paracrine signaling. To show this, we repeated the above simulations switching innate immunity off (kR = 0). Under these conditions, D outcompeted W owing to its higher virion production rate. Yet, again, b increased as the excess virion production of D became larger, even if no paracrine response was possible, confirming that b captured the effect of local competition for cells (S6 Fig).

[Figure omitted. See PDF.]

Fig 4. Indirect selection of costly IFN suppression (altruism).

Top: spatial structure of the infection at endpoint (36 h) as a function of the cost associated to IFN suppression. This cost was implemented as an excess virion production rate for IFN-stimulating viruses compared to IFN-suppressing viruses. Pure W infections and mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, cells infected with the D virus in green, and immunized cells in light blue. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as a function of the cost of IFN suppression (mean ± SEM values from ten replicate simulations). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.g004

Contribution of autocrine signaling to direct fitness

To explore the situation in which IFN blockade provides a direct benefit (c < 0), we allowed for autocrine IFN signaling. To implement this, infected cells responded to IFN by undergoing apoptosis (S2 Fig). This also allowed infected cells to respond to IFN secreted by other cells (post-infection paracrine signaling). Without changing the rate at which IFN made non-infected cells immune (kR), we increased the rate at which IFN triggered apoptosis (kA). This showed that, even for relatively strong post-infection IFN effects, autocrine signaling had a minimal impact on the direct fitness component c, which was always close to zero and much smaller than the indirect components b or rb (S7 Fig). The reason for this is that virion production preceded IFN responses, such that IFN signaling on already infected cells was relatively unimportant. This is in line with our recent experimental results, which indicated a low efficacy of IFN (autocrine or paracrine) in already infected cells [5].

In contrast, when we allowed IFN to be released before virions (5 h versus 6 h) and IFN-triggered apoptosis was strong (kA>kR), b decreased and–c increased by the same amount (Fig 5). Hence, the indirect benefits of blocking IFN secretion partially became direct benefits. Despite this drop in b, the rb fitness component was largely insensitive to changes in kA. This occurred because, as kA increased, D infections became more localized and thus produced fewer immunized cells. This increased spatial structure and allowed the W virus to be less severely affected by the presence of D, as shown by the greater rW and r values. As a result, the rb and–c terms were similarly large, indicating that spatial structure-dependent indirect fitness effects contributed significantly to the benefits of IFN suppression even in the presence of strong post-infection IFN-induced apoptosis.

[Figure omitted. See PDF.]

Fig 5. Effect of IFN-mediated post-infection apoptosis on the fitness of IFN supressing and IFN-stimulating virus variants assuming that IFN secretion precedes virion release.

Top: spatial structure of the infection at endpoint (36 h) as a function of the rate of IFN-mediated apoptosis (kA) relative to IFN-mediated immunization (kR). Pure W infections and mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, cells infected with the D virus in green, and immunized cells in light blue. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as the ratio between the rate of IFN-mediated apoptosis and the rate of IFN-mediated immunization kA/kR (mean ± SEM values from ten replicate simulations). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.g005

Discussion

In previous work, we showed that IFN shutdown can be considered a social process because the IFN-stimulating status of a given virus modifies the fitness of neighbor viruses [5]. Specifically, mutants that stimulated IFN production had a negative impact on the fitness of IFN-stimulating and IFN-suppressing members of the viral population, whereas IFN blockers provided little or no benefits to IFN stimulators. Hence, from the virus perspective, IFN functions as a diffusible “public bad”, and the ability of natural selection to promote IFN shutdown depends critically on spatial structure. Here, we have investigated this process quantitatively using a biophysically well-defined system in which a diffusion-reaction model was implemented in cell automaton simulations. We have shown that spatial structure in the form of infection and immunization foci can emerge from IFN and virion concentration gradients, which in turn depend on demography (e.g. transmission bottlenecks) and the properties of the medium (e.g. viscosity). Because IFN is smaller and hence diffuses faster than virions, immunization can reach larger areas than infection, even if cells release virions faster than IFN. Delayed but effective innate immune responses capable of successfully arresting virus progression by producing a ring of immunized cells around infection foci have been experimentally shown previously using IFN-stimulating virus variants [24].

We found that, whenever virion release precedes IFN secretion, the fitness associated to IFN evasion should not be strongly determined by direct effects exerted on the infected cell, but rather by indirect (neighborhood) effects. An important consequence of indirect fitness effects being larger than direct effects is that IFN evasion could be favored by natural selection even if costly for the actor, meaning that IFN evasion could evolve as an altruistic trait. Cooperation and altruism have long been studied by evolutionary biologists [31, 36–40]. In principle, natural selection should not favor traits that diminish individual fitness because cooperators can be invaded by cheaters that benefit from cooperative traits without reciprocating. However, cooperation breakdown can be avoided if cooperators tend to interact preferentially with other cooperators. Spatial structure allows for such assortment in a wide variety of organisms and is hence considered a fundamental factor driving the evolution of cooperation [41–44]. Spatial structure should actually be the main driver of cooperation in simple organisms for which other mechanisms of assortment such as learning and memory are not possible, as is the case of viruses. Here, we have shown that particle diffusion provides a simple and general physical basis for cooperator assortment in viruses. Within the context of innate immunity evasion, in mixed populations containing both IFN-stimulating and IFN-suppressing virus variants, IFN and virion concentration gradients allow IFN suppressors to be less severely affected by immunization than IFN stimulators.

The role of spatial structure in the evolution of innate immunity evasion should be particularly critical when shutting down IFN entails a direct fitness cost, since IFN suppression then becomes an altruistic trait. In future work, it would be interesting to elucidate the direct and indirect fitness effects of IFN suppression in different viruses, and the mechanisms involved in these effects. Yet, we found that in the absence of such direct costs, and even when IFN suppression has direct benefits, the evolution of IFN shutdown should depend strongly on spatial structure, provided that the ability of IFN to halt viral spread is mainly determined by paracrine signaling. As we have shown, the importance of spatial structure increases if virion release precedes IFN signaling. Currently available empirical evidence supports the view that this is the case for different viruses [24, 25]. Although further work is required to establish more general conclusions, virion release prior to IFN signaling should be frequent in viral infections because rapid replication and progeny release tend to be critical for the survival of many viruses.

Spatial structure is pervasive in viral infections. In addition to the basic diffusion-reaction process investigated here, spatial structure can result from other processes such as localized cell-to-cell virus transfer, subversion of cellular structures for enhancing viral spread, use of mobile cells for disseminating the infection, limited virus trafficking among organs or body compartments, inter-host transmission bottlenecks, and so on [45–50]. In future work it would be interesting to investigate the validity of our results in these more complex frameworks. Other complex scenarios, such as tridimensional structures and the spread of the virus and cytokines through vessels could be considered. A common qualitative observation, though, is that viral infections often exhibit metapopulation structures in which quasi-isolated demes are initiated by small number of founder viral particles, increasing genetic drift and intra-deme genetic relatedness. This should make the evolution of cooperative traits such as IFN evasion more likely. In some cases, though, spatial structure might be less evident, such as in blood-borne viruses at the intra-host level. Potentially, selection for IFN evasion might be less efficient in these viruses, or these viruses may utilize different IFN evasion strategies, such as antagonizing downstream antiviral responses triggered by IFN signaling in infected cells.

Supporting information

S1 Fig. Diffusion model.

A. Diffusion of particles with a hydrodynamic diameter of 180 nm (virions, left) and 7 nm (IFN, right). B. Mean particle concentration with time (in arbitrary units). C. Variance of the particle concentration with time (green: 180 nm; blue: 7 nm; dashed lines: intermediate sizes).

https://doi.org/10.1371/journal.pcbi.1007656.s001

(TIF)

S2 Fig. Scheme of the infection and immunization processes in a model allowing for post-infection IFN effects.

NS: susceptible (non-infected) cells. NE: infected cells in eclipse phase. NP: virion-producing cells. ND: cells killed by the virus. NR: immunized cells (blue). W-infected cells: red. D-infected cells: green. An eclipse phase between viral sensing (Nε) and IFN secretion (Nπ) is also considered. The relative speed of infection and IFN production can vary. Hence, D-infected cells can be in four possible stages (NEε, NPε, NEπ, NPπ). kv: virion infectivity (infection rate). kR: IFN immunization rate. τEP: viral eclipse half time. τεπ: half time between infection and IFN secretion. τPD: half time between virion production and cell death. Thus τPD + τEP, is the total duration of the infection cycle. rv: virion production rate. K = rvτPD is thus the number of virions produced per infected cell. δV: virion degradation/outflow rate. rβ: IFN production rate of immunized cells. δβ: IFN degradation/outflow rate. Infected cells also respond to IFN by undergoing apoptosis (NA cells). kA: IFN-induced apoptosis rate of infected cells.

https://doi.org/10.1371/journal.pcbi.1007656.s002

(TIF)

S3 Fig. Effect of medium viscosity on the fitness of IFN supressing and IFN-stimulating virus variants.

Top: spatial structure of the infection at endpoint (36 h) for increasing viscosities. Pure W infections and mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, cells infected with the D virus in green, and immunized cells in light blue. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as a function of medium viscosity (mean ± SEM values from ten replicate simulations are shown). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.s003

(TIF)

S4 Fig. Effect of virion size on the fitness of IFN supressing and IFN-stimulating virus variants.

Top: spatial structure of the infection at endpoint (36 h) for increasing virion sizes. Pure W infections and mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, cells infected with the D virus in green, and immunized cells in light blue. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as a function of virion size (mean ± SEM values from ten replicate simulations are shown). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.s004

(TIF)

S5 Fig. Dynamics of the infection-immunization process and social evolution parameters with a cost to IFN blockade (c > 0).

a. Spatial structure of the infection and immune response at different time points (left) and cell counts (right). The color legend is as in Fig 2. Parameter values are as in Table 1, except that virus D produces twice as much progeny virions, i.e. rv(D) = 0.56 virions cell–1 min–1 and rv(W) = 0.28 virions cell–1 min–1. b. Inferred fitness components. The twofold excess progeny produced by D results in a direct fitness advantage c = 0.026 ± 0.001 that counterbalances the indirect fitness advantage obtained by W (see also Fig 4).

https://doi.org/10.1371/journal.pcbi.1007656.s005

(TIF)

S6 Fig. Effect of local competition on indirect fitness.

In these simulations, the D virus had a fitness advantage over the W virus in terms of an increased virion production rate, but innate immunity was disabled.Top: spatial structure of the infection at endpoint (36 h) as a function of the cost imposed to the W virus. Only mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, and cells infected with the D virus in green. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as a function of the cost imposed to W (mean ± SEM values from ten replicate simulations). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.s006

(TIF)

S7 Fig. Effect of the IFN-mediated post-infection apoptosis on the fitness of IFN supressing and IFN-stimulating virus variants assuming that virion release precedes IFN secretion.

Top: spatial structure of the infection at endpoint (36 h) as a function of the rate of IFN-mediated apoptosis (kA) relative to IFN-mediated immunization (kR). Pure W infections and mixed (W+D) infections are shown. Susceptible cells are shown in dark blue, cells infected with the W virus in red, cells infected with the D virus in green, and immunized cells in light blue. Infected cells are shown in a cumulative manner, meaning that all cells infected with each variant throughout the progression of the infection are depicted. Thus, dead cells are not shown. Bottom: time-integrated fitness components calculated using the cumulative number of cells infected with each variant in mixed infections, shown as a function of the ratio between the rate of IFN-mediated apoptosis and the rate of IFN-mediated immunization kA/kR (mean ± SEM values from ten replicate simulations). Left: fitness of W and D relative to pure W infections. Center: direct and indirect fitness components. Right: descriptors of spatial structure.

https://doi.org/10.1371/journal.pcbi.1007656.s007

(TIF)

S1 File. A .rar file containing the MATLAB scripts that perform the simulations and associated files.

https://doi.org/10.1371/journal.pcbi.1007656.s008

(RAR)

Acknowledgments

We thank Pau Carazo and Stuart West for helpful discussions.

Citation: Segredo-Otero E, Sanjuán R (2020) The role of spatial structure in the evolution of viral innate immunity evasion: A diffusion-reaction cellular automaton model. PLoS Comput Biol 16(2): e1007656. https://doi.org/10.1371/journal.pcbi.1007656

References

1. Ivashkiv LB, Donlin LT. Regulation of type I interferon responses. Nat Rev Immunol. 2014; 14: 36–49. pmid:24362405

2. Fensterl V, Chattopadhyay S, Sen GC. No love lost between viruses and interferons. Annu Rev Virol. 2015; 2: 549–72. pmid:26958928

3. García-Sastre A. Ten strategies of interferon evasion by viruses. Cell Host Microbe. 2017; 22: 176–84. pmid:28799903

4. Coccia EM, Battistini A. Early IFN type I response: Learning from microbial evasion strategies. Semin Immunol. 2015; 27: 85–101. pmid:25869307

5. Domingo-Calap P, Segredo-Otero EA, Durán-Moreno M, Sanjuán R. Social evolution of innate immunity evasion in a virus. Nat Microbiol. 2019; 4: 1006–13. pmid:30833734

6. Rajani KR, Pettit Kneller EL, McKenzie MO, Horita DA, Chou JW, Lyles DS. Complexes of vesicular stomatitis virus matrix protein with host Rae1 and Nup98 involved in inhibition of host transcription. PLoS Pathog. 2012; 8: e1002929. pmid:23028327

7. Yanni D, Márquez-Zacarias P, Yunker PJ, Ratcliff WC. Drivers of spatial structure in social microbial communities. Curr Biol. 2019; 29: R545–R50. pmid:31163168

8. Yin J, McCaskill JS. Replication of viruses in a growing plaque: a reaction-diffusion model. Biophys J. 1992; 61: 1540–9. pmid:1617137

9. You L, Yin J. Amplification and spread of viruses in a growing plaque. J Theor Biol. 1999; 200: 365–73. pmid:10525396

10. Smith HL, Thieme HR. Persistence of bacteria and phages in a chemostat. J Math Biol. 2012; 64: 951–79. pmid:21656281

11. Amor D, Fort J. Cohabitation reaction-diffusion model for virus focal infections. Physica A. 2014; 416: 9.

12. Oyler-Yaniv A, Oyler-Yaniv J, Whitlock BM, Liu Z, Germain RN, Huse M, et al. A tunable diffusion-consumption mechanism of cytokine propagation enables plasticity in cell-to-cell communication in the immune system. Immunity. 2017; 46: 609–20. pmid:28389069

13. Payne P, Geyrhofer L, Barton NH, Bollback JP. CRISPR-based herd immunity can limit phage epidemics in bacterial populations. Elife. 2018; 7: e32035. pmid:29521625

14. Mobius W, Murray AW, Nelson DR. How obstacles perturb population fronts and alter their genetic structure. PLoS Comput Biol. 2015; 11: e1004615. pmid:26696601

15. Sato K, Matsuda H, Sasaki A. Pathogen invasion and host extinction in lattice structured populations. J Math Biol. 1994; 32: 251–68. pmid:8182357

16. Wodarz D, Hofacre A, Lau JW, Sun Z, Fan H, Komarova NL. Complex spatial dynamics of oncolytic viruses in vitro: mathematical and experimental approaches. PLoS Comput Biol. 2012; 8: e1002547. pmid:22719239

17. Akpinar F, Inankur B, Yin J. Spatial-temporal patterns of viral amplification and interference initiated by a single infected cell. J Virol. 2016; 90: 7552–66. pmid:27279621

18. Capitán JA, Cuesta JA, Manrubia SC, Aguirre J. Severe hindrance of viral infection propagation in spatially extended hosts. PLoS One. 2011; 6: e23358. pmid:21912595

19. Wodarz D, Levy DN. Effect of different modes of viral spread on the dynamics of multiply infected cells in human immunodeficiency virus infection. J R Soc Interface. 2011; 8: 289–300. pmid:20659927

20. Xiao X, Shao SH, Chou KC. A probability cellular automaton model for hepatitis B viral infections. Biochem Biophys Res Commun. 2006; 342: 605–10. pmid:16487938

21. Strain MC, Richman DD, Wong JK, Levine H. Spatiotemporal dynamics of HIV propagation. J Theor Biol. 2002; 218: 85–96. pmid:12297072

22. Kadolsky UD, Yates AJ. How is the effectiveness of immune surveillance impacted by the spatial distribution of spreading infections? Philos Trans R Soc Lond B Biol Sci. 2015; 370.

23. Seich Al Basatena NK, Chatzimichalis K, Graw F, Frost SD, Regoes RR, Asquith B. Can non-lytic CD8+ T cells drive HIV-1 escape? PLoS Pathog. 2013; 9: e1003656. pmid:24244151

24. Voigt EA, Swick A, Yin J. Rapid induction and persistence of paracrine-induced cellular antiviral states arrest viral infection spread in A549 cells. Virology. 2016; 496: 59–66. pmid:27254596

25. Howat TJ, Barreca C, O'Hare P, Gog JR, Grenfell BT. Modelling dynamics of the type I interferon response to in vitro viral infection. J R Soc Interface. 2006; 3: 699–709. pmid:16971338

26. Voigt EA, Yin J. Kinetic differences and synergistic antiviral effects between type I and type III interferon signaling indicate pathway independence. J Interferon Cytokine Res. 2015; 35: 734–47. pmid:25938799

27. Qi H, Chu V, Wu NC, Chen Z, Truong S, Brar G, et al. Systematic identification of anti-interferon function on hepatitis C virus genome reveals p7 as an immune evasion protein. Proc Natl Acad Sci USA. 2017; 114: 2018–23. pmid:28159892

28. Levy DE, Marie I, Smith E, Prakash A. Enhancement and diversification of IFN induction by IRF-7-mediated positive feedback. J Interferon Cytokine Res. 2002; 22: 87–93. pmid:11846979

29. Sanjuán R, Domingo-Calap P. Mechanisms of viral mutation. Cell Mol Life Sci. 2016; 73: 4433–48. pmid:27392606

30. Michalska A, Blaszczyk K, Wesoly J, Bluyssen HAR. A positive feedback amplifier circuit that regulates interferon (IFN)-stimulated gene expression and controls type I and type II IFN responses. Front Immunol. 2018; 9: 1135. pmid:29892288

31. Gardner A, West SA, Wild G. The genetical theory of kin selection. J Evol Biol. 2011; 24: 1020–43. pmid:21371156

32. Cureton DK, Massol RH, Whelan SP, Kirchhausen T. The length of vesicular stomatitis virus particles dictates a need for actin assembly during clathrin-dependent endocytosis. PLoS Pathog. 2010; 6: e1001127. pmid:20941355

33. Kreuz LE, Levy AH. Physical properties of chick interferon. J Bacteriol. 1965; 89: 462–9. pmid:14255715

34. Klasse PJ. Molecular determinants of the ratio of inert to infectious virus particles. Prog Mol Biol Transl Sci. 2015; 129: 285–326. pmid:25595808

35. Zimmer B, Summermatter K, Zimmer G. Stability and inactivation of vesicular stomatitis virus, a prototype rhabdovirus. Vet Microbiol. 2013; 162: 78–84. pmid:22995872

36. Nowak MA. Five rules for the evolution of cooperation. Science. 2006; 314: 1560–3. pmid:17158317

37. Fletcher JA, Doebeli M. A simple and general explanation for the evolution of altruism. Proc Biol Sci. 2009; 276: 13–9. pmid:18765343

38. West SA, Griffin AS, Gardner A, Diggle SP. Social evolution theory for microorganisms. Nat Rev Microbiol. 2006; 4: 597–607. pmid:16845430

39. Bourke AF. Hamilton's rule and the causes of social evolution. Philos Trans R Soc Lond B Biol Sci. 2014; 369: 20130362. pmid:24686934

40. Damore JA, Gore J. Understanding microbial cooperation. J Theor Biol. 2012; 299: 31–41. pmid:21419783

41. Nowak MA, May RM. Evolutionary games and spatial chaos. Nature. 1992; 359: 4.

42. Nowak MA, Tarnita CE, Antal T. Evolutionary dynamics in structured populations. Philos Trans R Soc Lond B Biol Sci. 2010; 365: 19–30. pmid:20008382

43. Nadell CD, Drescher K, Foster KR. Spatial structure, cooperation and competition in biofilms. Nat Rev Microbiol. 2016; 14: 589–600. pmid:27452230

44. Momeni B, Waite AJ, Shou W. Spatial self-organization favors heterotypic cooperation over cheating. Elife. 2013; 2: e00960. pmid:24220506

45. Zwart MP, Elena SF. Matters of size: genetic bottlenecks in virus infection and their potential impact on evolution. Annu Rev Virol. 2015; 2: 161–79. pmid:26958911

46. Gutierrez S, Michalakis Y, Blanc S. Virus population bottlenecks during within-host progression and host-to-host transmission. Curr Opin Virol. 2012; 2: 546–55. pmid:22921636

47. McCrone JT, Lauring AS. Genetic bottlenecks in intraspecies virus transmission. Curr Opin Virol. 2018; 28: 20–5. pmid:29107838

48. Pfeiffer JK, Kirkegaard K. Bottleneck-mediated quasispecies restriction during spread of an RNA virus from inoculation site to brain. Proc Natl Acad Sci USA. 2006; 103: 5520–5. pmid:16567621

49. Murooka TT, Deruaz M, Marangoni F, Vrbanac VD, Seung E, von Andrian UH, et al. HIV-infected T cells are migratory vehicles for viral dissemination. Nature. 2012; 490: 283–7. pmid:22854780

50. Doceul V, Hollinshead M, van der Linden L, Smith GL. Repulsion of superinfecting virions: a mechanism for rapid virus spread. Science. 2010; 327: 873–6. pmid:20093437