Introduction

Adaptive immune responses are relatively slow since they require pathogen-specific priming of immune cells (Sette and Crotty, 2021). For example, the time required for the body to activate adaptive immunity against the SARS-CoV-2 virus upon initial infection is around 10 days, comparable to the delay of immunization against SARS-CoV-2 after vaccination (Polack et al., 2020). Instead, the earliest infection dynamics are largely governed locally, by infected cells and their neighborhood. The innate responses including both interferon (IFN) mediated intercellular communication and expression of antiviral genes (ISGs) are determinants for confining the viral spread in the respiratory tract. Here, we address the spread of viruses within epithelial tissue, using SARS-CoV-2 as a model pathogen. The overall considerations are similar for other viruses, but the parameters governing infection may vary considerably due to the specific countermeasures of the virus in question, affecting its ability to bypass human antiviral defenses.

In terms of countermeasures, insufficient type I and III interferon secretion upon infection is a main immune signature feature of SARS-CoV-2 infection (Blanco-Melo et al., 2020; Hatton et al., 2021; Stanifer et al., 2020; Minkoff and tenOever, 2023). The failure to activate immediate antiviral responses with IFNs is also a pathogenic aspect of other viruses including Ebola (Mohamadzadeh et al., 2007), Marburg (He et al., 2019) and Herpes simplex (Barreca and O’Hare, 2004). Secretion of IFN relies on the cell’s ability to sense viral products during its replication. Despite the presence of sensors for DNA and RNA viruses in cells, many species of viruses partially evade detection. The SARS-CoV-2 virus is such a case: Only two of 16 putative RNA virus sensors, IFIH1 (MDA5) and DHX58 (LGP2) from the RIG-I-like receptor (RLR) family, play roles in inducing IFN upon SARS-CoV-2 infection (Yin et al., 2021) and IFIH1 is antagonized by SARS-CoV-2 (Liu et al., 2021).

Intriguingly, evidence shows that pre-activated innate immune states help combat the SARS-CoV-2 infection. The higher basal expression of viral sensors, IFIH1 and DDX58 (also from the RLR family), in the upper airway of children (relative to adults), reduces the severity of COVID (Loske et al., 2022). Furthermore, well-differentiated primary nasal epithelial cells derived from a donor with pre-activated IFN_γ do show resistance to SARS-CoV-2 infection (Broadbent et al., 2022). Thus, the extent to which innate immunity contributes to the observed heterogeneity in responses to SARS-CoV-2 between hosts (Schaller et al., 2021; Desai et al., 2020) is a compelling subject for investigation.

To address this question, we reanalyze single-cell RNAseq data (Fiege et al., 2021; Ravindra et al., 2021) providing gene expression profiles of virus sensors and antiviral genes in host cells during early SARS-CoV-2 infection. We propose a cellular automaton model based on a few transition rules suggested by observed cell states, to explain the heterogeneity in early disease progression as a consequence of criticality in the virus-host interaction system.

Cell states during infection

Directly observing cell responses and cell state transitions in a patient’s body upon viral infection is virtually impossible. However, human bronchial epithelial cells (HBECs) mimic the airway epithelium and have been used as a representative model for investigating the consequences of the viral invasion (de Jong et al., 1993, 1994; Davis et al., 2015). Single-cell RNAseq provides snapshots of the states of individual cells indicated by high-dimensional gene expression profiles at the mRNA level and can uncover the heterogeneity of cell responses obscured by aggregate measurement. Thus, by combining HBECs as a model and single-cell RNAseq data, one can in principle infer cell state transitions following viral infection. More importantly, single-cell RNAseq also captures copies of viral genes during sequencing, which allows us to simultaneously estimate viral replication inside cells.

To reconstruct the trajectory of cell state transitions during early SARS-CoV-2 infection, we re-analyze single-cell RNAseq data from experiments where HBECs are sampled before infection (0 h), as well as 24 and 48 hours post-viral infection (hpi) (Fiege et al., 2021). We focus on genes associated with antiviral responses and interferons from the host cells and detected viral genes. We project high-dimensional gene expression data onto a 2D plane using Uniform Manifold Approximation and Projection (UMAP) and obtain a low-dimensional visualization of single-cell expression patterns (Fig. 1a). On the UMAP plane, each dot represents a cell sample and the distance between dots correlates with the level of similarity of cellular states. The cells are not divided absolutely into discrete clusters and rather show continuous trajectories. For simplicity and convenience, we cluster the cells according to their gene expression.

Figure 1.

Different clusters on the UMAP indicate distinct cellular states during the progression of infection. For instance, there are three sub-clusters of susceptible cells (O₁, O₂, O₃). Neither viral genes nor IFNs are detected in these cells and only a few antiviral genes are expressed. The viral sensors (DHX58, DDX58, and IFIH1) are at their lowest level (Fig. 1b, Fig. S1). We refer to all of these cells as O cells due to their relatively similar gene expression profiles in terms of viral replication genes. The proportion of O cells decreases over time as the infection spreads (Fig. 1c).

We also observe three infected cell clusters where viral genes are primarily detected, V ⁱ, V ^r, and V . With the increasing counts of viral genes, we infer that the V ⁱ cluster is the earliest state after an O cell has been infected and the virus begins to replicate. Some but not all antiviral genes are activated in the V ⁱ cells (IFIT1/2/3 and OAS1/2/3) (Fig. 1b and Fig. S1), indicating that these cells are still vulnerable to viral invasion. This cluster is followed by two subsequent clusters, the V ^r cluster with pronounced viral replication and A cluster with barely any viral replication.

In the V cluster, the viral genes reach their highest level and antiviral genes are strongly inhibited, indicating that the virus has fully hijacked the cell. The antiviral genes are expressed most strongly in the A cluster and partially in the N cluster, indicating that the antiviral capability of the N cluster is weaker than the full antiviral state. Although the N cluster also shows a high level of viral genes, it severely lacks one of the viral genes (cov.E, Fig. S2) compared with the most highly expressed viral genes of the V cluster. This observation implies that viral replication and activation of the antiviral state coexist in the IFN-secreting cells (N cluster). We note the existence of a small subgroup of the V ^r cluster, close to the A cluster, that exhibits relatively high levels of both antiviral genes and viral genes but no appreciable IFN (Fig. 1d–f). The viral genes are also partially expressed in these cells, but different from the N cluster, the antiviral genes are fully expressed (Fig. S1 and S2). Thus, these cells are more likely to sustain the antiviral state.

At 24 hpi, some cells have switched from the pre-infection state (O) to other states. At 48 hpi, almost all cells have transitioned to other states and only a few cells remain in the O state (Fig. 1c and g). The aggregated gene expression of representative antiviral genes and detected viral genes indicates the cells move from the O state towards the three remaining terminal states on the considered timescale of 2 days: Antiviral state (A, Fig. 1d), Virus-conquered state (V, Fig. 1e), and IFN producing state (N, Fig. 1f). Central for the overall defense is the relatively few cells that reach the IFN-producing state (N ). These cells also express A and V genes.

When IFN is not expressed, the antiviral genes and viral genes exclude each other (Fig. 1d and e), except for a few cells around (UMAP1, UMAP2) ∼ (−2.5, 7.5) (green cells at hpi = 48, Fig. 1g). They represent cells where the virus succeeded in stopping IFN secretion, but could not fully hijack the cell. We still regard these cells as antiviral cells in our model.

The N state is associated with both high levels of virus sensors and viral genes, in agreement with the observation that IFN production is initiated after exposure to the virus (Lei et al., 2020) and that IFN can induce an antiviral state inside the same cell (Sanceau et al., 1987). Expression of the key SARS-CoV-2 sensitive sensors (IFIH1, DDX58, DHX58) is sparse in the O state (Fig. S1), indicating that a small fraction of cells have virus-sensing capacity prior to infection and are ready to mount a defense – and that this population increases with IFN tissue diffusion.

Model

We introduce a cellular automaton model to capture the cell state dynamics during early stages of SARS-CoV-2 infection in a sheet of epithelial tissue. At each simulation, we seed an infection site on a 2D square lattice and study how the infection spreads as the sites on the lattice switch between cell states following a set of simple rules derived from the observations of the single-cell RNAseq data.

In addition to the states corresponding to the dominant clusters observed in the single-cell data (Fig. 1a) (O, A, V and N states corresponding to O, A, V and N clusters), we introduce a transient pre-antiviral state (a) that can switch to the N state rapidly on viral exposure, considering the heterogeneity of viral sensing ability in susceptible cells.

It follows from this description, that those RNA viruses that can sensed by a large repertoire of sensors should be modeled with a larger fraction of cells in the a state.

The model is initialized with cells predominantly in the O state and a small fraction, P_a, in the pre-antiviral state a. The parameter P_a can also be understood as the probability that an O cell will switch to the N or A state when exposed to the virus or IFNs, respectively. As such, the value of P_a depends on both host and virus. In particular, a virus that is able to effectively interfere with the defense and signaling of host cells will be modeled by a low P_a value.

It is worth noting that the proportion of cells in the a state before the onset of SARS-CoV-2 infection is expected to be higher in hosts with pre-activated antiviral innate immunity (Loske et al., 2022; Broadbent et al., 2022), meaning that the value of P_a will, in general, depend on the exposure history of the host.

We propose the following transitions between five main discrete cell states (Fig. 2a):

• N, IFN-secreting cells. These arise from pre-antiviral cells (a state) that become infected (but not infectious). Here we ignore possible apoptosis (Wen et al., 1997; Tesfaigzi, 2006) of N cells since IFNs still have time to diffuse to the neighborhood.

• O, unaffected (susceptible) cells.

• V, infected cells, virus-producing cells. This state arises when a susceptible (O) cell is exposed to a virus from another V cell.

• a, pre-antiviral state. May develop into either the A or N state upon exposure to signals from N cells or virus from V cells.

• A, antiviral state. Immune to infection. Achieved when a pre-antiviral (a) cell is exposed to IFN. We do not consider the decay of the antiviral state as it may last more than 72 h (Gaajetaan et al., 2013).

Figure 2.

The dynamics are defined in terms of discrete time steps, representing the characteristic timescales of cellular viral infection. We explore the model for an extended time, keeping in mind that in reality other immune cells such as NK-killer cells and macrophages may migrate to the infected site and reduce viral spread (McNab et al., 2015).

The four rules of the model are (Fig. 2a):

where the notation X(Y ) = Z denotes a cell in state X acting on a cell in state Y and changing it to state Z in one time-step. Thus, cells in states O, a, and A are unable to influence their neighbors. The V state is the only directly self-replicating state.

Each site of the L × L lattice is assigned to either the O (probability: 1 − P_a) or the a state (probability: P_a). Infection is initiated by a single V cell, and we explore the percolation of the infection to larger scales. A time step consists of L² updates, in which a random site i is selected. If a V cell is selected, it interacts with its 4 nearest neighbors according to the rules V (O) = V and V (a) = N . If an N cell is selected, it interacts with all cells within a radius R, according to the rules N (O) = a and N (a) = A. The radius R thus quantifies the diffusion range of IFNs relative to the virus.

Results

At R = 1, the final number of infected cells depends strongly on the value of P_a. At a low P_a of 0.27, infections typically spread to the entire system, while at a higher P_a of 0.35, the propagation of the V state is inhibited (Fig. 2b).

We observe a threshold-like behavior of the final attack rate of the virus when the initial P_a changes continuously (Fig 2c). The virus spreads macroscopically for P_a < P_c ≈ 27.8%. At higher P_a, cells are sufficiently prone to convert to the antiviral state to prevent the infection from percolating.

The size distribution P (s) of infection clusters around the critical value of P_a = 27.8% obeys power-law decay (Fig. 3a). In Fig. 3b, the distribution P (s) ∝ 1/s^τ is further explored by re-scaling and the cluster size exponent is confirmed as τ = 1.83 ± 0.03 when P_a = 28%. Notably, this exponent is below the equilibrium 2D percolation yielding τ = 2.05 (Stauffer and Aharony, 2018). Further, our exponent τ ∼ 1.8 is above to the percolation-inspired cluster growth model for virus spread of ref. (Gönci et al., 2010) which has an exponent between τ = 1.58 and τ = 1.64 depending on the distribution of individual cells’ pre-defined ability to become infected.

Figure 3.

The actual critical value of P_a depends strongly on the choice of neighborhood. In particular, at R = 1, the V and N states have the same range in the tissue (a proxy for diffusivity), while a more realistic scenario is to allow IFNs to diffuse faster in the tissue (R > 1), facilitating the initiation of the antiviral state. The critical percolation threshold P_c decreases almost exponentially with the value of R (Fig. 4b), and viral propagation can be stopped for P_a as low as 0.4% when R ≥ 5. Such a small fraction of initial cells in the a state is consistent with the remarkably few N cells observed in experiments (Fig. 1b). Thus, a higher diffusivity of IFN provides a more than proportional decrease in the required number of antiviral cells. As revealed by the reanalysis of RNAseq data in Fig. 1, the fraction of IFN-positive cells is relatively low – around 1.7%. Comparing with simulations near the critical point, we find that, at R = 5, the ratio of N cells to all affected cells (N + A + V ) in the final state, lim_t→∞ N /(N + A + V ) ≈ 2%, i.e. it is of comparable magnitude to the experimental value. This holds in a wide range around the critical point, P_a ∼ P_c.

Figure 4.

The exponents for the cluster size distribution are the same at R = 1 and R = 5, while the structures of the clusters are different (Fig. 4a and c). Greater R leads to a different microscopic structure with fewer A and N cells in the final state (Fig. 4c).

To put the above findings in perspective we further explore a simplified version of our model with only 3 states (supplementary Fig. S3), the OVA model, which may be seen as a rephrasing of models for induced antiviral states suggested by Howat et al. (2006); Segredo-Otero and Sanjuán (2020); Michael Lavigne et al. (2021). In the OVA model, P_a is the probability that an infected cell produces interferons to warn neighbor cells within radius R. In the OVA model, one update consists of selecting a random cell. If the cell is in the V state then its neighbor cells may change by exposure to the virus, provided that they are susceptible (O). Each of the four neighbors is now chosen in random order, and if a neighbor cell i is in the O state, a random number ran_i ∈ [0, 1] is drawn. If ran_i ≥ P_a the neighbor is flipped to the V state. If, on the other hand, ran_i < P_a, all O cells within a radius R around the neighbor i are converted to the A state. Thus, for large R and moderate P_a, the spread of infection will be mitigated. We find that the OVA model has an “outbreak size” exponent τ ∼ 1.8, similar to the NOVAa model. However, the change in microstructure as a function of the IFN range R observed in the NOVAa model (compare Fig. 4 panel a and c) is not observed in the OVA model (Fig. S3), where the features instead scale proportionally with R. We also simulated standard percolation by randomly adding disks of radius R of blocking (“antiviral”) cells and checking for percolation of the infected state. While the critical behavior of the standard percolation model approximately resembles that of the OVA model (Fig. 4b), the antiviral state of the OVA model is somewhat less effective at blocking the spread (reflected in a higher threshold P_c).

Finally, we examined a version of the model where the discrete idealization of N -cells acting at all cells within a specific radius R is replaced by a probabilistic conversion with a diffusion-like profile. The algorithm for this is described in the supplement, with results in Fig. S4 to be compared to Fig. 4. We find that the probabilistic spreading of IFN is more effective, in terms of demanding lower R for obtaining similar limits on spreading of the infection.

Discussion

There are some preexisting models of the interplay between virus, host cells and triggered immune responses, with an antiviral state triggered by IFN signaling from neighbor cells (Graw and Perelson, 2016). Cellular automaton models of infection dynamics in epithelial tissue were explored by (Howat et al., 2006; Segredo-Otero and Sanjuán, 2020; Michael Lavigne et al., 2021), with the overall result that spreading depends on competition between the virus and an induced antiviral cell state. This competition is recapitulated in our model in terms of the two effective parameters P_a and R. Our model emphasizes the threshold dynamics, with a critical transition between effective confinement and unhindered spread that depends sensitively on the details of the relevant cell states. In particular, the presence of the specialized IFN-producing N cells allows for disease confinement at a much lower concentration of pre-antiviral cells (lower value of P_a) in the NOVAa model, than in the OVA model which lacks the N state (Fig. 4b). As a consequence of low P_a, the number of final N -state cells is also much lower.

The low concentration of ready-to-fight cells may seem perplexing, leading one to surmise that the organism could easily fight off an infection by only slightly increasing its investment in these primed cells. However, one should keep in mind that e.g. the human organism does indeed have ready-to-fight cells that are able to eliminate most foreign RNA, and only leave a few truly infectious viruses. As highlighted in the introduction, these select viruses often employ strategies to lower the P_a, for example by only being sensed by a small fraction of the RNA virus-sensitive receptors of our cells.

The parameter P_a can be interpreted as the probability that a cell is sufficiently antiviral to convert to a N state upon infection with a given virus. The relevant value of P_a will depend on the virus considered (and will be small for viruses that inhibit cell responses to infection) as well as on the host (e.g. on age (Kissler et al., 2020) and recent infection history). Dysregulated IFN responses are characteristic of the effective immunomodulatory strategies used by betacoronaviruses (Channappanavar et al., 2019; Acharya et al., 2020).

The parameter R reflects the signaling efficiency of an interferon-producing cell. Since R is measured in units of the typical distance that the virus spreads, it depends on viral properties, including its burst size, diffusion, and adsorption to host cells, with higher adsorption being associated with larger R values. For SARS-CoV-2 this suggests that lower ACE2 receptor counts would result in less adsorption to nearby cells, in turn allowing the virus to spread to more distant tissues (Bastolla, 2021) suggesting a lower value of R.

For viruses that do not delay the production of IFNs, P_a would be higher than for SARS-CoV-2, allowing neighbor cells around an infected site to form a kind of “ring vaccination” as the antiviral state dominates. In this sense, our model is consistent with the previous modeling of the roles of autocrine and paracrine interferon signaling suppression of viral infection (see e.g. (Michael Lavigne et al., 2021) for parallels between IFN response and ‘ring vaccination’).

We do not consider viral particles which enter the bloodstream and seed new infections non-locally. This may allow the virus to spread in the tissue at what would otherwise constitute sub-critical conditions in our model. Further, there may be tissue-specific variations in both P_a and R, adding larger-scale heterogeneity to the overall spreading. As the disease progresses one would expect additional heterogeneity to emerge, associated with variability in later host responses including macrophage activation and adaptive immunity (Wang et al., 2021).

The remarkable heterogeneity of disease progression in COVID-19, in the form of widely variable symptoms (Tabata et al., 2020) and transmission risk (Nielsen et al., 2021; Kirkegaard and Sneppen, 2021), has been widely observed. E.g. among university students, just 2% of SARS-CoV-2 positive hosts provided 90% of total respiratory viral load (Yang et al., 2021). In our formalism, we would understand such variability in terms of a P_a that is comparable to the critical value, but varying between hosts. A slight change of P_a then results in dramatic fluctuations in the outcome of an infection.

To be more quantitative, for SARS-CoV-2 the detected virus count on average grows by a factor of 3.5 (Kissler et al., 2020) in one infection generation of 8 hours (not to be confused with the between-host ‘generation time’ of the infection). This within-host reproductive number is far below the number of viruses produced from a cell, indicating severe restrictions from the innate immune system. On the other hand, 3.5 is still above the threshold for spreading, indicating that within-host amplification is super-critical. However, the measured amplification includes viruses that “jump” to other spots in an infected person, thereby suggesting a local spreading that is closer to the critical value than an amplification of 3.5 would suggest.

Our study finally compared the NOVAa model with the simpler OVA scenario that recapitulates earlier modeling of induced antiviral states (Howat et al., 2006; Segredo-Otero and Sanjuán, 2020; Michael Lavigne et al., 2021). These papers all build on a more homogeneous role of infected cells, each inducing some immunization of surrounding cells. Ref. (Howat et al., 2006) emphasizes the larger range of IFN signals compared to viral diffusion, while the focus of (Segredo-Otero and Sanjuán, 2020) is competition between viruses with different abilities to suppress IFN signaling. Ref. (Michael Lavigne et al., 2021) introduces a cellular automaton approach where the antiviral state leads to a type of ring vaccination that prevents the virus from spreading when more IFN is produced. Our OVA model may be seen as a simplified and more stochastic version of this last model. The NOVAa model then adds the additional benefits associated with the experimentally observed but low-abundance N state cells, which by their rarity adds to predicted randomness between the fate of individual infection centers during an early viral infection.

Acknowledgements

BFN acknowledges financial support from the Carlsberg Foundation in the form of an Internationalisation Fellowship (grant no. CF23-0173) and through the Carlsberg Foundation’s Semper Ardens programme (grant no. CF20-0046), as well as from the Danish National Research Foundation (grant no. DNRF170). Novo Nordisk Foundation Center for Stem Cell Medicine is supported by Novo Nordisk Foundation grant NNF21CC0073729.