key: cord-0821730-75odt5n8
authors: Sego, T. J.; Aponte-Serrano, Josua O.; Ferrari Gianlupi, Juliano; Heaps, Samuel R.; Breithaupt, Kira; Brusch, Lutz; Crawshaw, Jessica; Osborne, James M.; Quardokus, Ellen M.; Plemper, Richard K.; Glazier, James A.
title: A modular framework for multiscale, multicellular, spatiotemporal modeling of acute primary viral infection and immune response in epithelial tissues and its application to drug therapy timing and effectiveness
date: 2020-12-21
journal: PLoS Comput Biol
DOI: 10.1371/journal.pcbi.1008451
sha: 57f1aa2988a8f8e9ac0f130bfc74735030a51247
doc_id: 821730
cord_uid: 75odt5n8

Simulations of tissue-specific effects of primary acute viral infections like COVID-19 are essential for understanding disease outcomes and optimizing therapies. Such simulations need to support continuous updating in response to rapid advances in understanding of infection mechanisms, and parallel development of components by multiple groups. We present an open-source platform for multiscale spatiotemporal simulation of an epithelial tissue, viral infection, cellular immune response and tissue damage, specifically designed to be modular and extensible to support continuous updating and parallel development. The base simulation of a simplified patch of epithelial tissue and immune response exhibits distinct patterns of infection dynamics from widespread infection, to recurrence, to clearance. Slower viral internalization and faster immune-cell recruitment slow infection and promote containment. Because antiviral drugs can have side effects and show reduced clinical effectiveness when given later during infection, we studied the effects on progression of treatment potency and time-of-first treatment after infection. In simulations, even a low potency therapy with a drug which reduces the replication rate of viral RNA greatly decreases the total tissue damage and virus burden when given near the beginning of infection. Many combinations of dosage and treatment time lead to stochastic outcomes, with some simulation replicas showing clearance or control (treatment success), while others show rapid infection of all epithelial cells (treatment failure). Thus, while a high potency therapy usually is less effective when given later, treatments at late times are occasionally effective. We illustrate how to extend the platform to model specific virus types (e.g., hepatitis C) and add additional cellular mechanisms (tissue recovery and variable cell susceptibility to infection), using our software modules and publicly-available software repository.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

The current global pandemic of COVID-19, caused by the novel coronavirus Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), has motivated the study of beta coronavirus diseases at multiple spatial and temporal computational modeling scales [1] . The time course, severity of symptoms and complications from SARS-CoV-2 infection are highly variable from patient to patient [2] . Mathematical modeling methods integrate the available hostand pathogen-level data on disease dynamics that are required to understand the complex biology of infection and immune response to optimize therapeutic interventions [3] [4] [5] . Mathematical models and computer simulations built on spatial and ODE frameworks have been extensively used to study in-host progression of viral infection [6] , with a recent acceleration in the development of spatial COVID-19 viral infection models in response to the global pandemic [7, 8] .

Building multiscale models of acute primary viral infection requires integration of submodels of multiple biological components across scales (e.g., viral replication and internalization, immune system responses). Non-spatial, coupled ordinary differential equation (ODE) models can represent many aspects of pathogen-host interaction. In the context of viral infection dynamics, specialized ODE models can describe both the entire virus-host response at the tissue and organ levels and different stages of the viral replication cycle within cells, such as binding and internalization [9, 10] , viral genome replication and translation [11, 12] , assembly, packaging and release [13, 14] . By fitting ODE models to clinical or experimental data, researchers have been able to estimate important parameters, such as the turnover rate of target cells, average lifetimes of viral particles and infected cells and the rate of production of new viral particles by infected cells [15] . Other ODE models include pharmacokinetic models of drug availability [16] .

The simplest non-spatial models assume that the distribution of the modelled quantities (e. g., cells, viruses, chemical species) are uniformly distributed in space and time [17] . This assumption might not be realistic in solid tissues, where viruses and host immune cells are not usually distributed homogeneously and infection propagates locally [15] or in situations where transport limits the dynamics (e.g. the migration of antigen presenting cells to lymph nodes, the transmission of virus between organs or the microdosimetry of a drug therapy). By and attract immune cells within the tissue by chemotaxis [40] [41] [42] . The early innate immune response activates a number of cell types including dendritic cells, macrophages, neutrophils, mast cells, basophils, eosinophils, leukocytes, and natural killer (NK) cells [43] . Many of these immune cells themselves release both pro-and anti-inflammatory signals. Immune signals also recruit circulating neutrophils in the blood and, later, activate cytotoxic innate immune cells like NK cells within the tissue, which kill infected cells through release of diffusible factors and contact-mediated interactions, respectively.

The temporal dynamics of the concentration of extracellular virus varies greatly among virus families, tissues and host species [44] . However, for many viruses, including influenza and coronaviruses, once infected cells begin to release virus, the amount of extracellular virus increases exponentially over a period of a few days, reaching a peak during an early phase of infection [45] . As the viral load increases, immune signaling increases rapidly (this increase is associated with the onset of fever and other symptoms) recruiting more circulating cells of the innate immune system to the infection site [46] .

Immune signals from infected cells and innate immune cells help trigger the adaptive immune response. Macrophages and dendritic cells that have engulfed and degraded viral . Initial innate immune responses include phagocytosis of virus by neutrophils and macrophages, Type I interferon-induced antiviral resistance (IFN) (dark blue) and killing of infected cells by Natural Killer (NK) cells and other cell types (red). The black vertical dashed line denotes the transition between innate and adaptive immune responses. The adaptive immune response is triggered both by cytokine signaling to the lymph nodes and the migration of antigen-presenting cells from the tissue to the lymph nodes (not shown). In the later phases of infection, innate immune responses continue, but additional adaptive immune components come into play, including virus-specific cytotoxic T-cells (light blue) kill infected cells directly and also kill nearby cells through a variety of mechanisms. The orange vertical dashed line denotes the onset of the humoral adaptive immune response. B-cells produce virus-specific antibodies (orange line) which bind and inactivate virus directly and also allow its clearance and clearance of infected cells by other cell types. Tissue damage (shaded purple curve) accumulates due to cell death from direct responses to virus and from immune-cell killing by contact-mediated, diffusible factor-mediated and bystander-mediated mechanisms and eventually dissipates as cells proliferate to repair the damage (Adapted from [50, 51] ). The specific time course of all components varies among viruses, host tissues and host species, but the general sequence of events and immune response components are generally preserved.

https://doi.org/10.1371/journal.pcbi.1008451.g001 PLOS COMPUTATIONAL BIOLOGY pathogens or fragments of dead infected cells (i.e., phagocytosis) migrate over a period of days to nearby lymph nodes and serve as viral antigen presenting cells (APCs) to naive T-cells. Antigen presentation induces naive T-cell proliferation and differentiation into pathogen-specific memory and effector T-cells [47] . Effector T-cells migrate to the site of infection and induce apoptosis of infected cells by antigen recognition and contact killing. Both infected and uninfected cells in contact with dying cells can die through bystander-effect mechanisms. In acute infections, adaptive immune response leads to pathogen neutralization and clearance [48] . Viral loads usually decrease rapidly as adaptive immune cells like CD8+ T-cells enter the tissue and eliminate infected cells. Cells also begin to send out anti-inflammatory signals, shutting down the immune response as viral clearance proceeds. Antigen presentation within the immune system also induces activation of naive B-cell lymphocytes into antibody-producing memory B-cells and plasma cells, which leads to the production of antibodies. The adaptive immune response remembers its exposure to previous pathogens and provides the body with pathogen-specific effector cells and antibodies which neutralize and clear them, providing long term immunity [49] . Tissue damage results from virus and cytokine-induced cell death (which is first noticeable after 2 or 3 days) and from killing of infected and uninfected cells by immune cells, which increases steadily until the end of viral clearance (7-10 days). Tissue recovery and healing start around the time of viral clearance and may last for several weeks.

In this paper we develop a framework for the multiscale multicellular spatiotemporal simulation of the complex processes of infection and immune response in a small patch of epithelial tissue. The model provides a representation of the complex biology that reproduces key observed emergent behaviors of infection dynamics. We create representations of the main types of components and biological mechanisms associated with acute, primary viral infection and immune response, with a special emphasis on modularity of mathematical forms and computational implementation to support the development of more detailed models in future work (e.g., the creation of additional cell types, signals and detailed cell responses of various aspects of the immune response). We illustrate such an activity of development by integrating a detailed viral replication model for hepatitis C virus as an extension to the framework (see Integration of an explicit RNA synthesis model allows the spatiotemporal modelling of hepatitis C virus infection).

Our base model consists of three interconnected components (Fig 2A) : an epithelium component, an extracellular environment component and a lymph node component. The model represents the epithelium as a compact monolayer of initially identical immobile (which is appropriate for an epithelium) epithelial cells that it classifies by their current state of viral infection (i.e., uninfected, infected, virus releasing, dead, Fig 2C) . These epithelial cells are initially identical in their number of viral receptors (though we show how to include heterogeneity in Heterogeneous susceptibility inhibits spread of infection). The cells internalize extracellular virus, modulate their number of surface receptors, replicate virus, release virus and die in response to virus production, and release an extracellular cytokine signal when infected. Our model omits interferon-induced antiviral resistance, which can be implemented as a model extension using mechanisms demonstrated in this work (e.g., secretion of cytokine, modulation of internalization parameters, see S2 Text, Developing a Model Extension in CompuCell3D).

The model represents the extracellular environment as a space above the epithelium which provides the space in which immune cells are recruited and move, and into which cells release viruses and chemicals. We include a single type of immune cell that exhibits many key immune-cell behaviors associated with macrophage, neutrophil, NK cell and T-cell roles, including activation, chemotaxis, relaying and amplification of cytokine signals, contact killing, secretion of diffusible killing factors, and bystander killing, to represent the main types of Together the epithelial-cell, extracellular-environment and immune-cell components represent the epithelial tissue. Each model object has associated processes that dictate its states and behaviors. Epithelial-cell processes include viral internalization (E1), viral replication (E2), viral release (E3) and cell death (E4). Immune cell processes include activation (I1), chemotaxis (I2), contact cytotoxicity (I3) and oxidative cytotoxicity (I4). Activated immune cells participate in oxidative cytotoxicity (I4) and secrete oxidative agents into the oxidizing-agent field (T3). Activated cells become inactive after 1 hour. The virus field (T1), cytokine field (T2) and oxidizing-agent field (T3) describe spatially-varying densities of extracellular components. Field processes describe diffusive transport and clearance of material in the extracellular environment and activated transport to the lymph nodes. The lymph node is a single-compartment model whose pro-or anti-inflammatory state specifies the recruitment or removal rate (L1) of immune cells in the epithelial tissue. The transport of cytokines to the lymph node promotes its proinflammatory state. (B) Viral Life Cycle: Interactions in the viral internalisation, replication and release models. Schematic representation of inputs, outputs and interactions between stages of the viral replication model. Extracellular viral particles (represented as continuous fields) are internalized by the viral internalization model and initiate the viral replication model. The main stages of the viral replication model are: unpacking, viral genome replication, protein synthesis and viral assembly and packaging. The output of the viral replication model passes to the viral release model, which transfers newly assembled viral particles from the cells into the extracellular environment. (C) Cell types and transitions. Epithelial cells are of type uninfected if they have not yet internalized any virus (E1). They are of type infected if they have internalized virus, but are not releasing virus into the virus field (viral release E3 is inactive). They are of type virus releasing if they are releasing virus into the extracellular virus field (i.e., viral release E3 is activated). Immune cells are initially inactive and do not participate in the oxidative cytotoxicity (I4) or chemotax towards higher concentrations of the cytokine field (I2). They become activated when they experience activation (I1). In all panels, dashed arrows with barbed heads represent transformations, solid arrows with barbed heads represent transport and solid arrows with lollipop heads represent regulation.

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

immune-cell mediated cell killing, rather than any particular immune cell phenotype. We omit macrophages, neutrophils and their phagocytosis and signaling, which can be represented using simple extensions of the immune-cell type. We do not model contact interactions between immune cells (e.g., by T-cells and APCs). We also do not model tissue recovery, though we demonstrate examples of adding tissue recovery models in Model extensions.

We simulate extracellular-virus particle density as a continuum field and particle transport and clearance as continuous diffusion and decay. We approximate the discrete process of a cell's internalization of a virus particle by a stochastic virus internalization event (E1) determined by time elapsed, the local concentration of the virus field, and the number of available cell-surface receptors on the cell. We simplify the complexity of viral replication into four steps: unpacking, viral genome replication, protein synthesis and packaging/assembly (E2, Fig  2B) [7, [52] [53] [54] . The subcellular kinetics of viral replication determine the rate of release of new virus particles into the extracellular environment (E3). To represent the combined effect of the many types of virus-induced cell death, each infected epithelial cell has a probability of dying that depends on the number of assembled viral particles inside the cell per unit time (E4).

We simplify the complex biochemistry of the many molecular signals active in epithelial tissues, which include chemokines, interferons and viral fragments, into a single generic extracellular cytokine field that acts both as a tissue-local and systemic pro-inflammatory signal. Infected epithelial cells and immune cells both secrete cytokine (T2). The cytokine field produces local immune effects such as activation (I1) and chemotaxis (I2) of immune cells. Activated immune cells can revert back to inactive immune cells when the cytokine signal ceases. The cytokine field also recruits immune cells to the tissue through long-distance signaling via the lymphatic system (L1). We model recruitment of immune cells to the simulation domain using an ordinary differential equation for the immune signal (S), which represents the balance between pro-and anti-inflammatory signaling and the delay due to antigen-presenting cell transport from the tissue through the lymphatic system to the lymph node and due to the time required for T-cell amplification. In the absence of infection, the lymph node maintains a small resident immune cell population in the tissue. Immune cells can cause epithelial cell death (E4) by three mechanisms: 1) contact cytotoxicity; 2) the bystander effect; and 3) through the release of an oxidative agent. Immune cells kill infected epithelial cells by contact cytotoxicity, in which case neighboring uninfected, infected and virus-releasing epithelial cells can also die through a bystander effect (I3). In regions of the tissue with high cytokine levels, immune cells secrete a diffusive oxidative agent to the environment (T3) that kills uninfected, infected and virus-releasing epithelial cells (I4).

We begin by presenting our base multicellular model of viral infection in an epithelial tissue, along with a simulation for a baseline set of parameters and basic analyses. We then explore the simulation's dependence on critical parameters that are important to the dynamics of acute primary viral infection in airway epithelial cells. All simulations and spatial, population and system-level metrics presented in this section follow the specifications in Simulation Specifications. We performed simulations using the open-source modeling environment Compu-Cell3D [55] . Instructions on how to run these simulations are provided in our publicly available repository at https://github.com/covid-tissue-models/covid-tissue-response-models/ tree/master/CC3D/Models/BiocIU/SARSCoV2MultiscaleVTM.

We initialize the simulations with periodic boundary conditions parallel to the plane of the sheet and Neumann boundary conditions perpendicular to the plane of the sheet. Initially, the extracellular environment does not contain any extracellular virus, cytokines, oxidative agents or immune cells. We introduce infection by creating a single infected epithelial cell at the center of the epithelial sheet, comparably to (but less than) initial conditions of similar works that model discrete cellular distributions [18, 31] . To illustrate the full range of dynamics of viral infection in the presence of an immune response, we established a baseline set of parameters (Table 1) for which the immune response is strong enough to slow the spread of the infection, but insufficient to prevent widespread infection and death of all epithelial cells (Fig 3) . While we have adjusted the parameters for the viral replication model to agree with reported time scales for SARS-CoV-2 replication in vitro [56] , and we have selected parameter values in physiologically reasonable ranges, we have not attempted to match other model parameters to a specific tissue, virus or host species. Furthermore, this baseline parameter set is not unique with respect to its purpose, in that many sets of parameters can generate appreciable but insufficient inhibition of spread of infection (see S12-S17 Figs). Rather, as is shown in subsequent sections, this parameter set presents emergent dynamics of a theoretical virus and host immune response near, but not in, a regime of successful prevention of widespread infection, which is critical to showing the effects of underlying mechanisms on emergent dynamics and resulting outcomes.

The infected cell immediately starts releasing cytokines into the extracellular environment. After an incubation period (~150 minutes, 2 ½ hours), the first infected epithelial cell (green) transitions from infected to virus releasing (red), and starts releasing viruses into the extracellular environment. Initial release of extracellular virus causes additional epithelial cells to become infected. Release of cytokines leads to delayed addition of immune cells to the simulation domain (Fig 3D) . By 4000 minutes (67 hours, 2 ¾ days), the number of infected cells increases 10-fold and the epithelial cells start undergoing virally-induced death as the infection spreads radially outward from the initial site. The increase in the number of infected cells and the local cytokine concentration is accompanied by a similar increase in the local immune cell population. By 8000 minutes (133 hours, 5 ½ days), the number of dead epithelial cells around the initial infection site increases sharply. Following this phase of rapid cell death, the number of infected, virus-releasing and dead epithelial cells continues to increase exponentially but at a slower rate. This transition results in the formation of an annular region of infected cells spreading radially outwards from the initial infection site (Fig 3A) , analogous to the Fisher equation for deterministic front propagation [58] . Total extracellular virus and cytokine continue to increase exponentially. The increase in cytokine results in continued recruitment of additional immune cells. By 16000 minutes (267 hours, 11 days), the number of uninfected epithelial cells reaches zero and the number of infected and virus-releasing cells peaks. Despite the declining number of infected and virus-releasing epithelial cells, the delayed immune response continues to add immune cells to the tissue. After about 16000 minutes (267 hours, 11 days), the extracellular virus and the amount of cytokine decrease exponentially as the remaining virus-releasing epithelial cells die. By 20000 minutes (333 hours, 14 days), all epithelial cells die and many immune cells leave the tissue.

The rate at which infection propagates and the strength (speed and amplitude) of the immune response depend on multiple model parameters. Interplay between these rates leads to a wide spectrum of qualitatively-distinct spatiotemporal dynamics. The virus-receptor binding affinity k on and the immune response delay coefficient β delay are critical parameters affecting the rate of infection of epithelial cells and the strength of the immune response, respectively. Larger k on corresponds to a higher rate of infection propagation (increasing k on increases the rate of internalization of extracellular viral particles into epithelial cells, see Eq (3) in Models Contact coefficients J (all interfaces) 10 Selected comparably to [68] for low adhesion immune cell-immune cell and immune cell-medium interfaces 1 The diffusivity in water for a virus of radius 0.1 microns like SARS-CoV-2 according to Stokes-Einstein is about 3 microns 2 /s. The average steady-shear viscosity for lung mucus varies significantly and is shear thinning, but in the more viscous regions is found to vary for frequencies between 10 −4 and 102 Hz, spanning viscosity values as high as 103 Pa-s and as low as 10 −2 Pa-s. In general, at low shear rates, the viscosity of human mucus is as high as 104−106 times that of water [75] . Thus the minimal diffusion constant possible would be 3 x 10 −6 microns 2 /s and the maximal rate in water would be 3 microns 2 /s. 0.01 microns 2 /s is a reasonable geometric interpolation.

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

and methods). Larger β delay corresponds to weaker immune response (decreasing β delay increases the strength of immune-cell recruitment, see Eqs (12)- (14) in Models and methods).

Varying these two parameters around the baseline simulation values yields six patterns of spatiotemporal infection dynamics, ranging from unopposed infection to clearance (Fig 4) . We defined these classes based on the transient dynamics and the final state of the simulation at 20000 minutes (333 hours, 14 days). We terminated the simulations at 20000 minutes (333 hours, 14 days) because we assume that, in real tissues, additional adaptive immune responses at this time generally lead to rapid viral clearance. As a result, any epithelial cells uninfected at the end of the simulation would likely remain uninfected if we included these additional immune mechanisms. We define the six patterns (classes) of infection dynamics as follows:

No immune response: a limiting case (corresponding to in vitro and organoid culture experiments on viral infection, which lack immune cells) that serves as a reference simulation showing the spread of unopposed infection. When the cellular immune response is absent, an infection front travels across the epithelium until all epithelial cells have died due to intracellular virus (Fig 4A and S1 Fig) .

Widespread infection: when the immune response is weak (large β delay ) or the rate of infection propagation is large (large k on ), the immune cells kill enough infected epithelial cells to reduce both the concentration of extracellular virus and the propagation of the infection front. However, at the end of the simulation no uninfected epithelial cells survive (Figs 3 and 4B).

Slowed infection: for moderate immune response (moderate β delay ) and a moderate rate of infection propagation (moderate k on ), both uninfected and infected epithelial cells and some extracellular virus remain at the end of the simulation (Fig 4C) . These cases are functionally distinct from widespread infection, since even a single remaining uninfected epithelial cell could initiate tissue regeneration. In most cases of slowed infection, the numbers of infected cells and the extracellular virus continue to increase. A special case of slowed infection occurs when oscillations in the amount of virus (S2A Fig). Containment: for strong immune response (small β delay ) and low to moderate rate of infection propagation (moderate k on ), a few infected and virus-releasing cells are present for most of the simulation. However, the immune cells eventually kill all infected and virus-releasing epithelial cells. At the end of the simulation, no infected or virus releasing cells remain, while uninfected cells survive and some extracellular virus remains in the extracellular environment ( Fig 4D) .

Recurrence: for strong immune response (small β delay ) and a fast infection propagation (large k on ), relatively few epithelial cells become infected early in the simulation. All infected and virus-releasing epithelial cells die. However, the remaining extracellular virus infects additional epithelial cells later on, restarting the spread of infection ( Fig 4E) . Clearance: for strong immune response (small β delay ) and a slow infection propagation (small k on ), the number of infected and virus-releasing epithelial cells goes to zero without recurrence and the extracellular virus drops below a threshold (of 1/900 per cell in our analysis), rendering the frequency of recurrence negligible (Fig 4F) . A special case of clearance Table 1 and Fig 3 except for k on and β delay (S1 Table) . To quantitatively characterize simulation results, we measured the number of uninfected, infected, virus-releasing and dead epithelial cells, the total number of immune cells, the number of activated immune cells, the total amount of extracellular virus (integral over the virus field), the total diffusive cytokine (integral over cytokine field), the maximum total extracellular virus (over all simulation time) and the maximum total diffusive cytokine (over all simulation time). 

To explore the effects of the rate of virus internalization and the strength of the immune response, we performed a multidimensional parameter sweep of the virus-receptor association affinity k on and immune response delay coefficient β delay . Variations in virus receptor association affinity represent factors that affect the binding affinity of cellular viral receptors (e.g., ACE2 and TMPRSS-2 in the case of SARS-CoV-2) with a virus (e.g., mutations in viral coat protein or drugs to block viral entry). Variations in immune response delay coefficient represent factors that affect the strength of the systemic immune response (e.g., anti-inflammatory corticosteroids, IL-7 treatment or age, since older individuals tend to have slower adaptive immune responses) [36] .

We ran ten simulation replicas for each parameter set, increasing and decreasing the baseline parameter values 10-fold and 100-fold (Figs 5-7). For each simulation replica, we examined the number of uninfected epithelial cells (Fig 5) , the number of infected epithelial cells (Fig 6) , the total extracellular virus (Fig 7) , the number of dead epithelial cells (S3 Fig Table 1 ). The number of uninfected epithelial cells for each simulation replica for each parameter set, plotted on a logarithmic scale, vs time displayed in minutes. See 

Optimal therapeutic use of antiviral drugs requires considering the relationship between molecular efficacy (how effectively the drug blocks a particular aspect of the viral life cycle at saturation concentration), potency of therapy (the effect of the drug at a molecular level at a 

given dose) and clinical effectiveness (how well the drug reduces the severity or duration of the infection), as well as the tradeoff between side effects and bioavailability. One drug might have moderate efficacy but have few side effects. Another drug might have high efficacy, but have severe side effects at high doses that limit the maximum tolerated dose or use of even moderate doses in prophylaxis. A drug might also have a combination of beneficial and adverse effects (e.g., it might reduce viral replication early in infection, but also be immunosuppressive) [13, 37] . Antiviral drugs like Tamiflu retain their ability to block aspects of the viral life cycle (efficacy), but become much less clinically effective as the time before treatment increases: (in adults Tamiflu is most effective when given within 48 hours after exposure and thus is often used prophylactically) [38] .

In this section we use our model to show the trade-offs between time-of-use and potency of a drug that targets viral genome replication in a host cell. Several antiviral medications for RNA viruses reduce the net viral replication rate by inhibiting synthesis of viral RNA by the viral RNA polymerase. We focus on RNA-synthesis blockers in this paper because viral genome synthesis exponentially increases the production rate of viruses per cell, while the other stages of viral replication have linear amplification effects (see Eqs (6)- (9) in Models and methods).

To simulate the effects of treatment that targets RNA synthesis using different drug efficacies and times of administration, we generated a series of simulations in which we reduced r max , the replication rate of genomic material in the viral replication model (Eq (7) in Models and methods), by different amounts and at different times in the simulation. The "viral replication multiplier" represents the potency of the treatment, the factor by which r max is reduced 

(either a low dose with high efficacy, or a high dose with a less efficacy). The "time delay of application" is the simulation time at which r max is reduced, which corresponds to the time after infection at which the treatment is administered. To characterize therapeutic effectiveness, we distinguished three classes of simulation outcomes:

Positive outcomes: effective treatment, where at least 50% of the epithelial cells remain uninfected at the end of the simulation (green-shaded subplots).

Negative outcomes: ineffective treatment, where less than 10% of the epithelial cells remain uninfected at the end of the simulation (orange-shaded subplots).

Intermediate outcomes: partially effective treatment, where between 10-50% of the epithelial cells remain uninfected at the end of the simulation (unshaded or intermediate-shaded subplots).

To characterize how the potency and time of initiation of treatment affect the dynamics of the simulation and treatment effectiveness, we examined the time courses of the number of uninfected epithelial cells (Fig 8) , virus-releasing epithelial cells (Fig 9) , the total amount of extracellular virus (Fig 10) When the treatment is given early, while the level of extracellular virus is increasing rapidly and exponentially (before 6000 minutes, 100 hours, 4 days) after infection, most of the simulation replica outcomes are positive, showing effective treatment (Figs 9-11, green-shaded subplots). If the drug is administered prophylactically or very soon after infection (at 0 minutes) the treatment potency needs to be only 25% to achieve mostly positive outcomes (effective treatment). Increasing the time to treatment increases the potency required to achieve similar numbers of positive outcomes: the treatment is effective for a potency of at least 37.5% if administered by 4000 minutes (67 hours, 2 ¾ days), and at least 87.5% if administered by 6000 minutes (100 hours, 4 days). For all potencies greater than 12.5%, early intervention prevents significant increase in the number of virus-releasing cells (Fig 9, green-shaded subplots) , and a small number of immune cells suffices to stop the spread of infection (S9 Fig, green -shaded subplots). In this region, delaying treatment results both in a higher level of extracellular virus (Fig 10, green-shaded subplots) and more dead epithelial cells at the end of simulation (Fig 8, green-shaded subplots). With inhibited viral replication in the infected epithelial cells, the extracellular virus decays until it is mostly cleared by the end of simulation (Fig 10) . Variability between simulation replicas for a given parameter set increases with both decreasing potency and increasing time of initiation of treatment.

If the potency of the treatment is 12.5% (or less), most of the simulation replicas have negative outcomes (low effectiveness), even if the drug is administered prophylactically or soon after infection (at 0 minutes) (Figs 9-11, bottom row). In these cases, the time after infection at which the drug is given makes no significant difference to the treatment effectiveness. When the treatment is given late (time delay of application of at least 10000 minutes, 167 hours, 7 days), regardless of the potency of the drug, most simulation replicas have negative outcomes (Figs 9-11, orange-shaded regions). By the time of treatment, a significant number of epithelial cells have been infected (more than 10% in most cases- 

even when a significant number of virus-releasing epithelial cells remain, indicating that viral replication inside cells has been significantly reduced. Later intervention also increases variability between simulation replicas and, although most simulation replicas have negative outcomes, the same set of parameter values produced two distinct qualitative outcomes (some more and some less favorable) for higher potency (S11 Fig, orange-shaded regions) . Thus in a few cases, even late treatment can be effective.

When the treatment is given at intermediate times (times between 6000 and 10000 minutes, 100 to 167 hours, 4 to 7 days), most simulation replicas have intermediate outcomes. For potencies above 50%, the fraction of uninfected epithelial cells at the end of simulation is relatively high (around 50%) and the treatment is usually moderately effective (Fig 8) . For potencies below 50%, the number of virus-releasing epithelial cells remains approximately constant or continues to increase after treatment (Fig 9) and significant levels of extracellular virus remain at the end of the simulation, and so in most cases the treatment is ineffective (Fig 10) . In this regime, variability between outcomes for the same parameter values is higher than for potencies above 50%. 

A particular parameter set (time delay of application of 10000 minutes, 100% potency) produced simulation replicas that had instances of all three outcomes (Fig 11) . In a simulation replica with a positive outcome (Run 2, Fig 11A) , the first uninfected epithelial cell dies (as well as a few uninfected epithelial cells) before 4000 minutes (67 hours, 2 ¾ days), after which the total extracellular virus gradually decreases. At around 4000 minutes (67 hours, 2 ¾ days), an epithelial cell near the initially-infected cell becomes infected, causing a recurrence of infection, whose spread was stopped by the treatment. In contrast, simulation replicas with intermediate and negative outcomes (Runs 8 and 4, respectively, Fig 11A) have comparable, and significantly more, numbers of infected and virus-releasing epithelial cells at 4000 minutes (67 hours, 2 ¾ days), while total extracellular virus is greater in the replicas with a negative outcome than in the replicas with an intermediate outcome. 

In this section we demonstrate the deployment of extensions to the framework described in Models and methods, which we will refer to as the "main framework" when discussing extensions in this and subsequent sections, as well as particularization of the framework to specific biological problems like a different virus. We accomplish this through integration of an existing model in the literature of hepatitis C virus (HCV), and elementary examples of adding models of heterogeneous epithelia and tissue recovery, one model of which is constructed from another. We include schematic representations where appropriate of the software implementation according to the architecture we have developed (and will continue to develop) to support broad scientific use through rapid, parallel model development, flexible model integration, and model sharing through our publicly available online repository (see S2 Text for an overview of basic deployment, implementation, and public distribution of extensions and S1-S3 Code Snippets for specific examples).

Heterogeneous susceptibility inhibits spread of infection. To demonstrate basic extensibility of the framework to model additional complexity associated with viral infection, we introduced a basic notion of heterogeneity to the epithelium of simulated scenario with the premise that not all epithelial cells can be infected (i.e., some cells are unsusceptible, as in an actual heterogeneous respiratory epithelium) [22] . We implemented heterogeneous susceptibility by randomly selecting a fraction of the epithelial cell population at the beginning of simulation and setting the number of surface receptors of each to zero (i.e., these cells have no surface receptors for internalization of the virus). Otherwise, all mechanisms and model parameter values used in simulations of heterogeneous susceptibility were the same as those used to generate results shown in Fig 3 ( i.e., all mechanisms described in Models and methods using the baseline parameter set). Implementation of randomly selected heterogeneous susceptibility is available in the add-on modules library and is hosted on our repository for public use with the module name "RandomSusceptibility". For each set of replicas we simulated ten replicas using the RandomSusceptibility module while varying the fraction of unsusceptible cells (like different locations in the lungs), from 10% to 50% unsusceptible in intervals of 10% ( Fig  12) .

For all fractions of unsusceptible cells, infection spread throughout the epithelial sheet ( Fig  12A) . No replica produced the exact definition of widespread infection, since all replicas had at least some remaining uninfected cells at the end of simulation. However, only eight out of fifty total replicas had a few epithelial cells that were not either uninfected or dead at the end of simulation (for infected cells, 1 replica for 10% unsusceptible, and for virus-releasing 1 replica for 10% and 20% unsusceptible, and 2 replicas for 30%, 40% and 50% unsusceptible). The distributions of infected and virus-releasing cells were notably different from all previous simulations with appreciable infection in that both subpopulations were noticeably intermixed with uninfected (presumably unsusceptible) cells and increasingly so with increasing fraction of unsusceptible cells. Most of these intermixed, uninfected cells (particularly those nearer to the initial site of infection) also died due to immune response mechanisms (i.e., bystander effect and oxidative killing).

The final number of uninfected cells at the end of simulation increased with an increasing fraction of unsusceptible cells (Fig 12B) . The rate of spread of infection decreased with an increasing fraction of unsusceptible cells (as observed by inspection in a rightward shift in the total number of uninfected cells). One replica for 20% unsusceptible cells exhibited significantly delayed spread of infection, though this was not observed in other replicas for any other fraction of unsusceptible cells.

Integration of an explicit RNA synthesis model allows the spatiotemporal modelling of hepatitis C virus infection. The viral replication model described in Models and methods describes the viral life cycle, from internalization to release, as occurring over four stages and PLOS COMPUTATIONAL BIOLOGY in unitless quantities. As such, it is possible to integrate detailed models of various stages of the viral life cycle into the main framework through appropriate substitution of viral replication model terms and assignment of unit quantities as necessary (i.e., one unit in the original viral replication model corresponds to a physical unit in an integrated model that explicitly describes a lifecycle process). Integration of a detailed model of some viral process then particularizes the virus represented in the main framework to the level of biological information introduced by the integrated model.

As a demonstration of model integration, we particularized viral replication of the framework to hepatitis C virus (HCV) using an existing model of subgenomic HCV RNA synthesis in Huh-7 cells [59] . The model of HCV RNA synthesis describes various aspects of subgenomic HCV replication like translation of HCV polyprotein in the cytoplasm, replication kinetics in vesicular-membrane structures and availability of host ribosomes, and explicitly models quantities of HCV RNA molecules, translation complexes, HCV polyprotein molecules, and necessary enzymes in both the cytoplasm and vesicular-membrane structures.

Integration of the HCV replication model casts the viral genome taking part in genome replication R of the viral replication module described in Models and methods as the cytoplasmic plus-strand HCV RNA molecules of the HCV replication model (denoted R cyt P , see S1 Text). We assumed that internalized virus from the viral internalization module converts into cytoplasmic RNA, and that some of the decay of cytoplasmic RNA molecules described in the HCV model leads to the production of quantities produced by replicating viral genome from the viral replication module. To relate the unitless viral replication model described in E2-Viral Replication to the biological quantities of the HCV replication model, we assumed that one unit of replicating viral genome corresponds to 100 HCV RNA molecules, which was found to produce total infection dynamics comparably to the main framework (for comparison with results of this work, rather than for reproduction in silico of any specific HCV data) when using both the baseline parameter set demonstrated in Fig 3 and all model parameters reported in [59] when the virus-receptor association affinity coefficient k on was increased by a factor of 100 (S2 Table) . The integrated HCV model is implemented in the add-on modules library using the aforementioned software architecture and is hosted on our repository for public use with the module name "HCVIntegrated" (Fig 13A) .

Ten simulations were executed using the integrated HCV model for a simulation time of two weeks using initial conditions similar to some used in [59] , where the initially infected cell was seeded with 500 cytoplasmic viral RNA molecules (i.e., initial R = 5, rather than initial U =1 as in all other simulations). Spread of infection with the integrated HCV model produced comparable patterns of spread of infection, where an infection front radially advanced outward from the initial site of infection ( Fig 13B) . However, one notable difference was that the prominent band of virus-releasing cells during spread of infection only occurred in some simulation replicas (Fig 13B, replica "a") , while in other replicas only small, isolated groups of cells became virus releasing (Fig 13B, replicas " b" and "c"). Variability of outcomes was particularly notable among the ten simulation replicas. Some simulation replicas produced widespread infection at a comparable timescale to that of the baseline parameter set, between one and two weeks ( Fig  13C) . Such simulation replicas were those that produced comparable spatial distributions of infected and virus-releasing cells (specifically, with a prominent band of virus-releasing cells) to those of the main framework using the baseline parameter set. In other simulation replicas (e.g., Fig 13B, replica "c"), slowed infection occurred due to early elimination of many virusreleasing cells. In such cases, the epithelial sheet had very few virus-releasing cells and scattered infected cells around the region of dead cells. 

Many infected cells remained infected over a period of over a week, which was not observed using the viral replication model of the main framework. Such infected cells did not contribute to spread of infection, but rather diminished the likelihood of widespread infection by increasing the distance between uninfected and virus-releasing cells. By the end of simulation in replicas that produced slowed infection, total extracellular virus was negligible despite a significant number of infected cells, which presents an outcome not described in Classification of spatiotemporal infection dynamics that could be called "benign infection".

An extensible framework architecture enables the inclusion of tissue recovery. As a demonstration of modularity and extensibility, we developed, implemented and tested two models of tissue recovery, where dead cells are replaced over time with uninfected cells. Since epithelial cells in the simulated epithelial sheet are static, removal of dead cells and proliferation of uninfected cells was modeled as the changing of the type of a dead cell to uninfected (rather than explicitly modeling mitosis, [18] ). To generate a scenario of viral clearance and significant tissue damage, we simulated the baseline parameter set with the parameter variations in the top-right corner of Figs 5-7, where viral internalization is severely inhibited and the immune response is very strong and fast (i.e., parameters k on and β delay were reduced by a factor of 100 so that the virus is cleared but many uninfected cells die).

In the first model of tissue recovery, called "Simple Recovery", dead cells are replaced by an assumed layer of proliferative cells underneath the simulated epithelial patch, and so each dead cell has a fixed probability of recovering. We approximated the probability of recovery based on an assumed onset of tissue recovery of 7 days, in which case the probability of recovery for each dead cell over a 20-minute simulation step was 1.98×10 −3 . In the second model of tissue recovery, called "Neighbor Recovery", dead cells are replaced by nearby uninfected cells similarly to wound healing, and so each dead cell has a probability of recovery equal to the number of neighboring uninfected cells multiplied by a coefficient. For comparison of results to those from the Simple Recovery add-on module, we used the same probability coefficient value such that the probability of recovery according to Simple Recovery for each dead cell over a 20-minute simulation step was 1.98×10 −3 per unit of contact area with neighboring uninfected cell (measured in number of neighboring lattice sites).

Both models were implemented in CompuCell3D in the add-on modules library using the aforementioned software architecture and is hosted on our repository for public use. The Simple Recovery model is hosted with the module name "RecoverySimple", and the Neighbor Recovery model is hosted with the module name "RecoveryNeighbor". Since the only difference between the two recovery models is the criterion for cell recovery (i.e., whether a fixed probability, or a probability according to the neighborhood of a cell), the implementation of the Neighbor Recovery model inherits all features of the Simple Recovery model implementation using basic Python class inheritance, and required overwriting only one function that implements a criterion of cell recovery during development (Fig 14A, see S2 Text, Extending a model in CompuCell3D).

All simulation replicas for each recovery model began with tissue insult due to oxidative killing by a strong and fast immune response through the first few days (Fig 14B) . All simulation replicas experienced a loss of approximately 100-200 uninfected cells (Fig 14C) . 

Our spatial, multicellular model of primary acute viral infection of an epithelial tissue includes key aspects of viral infection, viral replication and immune response. By investigating sensitivity to model parameters and simulating drug therapies, we identified six distinct spatiotemporal classes of infection dynamics based on the model's transient behaviors and final configurations. Each of our simulation-defined classes corresponds to biologically or clinically observable factors and outcomes. The case of no immune response would be useful for analyzing in vitro experiments (e.g., organoids). Widespread infection corresponds to an initial infection that is likely to spread to surrounding tissue and cause major tissue damage. Slowed infection corresponds to an initial infection whose spread is more likely to be eliminated by the adaptive immune response. Containment corresponds to immune-cell elimination of all infected cells but where remaining extracellular virus could result in new sites of infection elsewhere. Recurrence corresponds to the situation when new lesions form within the observed tissue patch. Clearance corresponds to immune-cell-based elimination of all infected cells and extracellular virus (classical viral clearance).

We showed that key parameters of the model, such as those affecting viral internalization (i. e., virus-receptor association affinity k on ), can lead to both containment/clearance (e.g., small k on , Figs 5-7) and widespread infection (e.g., large k on , Figs 5-7). Multidimensional parameter sweeps showed how the interplay between immune response (e.g. immune response delay coefficient β delay ) and viral spread could lead to widespread infection (e.g., large β delay , large k on , Figs 5-7), rapidly cleared infection (e.g., small β delay , small k on , Figs 5-7) or containment/ clearance after substantial damage (e.g., small β delay , moderate k on , Figs 5-7). Some of these outcomes would be expected biologically (e.g., very fast internalization with a slow immune response is likely to lead to widespread infection; faster and stronger immune responses should control the spread of viral infection within the tissue [ ) and would also occur in deterministic non-spatial models. Others, like the coexistence of replicas with containment/clearance or failure to control for the same parameter set, are less expected, and could not occur in a deterministic non-spatial model (though they might occur in some stochastic non-spatial models). We have observed this interplay of parameters, as well as the potential for stochastic outcomes, in variations of other parameters of the model, whether related to spatial (e.g., viral and cytokine diffusion coefficients) or deterministic and stochastic cellular aspects (oxidative agent threshold for death and virally-induced apoptosis dissociation coefficient, see S12-S17 Figs).

We studied the influence of timing and potency of an RNA-polymerization inhibitor like remdesivir [60] on the spread of viral infection within tissue (Figs 9-11 ). As expected, in our model, drugs with this mode of action can improve viral control in tissue if administered prophylactically at high potency, and their effectiveness decreases the later they are administered. Less obviously, the lower-left region of Figs 9-11 shows how therapies with even reduced potency could control the infection when administered sufficiently early, consistently with predictions from a deterministic, non-spatial ODE model (though the mode of action is not explicitly described) [61] . While we expect prophylactic or early treatment at the same potency to be more effective than later treatment, our model suggests that, for antivirals, time of treatment is a more significant factor than potency in determining the effectiveness of the therapy. Our model thus suggests that drugs that interfere with virus replication are significantly more effective if used even at very low doses prophylactically or very soon after infection, than they would be if used even at a high dose as a treatment given later after exposure. Specifically, a prophylactic treatment in simulation which reduces the rate of viral RNA synthesis by only 35% (35% potency) is more effective than a treatment with 100% potency given two and a half days after infection, and has about the same efficacy as a treatment with 50% potency given one day after infection. Our model also showed that because of stochasticity in viral spread, later treatment at moderate to high potency may still be effective in a subset of individuals.

Both parameter sweeps had regions with little variation in outcome between replicas (e.g., the upper-right and lower-left corners of Fig 5) . In regions of the parameter space between these extremes (e.g., the unshaded areas in Figs 5 and 9 ), different replicas showed dramatically different outcomes. One such parameter set in our drug therapy simulations produced three distinct qualitative outcomes (i.e., positive, intermediate and negative outcomes, Fig 11) . For these parameters, replica outcomes were particularly sensitive to stochasticity early in infection when only a few cells were infected (Fig 11A) , with delayed spread of infection from the first infected cell producing more positive outcomes. Simulation replicas with negative outcomes (Fig 11A, Run 8 ) had higher extracellular virus levels at earlier times than those with intermediate outcomes (Fig 11A, Run 4) , even though the fraction of each cell type was similar. Since the viral replication module is deterministic, the primary cause of this difference is the spatial distribution of cells. Spatial structure (e.g., infection of neighboring cells), stochastic events (e. g., early cell death of infected cells before significant virus release) and cell-to-cell variation (e. g., difference in viral release between cells) all affect the variation between replicas.

Differences in spatiotemporal dynamics and variability of outcomes thus critically depend on the ability of the model to resolve the spread of virus and immune response spatially. The intrinsic stochasticity of many model processes makes the spatial patterns of the infection front and distribution of tissue damage nontrivial. The spectrum of outcomes in our parameter sweeps (Figs 5-7 and 9-11) depends not only on parameter values and model immune response, but also on the emergent spatial patterns of cytokine and virus fields (e.g., variations within the infection front expose different numbers of uninfected epithelial cells to the immune response). Such stochastic and spatial aspects can also introduce new considerations to ODE models that have been primarily employed in a non-spatial context. For example, as described in the original presentation of the HCV model that was integrated in Integration of an explicit RNA synthesis model allows the spatiotemporal modelling of hepatitis C virus infection, the subgenomic kinetics of the HCV model require a minimum number (seven) of cytoplasmic viral RNA molecules to reach a saturated state. When employing the HCV model in a multicellular, heterogeneous context, insufficient internalization of a spatially heterogeneous extracellular viral field for subgenomic replication to produce rampant viral production, leading to insufficient cytokine signaling to invoke further immune response, makes possible the so-called outcome of benign infection.

Our modeling framework can improve with the inclusion of additional cellular and immune mechanisms discussed in Fig 1. The modularity of model modules and built-in extensibility of the publicly available software implementation enables us, and other interested members of the scientific community, to accomplish such activities collaboratively or independently, concurrently, and even when in theoretical disagreement (see S2 Text). Modules accounting for viral clearance, tissue recovery and persistent adaptive immune response can be added to model later stages of disease progression (as demonstrated in Model extensions). The current immune model does not include important signaling factors (e.g., interferon-induced viral resistance in epithelial cells) and the different roles of tissue-local and systemic signals (e.g., various cytokines). It also omits many cell types associated with both innate and adaptive immune response and their roles (e.g., viral scavenging by macrophages, relaying and amplification of immune signals by dendritic cells). Of special interest to results like those presented in this work is the effect of specific roles by individual immune cell phenotypes on emergent dynamics and outcomes, considering that the timing of their activities during progression of events can be quite different (e.g., early neutrophil release of oxidative agent contrasted with later effector T-cell contact-mediated killing). Such details, which we are currently pursuing, are particularly important for using framework to interrogate the spatiotemporal details of the immune response, which are poorly understood. The model does not currently consider the production and role of antibodies in the humoral immune response or tissue recovery after damage (like the demonstration models presented in An extensible framework architecture enables the inclusion of tissue recovery). The model also greatly simplifies the structure of the epithelium and its environment, but could be easily generalized to a detailed, three-dimensional geometry, albeit at the cost of computational performance.

Our current results suggest priorities for improving the biological realism within existing modules, and for including modules representing additional biological components and mechanisms. We are currently implementing virus-scavenging by immune cells and local antiviral resistance due to Type 1-IFN paracrine signaling by epithelial cells. We are calibrating the virus replication module to existing experimental data for SARS-CoV-2 and influenza A. Because different tissues within the body have different responses to local viral infection, developing our framework to support the modeling and simulation of multi-organ disease progression (e.g., by identifying model parameters corresponding to specific tissues and physiological compartments) would allow us to understand the highly variable whole-body progression of many viral diseases.

The immune response to viral infection depends on locus of infection, degree of infection and patient immune state. Understanding the reasons for immune failure to contain infection, or pathological responses like cytokine storms or sepsis, requires models of immune response at multiple locations and scales. The same is true for understanding and predicting the possible protective or adverse effects of coinfection. The number of permutations of infection timing and combination of pathogens is too large to address purely by experiments, but could be addressed by simulations. Spatial modeling is also important because the spatiotemporal dynamics of coinfection within tissues may be important to the outcome (e.g., whether individual cells can be superinfected, whether viral lesions with a tissue are disjoint or overlap, whether the main foci of the pathogens are in the same or different tissues).

We can also study the systemic effects of possible therapies with known molecular modes of action (as seen in Results). Evaluating therapies in a simulated context prior to performing animal or human trials could lead to more effective and rapid drug discovery and to optimized dosage and timing of treatments. Understanding the origins of population variability in disease progression is crucial to providing optimal personalized treatment. While the simulations presented here begin with a single infected cell, a simulation which begins with multiple infected cells might better represent the infection dynamics of patients that have been subject to high level exposure, such as healthcare workers. Factors such as hypertension, immunosuppression and diabetes affect tissue state and immune response and could also be incorporated into our model. More detailed studies of these factors using our model could reveal more about the effects of population variability (due to age, genetic variation, prior drug treatment or immune status) on disease progression. Such computational studies could be accomplished using recently published concomitant, calibrated ODE-based simulations of COVID-19 treatment [61] .

We are working to implement validated non-spatial models of viral infection and immune response as agent-based spatial models (e.g., viral production, cytokine secretion, tissue damage). By starting with a validated model that uses ordinary differential equations and adding spatial components gradually, we can calibrate our spatial models and validate our results. In ongoing work, we have developed a formal method for spatializing ODE models and employing their parameters such that these analogous spatial models reproduce the ODE results in the limit of large diffusion constants. Using this method, we can combine the ability to do rapid formal parameter identification of ODE-based models and to leverage published ODE model parameters with the flexibility of spatial modeling. In these cases, any differences between the ODE results and the spatial model can be definitively attributed to spatialization (e.g. the local spread of virus or cytokine or the limited speed of movement of immune cells), or to additional factors which are difficult to include in an ODE model (e.g., the variability of individual cells or the complex time course of virus release by individual infected cells). We have developed formal spatializations of a number of interesting ODE models of COVID, such as [61, 62] , to explore the effects of stochasticity of outcomes, the effects of spatial mechanisms, and infection dynamics at a particular site of infection on the predictions of these models. An additional benefit of our approach is that we can easily and consistently combine and integrate ODE models which focus on different aspects of the complex process of infection, spread and clearance (e.g., combining published models of intracellular INF-induced viral resistance with spatial models of plaque spread in vitro [20, 63] ). We illustrate both of these strengths in Integration of an explicit RNA synthesis model allows the spatiotemporal modelling of hepatitis C virus infection, where we integrate a published HCV model of subgenomic replication into our framework. We can also conduct simultaneous, cross-platform validation of spatial models by building multiple implementations of the same conceptual and quantitative models on independent modeling platforms (here Chaste [64, 65] and Morpheus [66] ).

The COVID-19 crisis has shown that drug discovery and therapy development both require new predictive capabilities that improve their effectiveness and efficiency. We have developed our framework to explore the relationship between molecular, cellular-level and systemic mechanisms and outcomes of acute viral infections like SARS-CoV-2, and to support development of optimal, patient-specific treatments to combat existing and new viruses.

In this section we first present our model as a high-level conceptual model where we list each process included in an implementation-independent manner. We then detail the quantitative model and its computational implementation, which uses a Cellular Potts representation of cellular dynamics. All quantitative models are implemented in a modular, extensible simulation architecture built using the CompuCell3D simulation environment, which is publicly available for download and further development by interested members of the scientific community (see S2 Text).

As discussed in Introduction (Fig 2A) we consider viral propagation in an epithelial tissue and a lymph node. The tissue contains two interacting spatial components: an epithelium component (consisting of a monolayer of epithelial cells), and an extracellular environment component (containing immune cells, extracellular virus and chemicals). The lymph node component (whose state is affected by signaling from the tissue) adds immune cells to the extracellular space when in a proinflammatory state and removes them when in an antiinflammatory state. A set of processes and interactions govern how the states of these components evolve in time. We detail these components, processes and interactions in the following subsections and in Fig 2. Epithelium component. The epithelium component of the model represents the layer of epithelium in the tissue, and is composed of epithelial cells of four types: uninfected, infected, virus-releasing and dead (Fig 2C) . We assume that the epithelial cells are immobile. We implicitly model the ECM by considering its influence on all processes in the epithelium component. The epithelial cells contain modules that describe the viral life cycle and approximate the amount of virus as a continuous quantity (Fig 2B) , including: binding and internalization of viral particles from the extracellular environment (E1), intracellular replication (E2) and release (E3) of synthesized virus into the extracellular environment, as well as cell death caused by viral-replication-associated damage, immune-cell killing and oxidative agent killing (E4).

E1-Viral internalization. Module E1 models extracellular virus binding to epithelial cell receptors and internalization (including endocytosis-dependent and -independent routes). Internalization of viral particles involves binding of the viral spike protein to target cell-surface receptors, truncation by surface proteins and receptor-mediated endocytosis or fusion with the host plasma membrane. We assume the dynamics of internalization can be represented by describing the dynamics of virus-surface-receptor binding, determined by the amount of extracellular virus and target surface receptors, and by the binding affinity between them (T1!E1). We also consider the dynamic depletion of unbound target surface receptors on a cell when it internalizes a virus and superinfection of infected cells. Internalized viral particles initiate the viral replication process (E1!E2 and Fig 2B) .

E2-Viral replication. Module E2 models the viral replication cycle inside a host epithelial cell (Fig 2B) . Individual cells infected with many non-lytic viruses show a characteristic threephasic pattern in their rate of viral release. After infection and during an eclipse phase, a cell accumulates but does not yet release newly assembled viruses. In a second phase, the rate of viral release increases exponentially until the virus-releasing cell either dies or, in a third phase, saturates its rate of virus synthesis and release. Viral replication hijacks host synthesis pathways and may be limited by the availability of resources (amino acids, ATP, etc.), synthesis capability (ribosomes, endoplasmic reticulum, etc.) or intracellular viral suppression. A quantitative model of viral replication needs to be constructed and parameterized such that it reproduces these three phases.

We model viral replication based on processes associated with positive sense singlestranded RNA (+ssRNA) viruses. +ssRNA viruses initiate replication after unpacking of the viral genetic material and proteins into the cytosol (E1!E2). The viral RNA-dependent RNA polymerase transcribes a negative RNA strand from the positive RNA strand, which is used as a template to produce more RNA strands (denoted by "Viral Genome Replication" in Fig 2B) . Replication of the viral genome is the only exponential amplification step in the growth of most viruses within cells. Subgenomic sequences are then translated to produce viral proteins ("Protein Synthesis" Fig 2B) . Positive RNA strands and viral proteins are transported to the endoplasmic reticulum (ER) where they are packaged for release. After replication, newly synthesized viral genetic material is translated into new capsid protein and assembled into new viral particles ("Assembly and Packaging" in Fig 2B) . These newly assembled viral particles initiate the viral release process (E2!E3). We assume the viral replication cycle can be modeled by defining four stages: unpacking, viral genome replication, protein synthesis, and assembly and packaging. Fig 2B illustrates these subprocesses of replication and their relation to viral internalization and release.

E3-Viral release. Module E3 models intracellular transport of newly assembled virions and release into the extracellular environment (E3!T1 and Fig 2B "Release") . We conceptualize the virus being released into the extracellular fluid above the apical surfaces of epithelial cells.

Newly assembled virions are packed into vesicles and transported to the cell membrane for release into the extracellular environment (E2!E3). We assume that no regulation occurs after assembly of new virus particles, and that release into the extracellular environment can be modeled as a single-step process (E3!T1).

E4-Cell death. Module E4 models death of epithelial cells due to various mechanisms. Models the combined effect of the many types of virus-induced cell death (e.g., production of viral proteins interferes with the host cell's metabolic, regulatory and delivery pathways) as occurring due to a high number of assembled viral particles in the viral replication cycle (E2!E4). Models cell death due to contact cytotoxicity (I3!E4). Models cell death due to oxidizing cytotoxicity (T3!E4).

Extracellular environment component. The extracellular environment contains the immune cells, extracellular virus, cytokines and oxidative agent, and is the space where transport of viral particles (T1), cytokine molecules (T2) and the oxidizing agent (T3) occurs. We implicitly model the ECM in the extracellular environment by subsuming its geometrical, biochemical and biophysical influences on immune cell motility and virus/cytokine/agent spreading in the chosen rate laws and parameter set. Immune cells are mobile and can be either activated or inactive (I1). Inactive immune cells move through random cell motility and activated immune cells chemotax along the cytokine field (I2). The immune cell modules also account for cytotoxic effects of immune cells on contact due to antigen recognition (I3) and through the secretion of oxidizing agents (I4).

T1-Viral transport. Module T1 models diffusion of viral particles in the extracellular environment and their decay. Viral particles are transported by different mechanisms (e.g., ciliated active transport, diffusion) and media (e.g., air, mucus) at different physiological locations and through different types of tissue (e.g., nasopharyngeal track, lung bronchi and alveoli). Viral particles are eliminated by a variety of biological mechanisms. We represent these mechanisms by modeling transport of viral particles as a diffusive virus field with decay in the extracellular environment. We model transport in a thin layer above the apical surfaces of epithelial cells. Viral internalization results in the transport of a finite amount of virus from the extracellular environment into a cell and depends on the amount of local extracellular virus and number of cell surface receptors (T1-E1). Infected cells release viral particles into the extracellular environment as a result of the viral replication cycle (E3-T1).

T2-Cytokine transport. Module T2 models diffusion and clearance of immune signaling molecules in the extracellular environment. The immune response involves multiple signaling molecules acting upon different signaling pathways. We assume that the complexity of immune signaling can be functionally represented using a single chemical field that diffuses and decays in the extracellular environment. Once infected, epithelial cells secrete signaling molecules to alert the immune system (E2-T2). Locally, exposure to cytokine signaling results in activation of immune cells (T2-I1). Upon activation, immune cells migrate towards infection sites guided by cytokine (T2-I2). Lastly, activated immune cells amplify the immune signal by secreting additional cytokines into the extracellular environment (I1-T2). We model long-range effects by assuming that cytokine exfiltrates tissues and is transported to immune recruitment sites (T2-L1).

T3-Oxidizing agent burst and transport. Module T3 models diffusion and clearance of a general oxidizing agent in the extracellular environment. One of the cytotoxic mechanisms of immune cells is the release of different oxidizing agents, reactive oxygen species like H 2 O 2 and nitric oxide. The mechanism of action of such agents varies but we assume that we can generalize such effects by modeling a single diffusive and decaying oxidizing agent field in the extracellular environment. The oxidizing agent is secreted by activated immune cells after persistent exposure to cytokine signals (I4!T3). We assume that the range of action of the oxidizing agent is short. Cell death is induced in uninfected, infected and virus-releasing epithelial cells when sufficiently exposed to the oxidizing agent (T3!E4).

I1-Immune cell activation. Module I1 models immune cell maturation due to cytokine signaling. Immune cells mature at the recruitment site before being transported to the infection site as inactive immune cells (L1!Immune Cells). After infiltration, immune cells need to be exposed to local cytokine signals before activating (T2!I1). Once activated, immune cells chemotax along the cytokine field (I2) and amplify immune signaling by releasing cytokine molecules into the extracellular environment (I1!T2). Immune cells can also deactivate after a period of activation (I1 and Fig 2C) .

I2-Immune cell chemotaxis. Module I2 models activated immune cell chemotactic migration towards infection sites. We assume that upon activation (I1!I2), immune cells move preferentially towards higher concentrations of local cytokine (T2!I2).

I3-Immune cell direct cytotoxicity and bystander effect. Module I3 models immune cell cytotoxicity when immune cells (both activated and inactive) identify and induce death in epithelial cells with internal virus. Immune cells identify epithelial cells with internal virus on contact by antigen recognition and induce cell death by activating the caspase cascade (I3!E4). Uninfected, infected, and virus-releasing epithelial cells in contact with an epithelial cell that is killed by direct cytotoxicity can die through a bystander effect.

I4-Immune cell oxidizing agent cytotoxicity. Module I4 models activated immune cell killing of target cells through the release of a diffusive and decaying oxidizing agent into the environment. Cell death is induced in uninfected, infected and virus-releasing epithelial cells when sufficiently exposed to the oxidizing agent (T3!E4).

Lymph node component. The lymph node component models the net pro-or antiinflammatory state of the immune system. It responds to cytokines received from the tissue and adds or removes immune cells from the tissue (L1).

L1-Immune cell recruitment. Module L1 models immune cell recruitment and infiltration into the tissue in response to cytokine signaling by infected cells and activated immune cells. Infected cells secrete signaling molecules into the extracellular environment (E2!T3), which alerts resident immune cells and recruits new immune cells from the blood, distant lymph nodes and bone marrow. We model the local strength of the cytokine signal as causing an increase in the strength of the signal at the immune recruiting sites. We model long-distance signaling by assuming that cytokine molecules in the extracellular environment exfiltrate the infection site and are transported through the lymphatic system to the lymphatic system to lymph nodes and through the bloodstream to initiate immune-cell recruitment (T2!L1). A delay on the order of minutes to hours would represent semi-local recruitment (e.g., at the blood vessels). A delay on the order of days would represent long-range, systemic recruitment (e.g., the time required for a dendritic cell to reach a lymph node and an induced T cell to return). Recruited immune cells are then transported and infiltrate the infection site (L1!Immune Cell).

For model construction and integration we use the open-source multicellular modeling environment CompuCell3D (www.compucell3d.org) which allows rapid and compact specification of cells, diffusing fields and biochemical networks using Python and the Antimony language [55, 67] . CompuCell3D is specifically designed to separate model specification (conceptual and quantitative models) from the details of model implementation as a simulation and to make simulation specification accessible to biologists and others not specializing in software development. In this paper, we specifically designed the Python modules and their cross-scale integration to have clear and stable APIs, allowing modules to be rapidly swapped out by collaborating developers. CompuCell3D runs on Windows, Mac and Linux platforms without change of model specification, and allows cluster execution for parameter exploration.

Cellular Potts model (CPM). Cell types. Cells are divided into two broad groups, epithelial and immune cells, and have a type (see Fig 2C) which determines their properties, the processes and interactions in which they participate, and their events and dynamics. Epithelial cells can have one of four types (uninfected, infected, virus releasing and dead). Immune cells can have one of two types (activated and inactive). Cell types can change according to outcomes of various modules, and a module specifying such an event describes both the initial and final cell types of the transition. A cell type in the model is not a phenotype in the biological sense (e.g., epithelial cell), but an identifier for the various states that a particular cell can assume (e.g., dead epithelial cell). When an epithelial cell changes to the dead type, all epithelial modules are disabled and the cell is generally inactive.

Cellular dynamics. Cellular spatial dynamics is modeled using the Cellular Potts model (also known as CPM, or Glazier-Graner-Hogeweg model), which represents generalized cells and medium as occupying a set of sites in a lattice [68] . Random cell motility is modeled as the stochastic exchange of sites at intercellular and cell-medium interfaces. Configurations evolve to minimize the system's effective energy H,

Here σ is the integer identification of a cell and τ(σ) is the type of cell σ. v(σ) and V(σ) are the current and target volumes of cell σ, respectively, and λ volume is a volume constraint coefficient. N(x) is the neighborhood of site x, δ i,j is the Kronecker-delta, and J(τ, τ 0 ) is the effective contact energy per unit surface area between cells of types τ and τ 0 . The final term, H chemotaxis ; models chemotaxis-directed cell motility, and is prescribed by module I2. The cell configuration evolves through asynchronous lattice-site copy attempts. A lattice-site copy attempt starts by random selection of a site x in the lattice as a target, and a site x 0 in its neighborhood as a source. A configuration update is then proposed in which the value x 0 from the source site overwrites the value of x in the target site. The change in total effective energy DH due to the copy attempt is calculated, and the update is executed with a probability given by a Boltzmann acceptance function,

Here the intrinsic random motility H � controls the stochasticity of accepted copy attempts. Updates that reduce the system's effective energy are always accepted. The unit of simulation time is the Monte Carlo step (MCS)-taken to be 20 minutes in this work. One MCS corresponds to considering a number of copy attempts equal to the number of lattice sites.

Epithelial component modules. Processes E1-E4 describe epithelial cell functions as defined below. E1, E2 and E4 govern the cell-type transitions of epithelial cells (see Fig 15) . E1 transforms an uninfected epithelial cell into an infected epithelial cell. E2 transforms an infected epithelial cell into a virus-releasing epithelial cell. E4 transforms a virus-releasing epithelial cell into a dead cell.

E1-Viral internalization. To capture the stochasticity associated with internalization of discrete virus particles in terms of discrete binding events, we assign each uninfected, infected and virus-releasing epithelial cell a probability of absorbing diffusive viral particles from the extracellular virus field (T1). The uptake probability Pr(Uptake(σ)>0) for each cell σ is given by a Hill equation of the total amount of diffusive viral particles in the domain of the cell c vir (σ), the number of unbound cell surface receptors SR(σ) and the binding affinity between them.

Here h upt is a Hill coefficient, R o is the cell's initial number of unbound receptors, k on is the virus-receptor association affinity, k off is the virus-receptor dissociation affinity, α upt is a characteristic time constant of uptake and Δt is the time represented by one MCS. At each simulation time step, the uptake probability is evaluated against a uniformly-distributed random variable. When uptake occurs, the uptake rate is proportional to the local amount in the virus field (T1), and the probability of uptake is used to define the amount (Uptake) of virus taken up during the MCS,

The amount absorbed by each cell (Uptake) is uniformly subtracted from the virus field over the cell's domain and the cell's number of cell surface receptors is reduced by the same E2-Viral replication. Our viral replication model combines equations and parameters from several sources to represent the replication of a generic virus [7, 9, 52, 53] . The model contains four variables representing viral quantities in different states of the viral replication process: internalized virus U (Eq (6), the process of unpacking), viral genome taking part in genomic replication R (Eq (7), the process of viral genome replication), synthesized protein P (Eq (8), the process of protein synthesis), and assembled and packaged virions A (Eq (9), the process of assembly and packaging). Biologically, the only process which can exponentially increase the rate of virus production by a single cell is viral genome replication, so the equations include the positive feedback by R in Eq (7) . Biologically, factors like the cell's metabolism, limited number of ribosomes, maximum rate of endoplasmic reticulum function and activation of intracellular viral suppression pathways all limit production of viral components, so we include a Michaelis-type saturation term for the rate of replication in Eq (7). Each uninfected, infected and virus-releasing cell in the simulation contains an independent copy of the system of ordinary differential equations modeling the viral replication process,

dP dt ¼ r t R À r p P; ð8Þ

Here r u is the unpacking rate, r max is the viral replication rate, r t is the translating rate (rate at which viral genomes turn into RNA templates for protein production) and r p is the packaging rate. Uptake is defined in E1 and Release is defined in E3. The saturation of the rate of viral genome replication is represented by a Michaelis-Menten function, r half Rþr half , where r half is the amount of R at which the viral genome replication rate is reduced to r max 2 (and the flux is reduced to 1 2 r max r half Þ. The duration of the eclipse phase of single-cell infection (the time between the first entry of the virus into the cell and the first release of newly synthesized virus) is approximately t eclipse � 1 r u þ 1 r max þ 1 r t þ 1 r p (11.7 hours for the reference parameter set in Table 1 ), with the additional complication that in our model, an epithelial cell does not release virus until A reaches a threshold value. The timescale for tenfold increase of virus release when viral replication is maximal is t 10 � logð10Þ r max (7.7 hours for the reference parameter set in Table 1 ). The number of newly assembled virions is passed to the cell's instance of the viral release module (E3). See Fig 2B for 

Here r s is the release rate of viral particles and A is the level of assembled virus in the cell (defined in E2). The total amount released by each cell r s AΔt is subtracted from the cell's state variable for assembled virions A and passed to the source term of the extracellular virus field (T1) to maintain mass balance.

E4-Cell death. For cell death due to virally-induced apoptosis, each infected and virusreleasing cell can die due to the amount of intracellular virus. The rate of death is defined as a stochastic function of the state variable for assembled new virions from the viral replication module (E2). If a virus releasing cell dies then it changes its cell type to dead and the cell's instances of the viral internalization, replication and release modules are disabled. The probability of virus-induced apoptosis per simulation step is a Hill equation of the current load of assembled virus,

where A(σ) is the number of assembled virions in cell σ, h apo is a Hill coefficient, V apo is the amount of assembled virions at which the apoptosis probability is 0.5 per unit time and α apo is a characteristic time constant of virally-induced apoptosis. For modeling of cell death due to contact cytotoxicity, see I3-Immune cell direct cytotoxicity and bystander effect. For modeling of cell death due to oxidizing cytotoxicity, see I4-Immune cell oxidizing agent cytotoxicity. 

Regardless of the death mechanism the internally assembled virions are not released to the environment and do not take action in further infection. We assume that in the process of death the assembled virus and viral particles are either damaged or deactivated. Lymph node modules. L1-Immune cell recruitment. The total immune cell population is governed by an ordinary differential equation of a dimensionless state variable S that represents immune response due to local conditions and long-distance signaling. Our convention is that when S>0, immune cells are recruited to the simulation domain; likewise, immune cells are removed from the simulation domain when S<0. Probability functions of S describe the likelihood of immune cell seeding and removal,

Prðremove immune cellÞ ¼ erfðÀ a immune SÞ; S < 0:

Here the coefficient α immune is the sensitivity of immune cell addition and removal to the state variable S. The dynamics of S are cast such that, in a homeostatic condition, a typical number of immune cells can be found in the simulation domain, and production of cytokine (T2) results in additional recruitment via long-distance signaling (i.e., with some delay). We model this homeostatic feature using the feedback mechanism of the total number of immune cells N immune in the simulation domain. Cytokine signaling is modeled as perturbing the homeostatic state using the term α sig δ. Here δ is the total amount of decayed cytokine in the simulation domain and α sig >0 models signaling by transmission of cytokine to some far-away source of immune cells. We write the rate of change of S as

Here β add and β sub control the number of immune cells in the simulation domain under homeostatic conditions. β delay controls the delay between transmission of the cytokine to the lymph node and corresponding immune response by adjusting the rate of recruitment due to total cytokine (i.e., increasing β delay increases the resulting delay). β decay regulates the return of S to an unperturbed state (i.e., S = 0, increasing β decay increases the rate of return to S = 0). To determine the seeding location, the simulation space is randomly sampled n seeding times, and an immune cell is seeded at the unoccupied location with the highest amount of the virus field. If no location is unoccupied, then the immune cell is not seeded. The removal probability is evaluated for each immune cell at each simulation step. Immune cells are removed by setting their volume constraint to zero.

Immune cell modules. The four processes I1-I4 capture immune cell functions which are defined below. These processes control how immune cells are activated, translocate, and kill other cells. Their interactions with epithelial cells and other model components are illustrated in Fig 17. I1-Immune cell activation. Inactive immune cells become activated with a probability according to a Hill equation of the total cytokine bound to the cell B cyt (σ,t), After ten hours, an activated immune cell becomes inactive, in which case evaluations of activation (Eq (15)) recommence. The immune cells "forget" a percentage (1−ρ cyt ) of the bound cytokine each time step while taking up an amount of cytokine from the environment (ω cyt (τ(σ),t) defined in T2), B cyt ðs; tÞ ¼ r cyt B cyt ðs; t À DtÞ þ o cyt ðtðsðxÞÞ; tÞ:

I2-Immune cell chemotaxis. Activated immune cells experience a motile force as a response to a signaling field. Immune cells chemotax along concentration gradients of the cytokine field. The chemotactic effective energy H chemotaxis associated with the gradient is calculated according to a chemotactic sensitivity parameter λ chemotaxis and calculated chemotactic force F chemotaxis . The contribution of H chemotaxis to the change in the system's total effective energy is calculated using F chemotaxis when considering copy attempts. The chemotactic force at a location x is saturated by normalizing the chemotactic sensitivity parameter by the concentration of cytokine at the center of mass of the cell at x, c cyt,CM (σ(x)),

I3-Immune cell direct cytotoxicity and bystander effect. Immune cells, whether activated or not, have the ability to kill infected cells by direct contact. At each simulation step, immune cells trigger cell death in the infected and virus-releasing epithelial cells with which they come in contact. When an infected cell is killed by direct cytotoxicity, each of its first-order neighbors is evaluated for cell death by a bystander effect model. Each of those neighbors σ 0 2N(σ) in the first-order neighborhood N(σ) of a cell σ killed by direct cytotoxicity has a probability k bystander of dying from the bystander effect given by, Prðtðs 0 ; tÞ ! DeadjDirect CytotoxicityðsÞ ¼ TrueÞ ¼ k bystander 8s 0 2 NðsÞ: ð18Þ 

I4-Immune cell oxidizing agent cytotoxicity. Immune cells release a short-range, diffusive oxidizing agent when exposed to high cytokine concentration (T3). The oxidizing agent kills an epithelial cell of any type when the total amount of oxidizing agent in the domain of the cell c oxi (σ) exceeds a threshold for death t death oxi , Prðtðs; tÞ ! Deadjc oxi ðsÞ > t death oxi Þ ¼ 1:

Extracellular environment modules. T1-Viral transport. The change in concentration of the virus field c vir is calculated at each location in the simulation domain by solving the following reaction-diffusion equation,

Here D vir is the diffusion constant of the extracellular virus and γ vir is the decay rate is the decay rate. Uptake and release by a cell at each location are determined using the viral internalization (E1) and the viral release (E3) modules, and are uniformly applied over all sites of the domain of the cell.

T2-Cytokine transport. The change in concentration of the cytokine field c cyt is obtained by solving a reaction-diffusion equation of the following general form,

The decay term γ cyt c cyt represents cytokine leaving the simulation domain (e.g., in immune recruitment). To model immune signaling, the rate of cytokine secretion is described by an increasing Hill function of c sig (σ(x)) with Hill exponent h cyt = 2. The meaning of c sig (σ(x)) depends on the cell type and the Hill exponent can differ for other cell types and states. The rate of cytokine secretion s cyt is written as, Here σ cyt (τ(σ(x),t)) is the maximum cytokine secretion rate for the cell type at x, c sig (σ(x)) is a quantity that affects the cells cytokine secretion, ω cyt (τ(σ(x),t))is the cytokine uptake rate of the cell type at x and V cyt (τ(σ(x),t)) is a dissociation coefficient of cytokine secretion for the cell type at x. σ cyt is nonzero for infected epithelial cells, virus-releasing epithelial cells and activated immune cells. For infected and virus-releasing epithelial cells c sig is the amount of assembled virus A in the viral replication module, and for activated immune cells c sig is the total amount of cytokine in the domain of the cell. Similarly, for epithelial cells V cyt is the amount of assembled virus, and for immune cells V cyt is the amount of cytokine in the domain of the cell. ω cyt (τ (σ(x),t)) is constant and nonzero for activated and inactive immune cells.

T3-Oxidizing agent transport. The oxidizing agent field diffuses according to the reactiondiffusion equation,

Bursts of oxidizing agent are implemented as a source term for one time step according to a rate coefficient σ oxi , which is uniformly mapped onto the source term s oxi over the domain of each activated immune cell. An oxidizing burst occurs in immune cells with an activated state when the total cytokine in the immune cell's domain exceeds a threshold t sec oxi . Initial and boundary conditions. The domain of all simulations had dimensions of 90 x 90 x 2 lattice sites. The initial cell configuration consisted of a 30 x 30 sheet of uninfected epithelial cells, each of size 3 x 3, on the lower layer of lattice sites (see S18 Fig for a demonstration of the negligible effects of a non-uniform arrangement of epithelial cells). Epithelial cells were "frozen", in that they were immobile, leaving the remaining 90 x 90 subdomain for occupancy by recruited immune cells. For cellular dynamics and mass transport, periodic boundary conditions were applied in the plane of the epithelial sheet, and Neumann conditions were applied along the direction orthogonal to the epithelial sheet. All field solutions for the diffusive viral, cytokine and oxidizing agent fields were initialized as zero everywhere.

At each first simulation step, the epithelial cell in the center of the sheet was set to infected, and the amount of internalized virus U of the viral replication model was set to a value of one. All epithelial cells were initialized with a number of unbound surface receptors SR = R o . All immune cells, when introduced to the simulation by recruitment, were initialized in an inactive state, and with a bound cytokine value equal to zero (B cyt = 0). During transition of an uninfected epithelial cell to the infected type, all state variables of the viral replication model were initialized with a value of zero.

Simulation specifications. Model implementation and all simulations were performed using CompuCell3D, which uses a non-dimensional lattice for CPM-based cellular dynamics and non-dimensional explicit time integration of reaction-diffusion field solutions. As such, a baseline parameter set was constructed for all CPM parameters and modules developed in this work (Table 1) . Non-dimensionalization was performed on model parameters for a lattice dimension of 4 μm per pixel along each dimension, at 20 minutes (1/3 hours) per MCS. All replicas were simulated for ten trials, each 1,000 MCS (20000 minutes, 333 hours, 14 days) long. Simulation data was collected at a frequency of 10 MCSs (200 minutes, 3 hours) for all simulations. Table. Varying parameters in simulations shown in Fig 4. Virus-receptor association affinity and immune response delay coefficient shown for no immune response (Fig 4A) , widespread infection (Fig 4B) , slowed infection (Fig 4C) , containment (Fig 4D) , recurrence ( Fig  4E) and clearance ( Fig 4F) . All other parameters are as in Table 1 

Remdesivir for severe acute respiratory syndrome coronavirus 2 causing COVID-19: An evaluation of the evidence

Clinical course and outcome of 107 patients infected with the novel coronavirus, SARS-CoV-2, discharged from two hospitals in Wuhan

A systems approach to infectious disease

Big Data's Role in Precision Public Health. Front Public Health

A primer on multiscale modelling of infectious disease systems

Within-host infectious disease models accommodating cellular coinfection, with an application to influenza †

Rapid community-driven development of a SARS-CoV-2 tissue simulator

-host Modelling of COVID-19 Kinetics in Humans

A model of the binding, entry, uncoating, and RNA synthesis of Semliki Forest virus in baby hamster kidney (BHK-21) cells

Brownian adhesive dynamics (BRAD) for simulating the receptor-mediated binding of viruses

The transcriptional cycle of HIV-1 in real-time and live cells

Modeling rotavirus-like particles production in a baculovirus expression vector system: Infection kinetics, baculovirus DNA replication, mRNA synthesis and protein production

Forces during bacteriophage DNA packaging and ejection

A minimal kinetic model for a viral DNA packaging machine

Spatial aspects of HIV infection

Physiologically based pharmacokinetic models of small molecules and therapeutic antibodies: a mini-review on fundamental concepts and applications

Agent-based modeling of host-pathogen systems: The successes and challenges

Probing the effects of the well-mixed assumption on viral infection dynamics

Dose Optimization of Chloroquine by Pharmacokinetic Modeling During Pregnancy for the Treatment of Zika Virus Infection

Design considerations in building in silico equivalents of common experimental influenza virus assays

Early chest computed tomography to diagnose COVID-19 from suspected patients: A multicenter retrospective study

COVID-19 severity correlates with airway epithelium-immune cell interactions identified by single-cell analysis

Imaging features and evolution on CT in 100 COVID-19 pneumonia patients in Wuhan, China

The pathology of influenza virus infections

Spatiotemporal Dynamics of Virus Infection Spreading in Tissues

Spatiotemporal Dynamics of a Generalized Viral Infection Model with Distributed Delays and CTL Immune Response

Towards whole-lung simulations of aerosol deposition: A model of the deep lung

Spatiotemporal Dynamics of HIV Propagation

Agent-Based Modeling of Endotoxin-Induced Acute Inflammatory Response in Human Blood Leukocytes

Characterizing the Dynamics of CD4+ T Cell Priming within a Lymph Node

A hybrid stochastic-deterministic computational model accurately describes spatial dynamics and virus diffusion in HIV-1 growth competition assay

Dynamics of HIV Infection: A cellular automata approach

Experimental and computational analyses reveal that environmental restrictions shape HIV-1 spread in 3D cultures

Dynamically linking influenza virus infection with lung injury to predict disease severity. bioRxiv

Spatial models of virus-immune dynamics

COVID-19: unravelling the clinical progression of nature's virtually perfect biological weapon

Global landscape of mouse and human cytokine transcriptional regulation

Pathogen Recognition and Innate Immunity

Host Factors in Positive-Strand RNA Virus Genome Replication

Complement fragment C5a and inflammatory cytokines in neutrophil recruitment during intramammary infection with Escherichia coli

Cytokines. 2. Cytokines and lung inflammation: mechanisms of neutrophil recruitment to the lung

Lipid-Cytokine-Chemokine Cascade Drives Neutrophil Recruitment in a Murine Model of Inflammatory Arthritis

Epithelial cells: structure, transport, and barrier function

Validated Models of Immune Response to Virus Infection

Influenza A virus infection kinetics: quantitative data and models

Viral Load Affects the Immune Response to HBV in Mice With Humanized Immune System and Liver

Role of Regulatory T Cells during Virus Infection

T cell responses to viral infections-opportunities for Peptide vaccination

Interferon activation and innate immunity

Cellular and molecular immunology

Ebola Immunity: Gaining a Winning Position in Lightning Chess

Optimal Replication of Poliovirus within Cells. The American Naturalist

Kinetic Modeling of Virus Growth in Cells

Validation of Cycloserine Efficacy in Treatment of Multidrug-Resistant and Extensively Drug-Resistant Tuberculosis in Beijing, China

Multi-scale modeling of tissues using CompuCell3D

Severe Acute Respiratory Syndrome Coronavirus 2 from Patient with Coronavirus Disease, United States

Heritage and early history of the boundary element method

Fisher waves in an epidemic model

Mathematical modeling of subgenomic hepatitis C virus replication in Huh-7 cells

Therapeutic efficacy of the small molecule GS-5734 against Ebola virus in rhesus monkeys

Timing of Antiviral Treatment Initiation is Critical to Reduce SARS-CoV-2 Viral Load

Potency and timing of antiviral therapy as determinants of duration of SARS CoV-2 shedding and intensity of inflammatory response

Mathematical Modeling of RNA Virus Sensing Pathways Reveals Paracrine Signaling as the Primary Factor Regulating Excessive Cytokine Production

Chaste: an open source C ++ library for computational physiology and biology

Chaste: A testdriven approach to software development for biological modelling

Morpheus: a user-friendly modeling environment for multiscale and multicellular systems biology

Antimony: a modular model definition language

Simulation of biological cell sorting using a two-dimensional extended Potts model

Apoptosis by influenza viruses correlates with efficiency of viral mRNA synthesis

A tunable diffusion-consumption mechanism of cytokine propagation enables plasticity in cell-to-cell communication in the immune system

Competition for IL-2 between regulatory and effector T cells to chisel immune responses

Distinct T cell dynamics in lymph nodes during the induction of tolerance and immunity

Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Micro-and macrorheology of mucus

We would like to thank Ramray Bhat, Ryan Danehy, Javier Toledo, Bard Ermentrout, Geoffrey Fox, Yi Jiang, Paul Macklin, Jairaj Mathur, Roeland Merks, Ericka Mochan, Jason Shoemaker and Amber Smith for useful discussions and ideas and for comments on the manuscript.