A Modular Framework for Multiscale Spatial Modeling of Viral Infection and Immune Response in Epithelial Tissue

The COVID-19 crisis has shown that classic sequential models for scientific research are too slow and do not easily encourage multidisciplinary scientific collaboration. The need to rapidly understand the causes of differing infection outcomes and vulnerabilities, to provide mechanistic frameworks for the interpretation of experimental and clinical data and to suggest drug and therapeutic targets and to design optimized personalized interventions all require the development of detailed predictive quantitative models of all aspects of COVID-19. Many of these models will require the use of common submodels describing specific aspects of infection ( e.g. , viral replication) but combine them in novel configurations. As a contribution to this development and as a proof-of-concept for some components of these models, we present a multi-layered 2D multiscale, multi-cell model and associated computer simulations of the infection of epithelial tissue by a virus, the proliferation and spread of the virus, the cellular immune response and tissue damage. Our initial, proof-of-concept model is built of modular components to allow it to be easily extended and adapted to describe specific viral infections, tissue types and immune responses. Immediately after a cell becomes infected, the virus replicates inside the cell. After an eclipse period, the infected cells start shedding diffusing infectious virus, infecting nearby cells, and secretes a short-diffusing cytokine signal. Neighboring cells can take up the diffusing extracellular virus and become infected. The cytokine signal calls for immune cells from a simple model of the systemic immune response. These immune cells chemotax and activate within the tissue in response to the cytokine profile. Activated immune cells can kill underlying epithelial cells directly or by secreting a short-diffusible toxic chemical. Infected cells can also die by apoptosis due to the stress of viral replication. We do not include direct cytokine mediated protective factors in the tissue or distinguish the complexity of the immune response in this simple model. Despite unrealistically fast viral production and immune response, the current base model allows us to define three parameter regimes, where the immune system rapidly controls the virus, where it controls the virus after extensive tissue damage, and where the virus escapes control and infects and kills all


Introduction Motivation
Viruses are obligate intracellular parasites that exist almost everywhere on the planet. Viruses require a host to replicate and have the ability to infect humans, animals, plants, fungi and bacteria. They are considered to be the most abundant biological entity on the planet. Viruses cause a number of cancers ( e.g. , human papillomavirus (HPV) causes cervical cancer) and infectious diseases such as rabies, herpes, flu, common cold, severe acute respiratory syndrome ( 19). COVID-19 has now altered the environment, and socioeconomic status of human life. The current global outbreak of COVID-19 has reinforced the need for rapid and sufficient characterization and prediction of the pathways for highly contagious viral related disease.

Biological Background
Betacoronavirus Life Cycle Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an enveloped virus with a positive sense, single stranded RNA genome (ssRNA+). ssRNA + viruses have mutation rates that are millions of times higher than that of their hosts [1] .This characteristic allows ssRNA + viruses to have enhanced virulence and evolution rates, which makes this particular category of virus complex [2,3] . The genomic evolutionary history of SARS-CoV-2 shows lower mutation rates than other ssRNA + viruses such as influenza [3][4][5][6][7] . Its genome encodes 16 non-structural proteins, 4 structural proteins and 7 accessory proteins (NCBI accession #: NC_045512) [8] . SARS-CoV-2 displays two key spike proteins (S1 and S2) on its surface that are vital for viral attachment and internalization [9,10] . The S1 spike protein is responsible for the initial binding to host cells, whereas the S2 spike protein is activated by proteolytic cleavage by host cell proteases that mediate host cell-viral membrane fusion [11] . SARS-CoV-2 shares 80% sequence identity with the previously identified SARS-CoV, the virus identified as the causative agent of SARS in 2003 [2,8] . Although inherently different, both SARS-CoV and SARS-CoV-2 retain similar key amino acid residues which lead to similar infection pathways such as primary binding domain to host cell surface receptor angiotensin-converting enzyme 2 (ACE2) [12,13] .

Mode of Entry
Viral transmission and attachment are the first steps of a viral infection [14] . Principal modes of transport and transmission of RNA + sense strand virus particles include cilia, air, mucus, and water [15] . Zoonotic origin, as well as, high transmission rates of SARS-CoV-2 through human-human interactions has led to widespread global infection [16,17] . Both SARS-CoV and SARS-CoV-2 can infect host cells by two different modes of entry 1) binding to the host cell surface receptor angiotensin-converting enzyme 2 (ACE2), triggering pH-dependent proteolytic activation of S2 spike protein by cathepsin L resulting in endocytosis of the viral particle or by 2) binding to the ACE2 receptor and triggering pH-independent proteolytic activation of the S2 spike protein by the transmembrane serine protease TMPRSS2 that often resides near the ACE2 receptor on the host cell surface resulting in host cell-viral membrane fusion [12,[18][19][20] . Viral internalization occurs in the receptor-binding domain of host cells via endocytosis and/or fusion [21][22][23] .
Target Host Epithelial Cells SARS-CoV-2 follows a disseminated infection pathway that has shown the ability to locally infect mammalian respiratory tissue [17, [24][25][26] . Epithelial tissue lines the outer surface of organs and blood vessels throughout the body [27] . As well as the inner surface of cavities in many internal organs including the lungs [28] . Mammalian respiratory passageways from the naval cavity through the bronchi are lined by ciliated, columnar epithelium [28] . While lung alveoli are lined by a thin layer of simple squamous epithelium [29] . In the lung, epithelial cells seperate the airways and potentially harmful materials within them from the bloodstream, while allowing for the free diffusion of carbon dioxide and oxygen [28,29] . Viral entry occurs in the receptor-binding domain of host epithelial cells via endocytosis and/or fusion [21][22][23] . Establishing that transport of viral proteins and internalization of SARS-CoV-2 must occur first at the epithelial tissue layer in the respiratory system before spreading to inner tissue. / Viral Replication Viral replication of SARS-CoV-2 and other enveloped single stranded positive sense RNA viruses occurs within the cytoplasm of host epithelial cells after unpacking of viral proteins from the internalized particle [30] . An RNA-dependent RNA polymerase (RdRp) encoded within the virus genome produces a negative strand from positive host RNA strand by hijacking host epithelial cell translation machinery. The negative viral RNA strand is then used to produce more positive RNA strands and smaller positive strand subgenomic sequences [31] . Subgenomic sequences are translated to produce viral proteins [31] . After replication inside the host, viral positive RNA strands and viral proteins are transported to the Endoplasmic Reticulum (ER) where they are assembled to form new virions. New virions are then packaged into vesicles and transported to the cell membrane for release into the extracellular environment through the continuous process of viral budding [32][33][34] . The yield for newly produced virions is not clearly defined in SARS-CoV-2 as it is with lytic viruses [35] . The main regulation of the viral replication process occurs at the replication stage, because the balance between replication and translation must be carefully maintained in order to achieve maximum rates of replication and translation without exhaustion of host cell resources [33,34] .
Triggering Innate Immune Response The innate immune system is the first response to encounter SARS-CoV-2 infection. Production of SARS-CoV-2 viral proteins interfere with various host cell metabolic, regulatory and delivery pathways [36,37] . The detection of foreign viral proteins and the hijacking of host cell translation machinery triggers immune system gene expression [38] . This triggered genetic expression signals the production of type I interferons (IFNs), cytokines, chemokines, interleukins, myeloids, and small molecule RNA signals through a series of toll-like receptors (TLRS) [39] . This process occurs in order to alert local cell populations of infection through small molecule signaling and to stimulate the growth and recruitment of immune cells [40] . The immune cells that are targeted for early activation in the innate immune response are dendritic cells, macrophages, neutrophils, mast cells, basophils, eosinophils, leukocytes, and natural killer (NK) cells [28] . The stimulated activation of immune cells such as macrophages and monocytes interplays the recruitment, initiation, and maturation of NK cells and naive T cells from stem cells in bone marrow to the site of local infection [40] . Simultaneously, programmed cell death (apoptosis) is initiated in cells with hijacked translation machinery through caspase activation and through contact with NK cells [41][42][43] . After immune cell recruitment, mature monocytes, macrophages, and T cell populations begin to deal with cells infected by virus by direct killing [44][45][46] . This initial innate immune response that takes place within hours and days of infection involves multiple signaling molecules acting upon multiple complex signaling pathways [47,48] . Cytokines produced in this time period are associated with early clinical symptoms of viral infection such as fever, fatigue, and cough [49] .
Existing kinetic viral replication models address how viruses play a significant role in human health and have the ability to change population trends and dynamics [50,51] ; however, it remains a challenge to predict how novel viruses interact with their host cellular environments. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent of the current global pandemic and existing kinetic viral replication models do not address the need for a specific multicellular model representing the infection pathway.

Modeling Background
Context Mathematical models and computer simulations have been extensively used to study in-host progression of viral infection. Kinetics approaches are commonly used to model different stages of viral replication cycle such as binding and internalization [52,53] , replication and translation [54,55] , assembly, packing and release [56,57] [70,71] have been modeled using agent-based spatial approaches. With respect to the family of beta coronaviruses spatial models have been developed motivated by the offset of the ongoing pandemic [72] .
Our Approach to an Initial Proof-of-Concept Model In this paper we propose a modular framework to model interactions between generalized epithelial, immune cells and their extracellular environment during viral infection. The model can be used to develop and interrogate hypotheses related to the spatiotemporal dynamics of viral SARS-CoV-2 infection of critical nasopharyngeal and lung tissue and model Covid-19 progression. The framework is intended to serve as a base model for constructing and implementing more advanced models of targeted cellular-and intracellular-level phenomena in tissue after initial exposure. In its current state, it has not been formally peer-reviewed, and should not be used for patient diagnostics or predicting clinical outcomes.
The model and its implementation can be employed and developed to interrogate questions and mechanistic hypotheses about the spread of a virus and how the interplay between viral spreading and immune response determine the outcome of the disease, such as: • Why does the progression of the disease seem to be dependent on the initial viral exposure level? • Why is the start time of symptoms and immune response so variable? • What is the role of cytokine signaling in explaining immune response variability?
• What are the specific factors and key players determining the offset of early immune response? The model includes a representation of extracellular virus in the mucus, epithelial cells and immune cells. It also includes the processes of epithelial cell infection by extracellular virus, viral replication and cell damage inside epithelial cells, release of viruses by infected epithelial cells, immune cell response to infected epithelial cells and immune cell killing of infected and non-infected epithelial cells.

Results
In this section we present results from simulations intended to interrogate the parameter space of the model framework (see Models and Methods ) with respect to select, critical parameters of interest to understanding SARS-CoV-2 and viral infection, in general. All simulations were performed according to specifications described in Simulation Specification , and results described in that same section were recorded concerning spatial, population, and system-level metrics for each simulation trial of each parameter set presented here. All simulations were performed using CompuCell3D [73] .

Model Captures Spatial Dynamics of Viral Infection and Spread in Epithelial Tissue
We first determined a baseline set of parameters (Table 1) for which widespread infection occurs over the course of simulation time. We initiate our simulation with a sheet of uninfected epithelial cells, no immune cells, diffusive virus, or diffusive cytokine, and a sheet of one initially infected cell at the center of the simulation domain ( Figure 1A). As viral replication progresses, the initial infected cell starts secreting virus into the extracellular environment. Neighboring cells are infected by viruses in their extracellular environment and initiate their own viral replication cycles. Infected cells also secrete cytokines, initiating the recruitment of immune cells by exfiltration of cytokine from the simulation domain ( Figure 1B).
The reference simulation shows significant immune cell recruitment between 400 and 800 Monte Carlo steps (MCS), but the virus spreads uncontrolled, leading to the eventual infection and death of all epithelial cells, At time 200 MCSs infection has already spread to neighboring cells. Immune cells have already been recruited to the tissue, and dead cells not at the initial site of infection indicate likely cytotoxic killing. At time 400 MCSs massive cell death has occured centered about the initial site of infection. Immune cell recruitment from increased levels of diffusive cytokine has more than doubled the number of immune cells in the simulation.
By 600 MCSs infection has reached the simulation domain, indicating total spread of the infection. Immune cells continue to be recruited even to 800 MCSs, even though most cells were dead, due to the delay in immune response. By time 1,000 MCSs, all cells were dead and many immune cells had exited the simulation. The number of infected cells was maximal at 700 MCSs, the amount of extracellular virus was maximal at 860 MCSs. The extracellular cytokine level was maximal at 740 MCSs. As the number of live cells decreases, the amount of extracellular virus and the amount of extracellular proinflammatory cytokine decreases. The immune recruitment signal peaks at 740 MCSs (1480 minutes), leading to infiltration of immune cells too late to contain the spread of infection (maximal number of immune cells at 870 MCSs, 1740 minutes). For all simulation parameters see Table 1. Code to execute this specific simulation for these parameters in CompuCell3D is available at https://github.com/covid-tissue-models/covid-tissue-response-models/tree/cc3d_first_model_v0 _cand/CC3D/Models/BiocIU/SARSCoV2MultiscaleVTM/interesting-results/Preliminary%20Set% 201/simimg .  Critical parameters control outcomes of the simulation. Increase in the virus-receptor association affinity , a parameter that controls the internalization of extracellular viral k on particles into epithelial cells, drives the system toward quantitative distincts simulation outcomes ( Figure 2). We generated results using an established baseline value and a 100-fold k on increase from the baseline. Baseline shows a low number of infected and dead cells at the k on end of the simulation (Figure 2A Variations in critical and parameters drive the system toward two distinct β delay k on simulation outcomes: containment and widespread infection (red and shaded regions). , a β delay parameter that controls the responsiveness of immune cell recruitment to local diffusive cytokine in the simulation domain. We generated simulations using established baseline values for k on and and compared results with simulations generated using 10-fold and 100-fold β delay increase and decrease from the baseline values ( Figure 3). Parameter variation produces distinctive qualitative outcomes in the number of dead cells ( Figure 3A). Baseline values of parameters result in propagation of the infection through the tissue and widespread death of epithelial cells ( Figure 3A, center panel). Increasing and decreasing produces drives k on β delay the system towards widespread infection and cell death with small variability between simulation replicated ( Figure 3A, red shaded panels). In this region, decreasing decreases the rate β delay of cell death as shown by the inflection points on the number of dead cells occurring earlier for smaller values of . β delay Decreasing and increasing drives the system towards containment of the k on β delay initial infection, as shown by early saturation in the number of dead cells ( Figure 3A, green shaded panels). In this region, increasing controls the rate at which containment is β delay achieved, as shown by saturation being arrived at earlier for bigger values of . Also in this β delay region, controls the total number and variability of infected and dead cells before k on containment whenever is not too large, as shown by the saturation value being higher for β delay increasing values of ( Figure 3A, first green shaded column). For the remaining regions, the k on final outcome of the simulation is undetermined, since the number of dead cells is still changing at the end of the simulation.
The distinctive qualitative outcomes of the simulation can also be observed in the dynamics of viral diffusion in the extracellular environment ( Figure 3B). In the widespread infection region, there is a regime of fast viral growth followed by a regime of viral clearance (FIgure 3B red shaded region). In this region, controls how fast the maximum concentration k on of virus in the extracellular environment is achieved during the growth regime, as shown by the maximum value being reached at earlier for increasing values of . The parameter k on β delay controls the amplitude of the rebound after a short clearance phase, as shown by the distance between two maximums in the viral concentration curve increasing for increasing values of . β delay In the containment region, the same two qualitative regimes of viral growth can be observed but with quantitative different properties ( Figure 3B red shaded region). An initial growth of virus is followed by clearance of virus. The maximum concentration of the viral field is / smaller in the containment region than the widespread infection region. In these region, β delay controls the clearance rate of the extracellular virus as shown by the rapid decay with increasing values of . The time scales of viral growth in the extracellular environment is faster than β delay the rate of cell death, as shown by maxima in the extracellular environment occurring before maxima in the number of dead cells.
Increase in the Rate of the Immune Response and Decrease in Internalization Drive the System Towards Containment of Infection The phase space for key output metrics reveal that decreasing viral internalization and increasing immune response consistently drives the system towards containment of the infection (Figure 4). High viral internalization is generally correlated with high number of final dead cells ( Figure 4A), high number of maximum viruses ( Figure 4B) and high number of maximum cytokine levels ( Figure 4C). In general, values for each of these metrics increases proportionally with increased internalization, with the lowest value providing the lower bound and the highest value providing the upper bound. A noticeable exception is the max cytokine level, where at least one of the lower viral internalization values (0.1) produces a high number of cytokines ( Figure 4C).
Rapid immune response is generally associated with lower number of final dead cells, maximum virus and maximum cytokine in the extracellular environment. For fixed viral internalization values, each of these output metrics decreases with increase in the rate of immune response. In general, the output metrics behave as decreasing saturation functions on the immune response, with the exception of the cytokine levels for one of the lower viral internalization values (0.1), where there is a rapid growth in cytokine levels followed by a long decay ( Figure 4C).
Clearance time measures the earliest time from the end of the simulation at which no infected cells are found inside the simulation domain ( Figure 4D). Higher viral internalization is generally associated with smaller clearance times for immune response rates below the baseline and with higher clearance times for immune response rate above the baseline. With respect to the immune response, clearance times show varying regimes, but in general there is a slow increase in clearance times for smaller immune response rates, and a sharp decrease in clearance times for higher immune response rates. Viral internalization controls at which level of immune response rate the sharp decline in clearance times starts. The decline in clearance times occurs at slower immune response rate for viral internalization below the baseline and higher immune response rates for values above the baseline.
Variation in the Rate of Viral Replication Show Difference Outcomes of the Diseases showed a bounded range of rmax parameter values that produces spread of infection. For low viral replication rates, immune cells were recruited that successfully eliminated the infected cell before viral secretion could spread infection to neighboring cells. For high viral replication rates, excessive internalized viral particles due to fast replication induced apoptosis before spread of the virus to neighboring cells. A: Number of dead cells over simulation time in ten trials of each variation (top row) and the mean value of results (bottom row) with blue shading demarcating one standard deviation in each direction. Large variance in the number of dead cells for trials with the variation was due to that all trials either produced 0r 1 max complete cell death or nearly none at all, reflecting that the parameter set with this variation is near a bifurcation of the regimes of certain massive cell death, and viral containment due to excessive viral replication rate. B: Total diffusive virus over simulation time for ten trials of each variation. Marginal diffusive virus was measured for downward variations. C: Number of infected cells over simulation time for ten trials of each variation. D: Mean number of activated immune cells over simulation time for ten trials of each variation, with blue shading demarcating one / standard deviation in each direction. No activated immune cells were measured for the most upward variation, indicating death of the initially infected cell was due to viral apoptosis, rather than by cytotoxicity.
In our second parameter sweep, the replication rate ( , from the viral replication r max model) was varied two orders of magnitude upward and downward by factors of ten while holding all other values of the baseline parameter set constant. In both variations of slower replication rate ( i.e. , and ), marginal amounts of diffusive viral particles were r 10 1 − max r 10 2 − max measured, followed by the activation of a few immune cells (Fig. 5). Once no immune cells in these simulations were in the activated state, no additional activity in the simulations were measured concerning viral infection.
For the first upward variation in replication rate ( ), all cells in simulation died within 0r 1 max 300 MCSs in most simulations. However, in two simulations of the first upward variation, the initially infected cell died before immune cells were even recruited to the simulation domain.
This same phenomenon was observed for all trials of the second upward variation ( ). r 10 2 max Without immune cells, the only mechanism of the model framework available to kill a cell is that of viral apoptosis.

Discussion and Future Perspectives
We have produced a spatial model and computer simulation of epithelial tissue that includes key aspects of viral infection, replication and immune response. This model, while highly simplified, does provide a number of very suggestive results, including showing that the viral production peaks before tissue damage as seen in diseases like influenza and that the control of immune response and viral spread can lead to conditions in which the virus is rapidly controlled, controlled after substantial damage or uncontrolled. We have studied the possible effect of a drug that reduces the rate of RNA polymerization and shown that it can lead to improved viral control if given early. A key next step would be to determine how its effectiveness decreases when it is administered progressively later after infection.
We also proposed a framework for collaborative rapid development of models of viral infection and immune response. We recognize that our overall model architecture is incomplete in many respects, from the number and types of immune cells and cytokines to the processes we are including ( e.g. , the cytokine signals from an infected cell can alert uninfected cells to viral challenge and increase their probability of apoptosis or reduce their rate of virus production after infection, macrophage can scavenge extracellular virus, immune cells can be infected by virus and change their cytokine relaying). We are working with a wide variety of biological researchers to determine the additional components we should prioritize to add to our high-level model structure Our model is open-source and organized in modules that are free to extend, reuse and adapt. We would like to develop model sharing workflows and tools to facilitate parallel, independent submodel development. The key to such collaborative development is to adopt guidelines for model specification and sharing: specify models and submodels at the conceptual and quantitative levels with well-defined inputs, outputs and validation metrics and data and to implement them as independently executable callable components, which can be independently validated. Some of the directions for future collaborative developments are: • Improvement of the biological realism of core submodels. In our preliminary model, many subcomponents are proof-of-concept versions, which we can be replaced with more supportable models. Accessible areas of improvements will be more realistic models of viral internalization and replication (including timing) and models immune cell response ( e.g. , scavenging of virus by macrophages, local IFN-gamma signaling from infected epithelial cells to non-infected cells to reduce their rate of viral production). • Identification of model parameters corresponding to specific critical tissue (nasopharyngeal, alveolar) and physiological compartments (throat, upper respiratory and lower respiratory tracts). Building compartment models that can be connected using models of transport of virus, cytokines and immune cells between through the body. • Use the models to do basic disease mechanism studies to explore the ability of the model to simulate immune clearance under a variety of different situations (effect of initial viral load and locus on disease progression in a patch, immune excitability, rate of delay in signaling to lymph nodes). Identify key periods when infection progression is susceptible to particular types of perturbation ( e.g. , slowing viral replication after 3 days generally has little effect on outcome in influenza models). • Studying the systemic effects of possible therapies with known molecular modes of action, e.g. , remdesivir is an analog of the nucleoside adenosine and targets the RNA-dependent RNA polymerase and blocks viral RNA synthesis [75]. This would reduce the rate of viral replication but also might affect the intracellular infection response mechanisms in apparently paradoxical ways. An antibody therapy would introduce humoral response mechanisms earlier in the infection. Model the effect of immune stimulation or repression at different stages of infection to optimize dosage and timing. • Studying the origins of population variability in disease progression, by modeling the effects of hypertension, immunosuppression, diabetes and by considering how they change the timing of critical outcomes. • Benchmark against established nonspatial models of infection, immune response and clearance. We are working to implement established and validated nonspatial models of viral infection and immune response and replace specific boxes in these models ( e.g. , viral production, cytokine secretion or tissue damage) with the appropriate spatial components of our model. By starting with a validated ODE and adding spatial components gradually, we can calibrate our spatial models and validate our results. • Conduct simultaneous validation by building multiple implementations of the conceptual and quantitative models independently and simultaneously. We are currently working with teams in Australia and Germany to develop independent multicellular implementations of our underlying conceptual models to validate them in parallel with our own implementations. / The COVID-19 crisis has shown that classic sequential models for scientific research are too slow and do not easily encourage transdisciplinary scientific collaboration. Our open-source project is designed to allow simultaneous component development and deployment by many groups and to minimize the overhead required for collaborative development. Model specification is compact and uses simple Python scripts to make model development accessible to those with critical disciplinary knowledge but without extensive computational experience. Transparency and validation in model development are particularly critical in a disease response model especially when that model is being developed rapidly and in parallel. We are already working closely with the CHASTE and Morpheus multicell modeling teams to develop validating replications of our base conceptual model. This near-simultaneous model replication and its lessons for model design and implementation will be the subject of a future paper.
We are reaching out to any community members to encourage them to extend, replace, improve and develop the entire model or any model components. Multiple teams using multiple methods to study a single topic has proven powerful in other areas of science and we want to replicate that effective approach to rapid scientific progress with modular virtual tissue modeling.
We are specifically eager to collaborate and support groups modeling viral replication, cell death due to viral replication, local cytokine signalling effects and systemic immune response. We are also eager to help experimental and drug discovery and therapy development teams adapt and refine this base model for their specific applications and to work with groups with relevant experimental validation data for the integrated model and model components. We would also be happy to discuss approaches to integrating this model framework as a component of teams developing whole body and population level models. Many submodels and model parametrization may be appropriate for replacement with AI surrogates. Finally, we are eager to work with groups who want to replicate these results using different quantitative and computational methodologies or develop improved approaches to conceptual and quantitative model specification.

Models and Methods
Here we model multiscale cell-cell transport interactions in a local infected population of cells. We generalize these mechanisms by representing transport of viral particles as a diffusive chemical field in the extracellular environment. Viral transport and internalization of viral particles is the necessary first step of our viral infection model. In our model, we approximate these discrete processes using a continuous kinetics approach, determined by association and dissociation constants, the number of available cell surface receptors and the amount of viral particles in the extracellular environment.Therefore, in our model we use the relevant aspects of transport on the thin layer above the apical surface of epithelial cells as the initial site of viral transport and attachment to represent an early local respiratory tissue infection. Where the infection of susceptible epithelial cells occurs when the diffusive viral field comes in contact with the cell surface. In our model, we represent the complexity of viral replication by defining 4 broad stages: unpacking, replication, translation and assembly. Therefore, in our model the kinetics of viral replication determine the release of new viral particles to the extracellular environment that contribute to the further spreading of the virus in host tissue.
For virally-induced apoptotic processes in our model, each cell is given a probability of dying associated with the number of assembled viral particles inside the cell. Cytokines were chosen to represent the larger group of small molecule signals that include chemokines, interferons, and RNAi. For the purpose of this early infection model, we functionally represented the complexity of immune signaling by using a single chemical field diffusing in the extracellular environment. The extracellular field can produce local immune effects such as activation of immune cells after immune cells are exposed to a cytokine signal for a period of time. The extracellular field can also produce long range immune effects, by correlating the concentration of cytokines in tissue to the strength of signal at the immune cell recruiting sites. We represented the transport of cytokines through the lymphatic system and bloodstream by introducing delay terms in the appearance of immune cells at the infection site. These delay terms may also correlate with the variation in response time between host immune cells and programmed cell death. Cell death also occurs in our model due to various mechanisms associated with host immune response.

Conceptual Model: Biological Hypotheses and Assumptions
Epithelial Cell Level At the epithelial cell level, the model accounts for binding and endocytosis of viral particles, intracellular replication and exocytosis to the extracellular environment, as well as for induced apoptosis from viral replication associated damage.
E1 -Viral internalization : model of extracellular virus binding to epithelial cell receptors, endocytosis and release of viral genetic material into the cytoplasm. 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 captured by focusing on the dynamics of virus-surface receptor binding, determined by the densities of extracellular virus and target surface receptors, and the binding affinity between them ( T1-E1 ). Internalized, viral particles initiate the viral replication process ( E1-E2 ).
E2 -Viral replication : model of viral replication cycle inside the host cell. Single stranded positive RNA viruses can initiate replication after unpacking of viral genetic material and proteins into the cytosol. The viral RNA-dependent RNA polymerase transcribes a negative RNA strand from the positive RNA strand. This negative strand is used as a template to produce more positive RNA strands and smaller positive strand subgenomic sequences. Subgenomic sequences are then translated to produce viral proteins. Positive RNA strands and viral proteins / are transported to the ER where they are packed for release. We assume the viral replication cycle can be modeled by defining four replication stages: unpacking, replicating, translating and assembling.
Internalized viral particles are disassembled at the unpacking stage ( E1-E2). Viral replication hijacks some of the host metabolic pathways and is limited by the availability of resources in the host cell. We assume we can model the rate limiting effect of resource availability as regulation at the replicating step. After replication, newly synthesized viral genetic material is translated into new capsid protein and assembled into new viral particles. These newly assembled viral particles initiate the viral release process ( E2-E3 ).
E3 -Viral release : model of intracellular transport of newly assembled virions and exocytosis into the extracellular environment. After assembly inside the host, newly packed virions are transported to the ER where they 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 assembling of new virions particles and that exocytosis into the extracellular environment can be modeled as a single-step process ( E3-T1 ).

Immune Cell Level
At the immune cell level, the model accounts for activation and chemotaxis of immune cells due to cytokine signaling and the cytotoxic effects of immune cells on infected epithelial cells due to antigen recognition or oxidative agents.

I1 -Immune cell activation : model of immune cells maturation due to cytokine signaling.
Immune cells mature at the recruitment site before being transported to the infection site. We assume however that upon infiltration, immune cells need to be exposed to local cytokine signals before exhibiting active immune cell behavior ( T2-I1 ). Once activated, immune cells amplify immune signaling by releasing cytokine molecules into the extracellular environment ( I1-T2 ).

Extracellular Environment Level
At the tissue level, the model accounts for the extracellular transport of viral particles, cytokine signaling molecules, and an oxidative burst agent.
T1 -Viral transport : model of diffusion and spreading of viral particles in the extracellular environment. Viral particles are transported by different mechanisms (ciliated active transport diffusion) and mediums (air, mucus) at different physiological locations and through different types of tissue (airway, nasopharyngeal track, lung). We assume that we can generalize these mechanisms by representing transport of viral particles as a diffusive chemical field in the extracellular environment. We model transport on a thin layer above the apical surfaces of epithelial cells where viral particles are deposited and transported. Infection of susceptible cells occurs when the diffusive viral field comes into contact with the cell surface ( T1-E1 ). Infected cells release viral particles to the extracellular environment as a result of the viral replication cycle ( E3-T1 ).
T2 -Cytokine transport : model of diffusion of small immune signaling molecules in the extracellular environment. Immune response involves multiple signaling molecules acting upon different signaling pathways, but we assume that the complexity of immune signaling can be functionally represented using a single chemical field diffusing in the extracellular environment. Once infected, epithelial cells secrete signaling molecules to alert the immune system ( E3-T2 ). Cytokine signaling has both local and distant effects. Locally, exposure to cytokine signaling results in activation of newly recruited immune cells occurs ( T2-I1 ). Upon activation, immune cells further infiltrate the tissue towards infection sites guided by the cytokine molecules ( T2-I2 ).
Lastly, active immune cells amplify the immune signaling by further secreting cytokines into the extracellular environment ( I1-T2 ). We model long range effects by assuming cytokine exfiltrate tissue and is transported to immune recruitment sites ( T2-L1 ). We assume that the local strength of the cytokine signal tissue correlates to the strength of the signal at the immune recruiting sites. We model transport of cytokines through the lymphatic system and bloodstream with delays to account for exfiltration and recruitment.
T3 -Oxidative Burst Agent: model of diffusion of a general oxidative agent. One of the cytotoxic mechanisms of immune cells is the release of different oxidative agents, reactive oxygen species (ROS) like H 2 O 2 and nitric oxide (NO). The mechanism of action of such agents vary depending on the agent but we assume we can generalize such effects by modeling a single oxidative diffusing field in the extracellular environment. The oxidative agent is secreted by active immune cells after persistent exposure to cytokine signals ( I4-T3 ). We assume that the range of action of the oxidative agent is short. Cell death is induced in epithelial cells when they come into contact with the oxidative agent ( T3-E4 ). Figure 6: Conceptual model. Schematic representation of the model objects, compartments, processes and interactions. Epithelial and immune cells refer to the two agents in the model. Extracellular environment and lymph node refer to compartments where interactions between agents and fields take place. Each agent has associated submodels that dictate their behaviors. Epithelial cells have viral internalization (E1), viral replication (E2), cell death (E4) and viral release (E4) submodels. Immune cells have activation (I1), chemotaxis (I2), contact cytotoxicity (I3) and oxidative cytotoxicity (I4) submodules. Fields describe transport of material in the extracellular environment and to the lymph nodes. Three fields characterize the model: viral field (T1), cytokine field (T2) and oxidative agen field (T3). At the lymph node compartment, transport of cytokines feeds into an immune recruitment submodel (L1).

Quantitative Model and Implementation
For our model construction and integration we use the open-source multicell 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 [73,74] . 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 specialized in software development. In this work we have specifically designed the Python modules and their cross-scale integration to have clear APIs, allowing the model elements to be rapidly swapped out by collaborating developers. CompuCell3D runs on Windows, Mac and Linux platforms without change of model specification. Recent versions allow cluster execution for parameter exploration.

Cellular Potts Model (CPM) Cell Types
Cells are divided into two classes-epithelial and immune-and assigned a phenotype by which various submodels behave. These phenotypes can change according to outcomes of various submodels, and a submodel specifying such an event describes both the initial and final phenotypes of the transition, as well as the conditions of its occurrence. As such, a cell phenotype in the model framework is not a phenotype in the biological sense ( e.g. , epithelial cell), but rather serves as an identifier for the various states that a particular cell class can take ( e.g ., dead epithelial cell) due to events defined by the submodels. Epithelial cells can adopt one of three different phenotypes: uninfected, infected and dead. The specific behaviors of each cell phenotype is defined per submodel as relevant to their purpose. When a cell changes to a dead type, all epithelial submodels are disabled the cell is generally inactive.

Cellular Dynamics
Cellular spatial dynamics is modeled using the Cellular Potts model (a.k.a., the CPM, Glazier-Graner-Hogeweg model), which represents generalized cells and medium as occupying a set of sites in a lattice [75] . Cell random motility is modeled as the stochastic exchanging of sites at the interface of cells and medium so to minimize the system's effective energy that ℋ governs various behaviors, Here is the identification of a cell and is the type of cell . and are the σ τ (σ) σ v (σ) V (σ) current and target volumes of cell , respectively, and is a volume constraint coefficient. σ λ volume is the neighborhood of site , is the Kronecker-delta, and is the contact N (x) x δ i, j J (τ , τ ) ′ effective energy between types and . The final term models directed motility by τ τ ′ ℋ chemotaxis chemotaxis, and is prescribed by submodels. For every spin flip attempt, a site in the lattice is x randomly selected, as is a site in its neighborhood. The change in the effective system x′ ℋ Δ energy is calculated due to the identification at being changed ( i.e. , "flipped") to the ℋ x identification at the neighborhood site , and the spin flip occurs with a probability according to x′ a Boltzmann acceptance function, Here the intrinsic random motility controls the stochasticity of spin flips, and spin flips that ℋ * reduce the effective system energy are always accepted. The unit of simulation time is the Monte Carlo step (MCS), which demarcates the accomplishment of having considered a number of spin flips equal to the number of lattice sites.

Epithelial Submodels Viral Internalization
Internalization is a discrete process by which a viral particle binds to one or more cell receptors. To capture the stochasticity associated with discrete binding events, we assign each uninfected and infected cells a probability of absorbing diffusive viral particles from the extracellular viral field. The uptake probability for each cell occurs according (U ptake(cell) ) Pr > 0 to a Hill equation of the total amount of diffusive viral particles in the domain of the cell , (cell) c vir the number of unbound cell surface receptors and the binding affinity between them. R(cell) S Here is a Hill coefficient, is the initial number of unbound cell receptors, is the h upt R o k on association constant between virus and cell surface receptors, is the dissociation constant k of f and is the cell volume. At each simulation time step the uptake probability is evaluated ol(cell) V against a uniformly distributed random variable. When uptake occurs, the uptake rate is proportional to the local amount of the viral field, and the probability of uptake is used to describe the efficiency by which uptake occurs, ptake. dt dSR(cell) =− U The amount absorbed by each cell is uniformly subtracted from the viral field over the cell's domain and the number of unbound receptors and passed to the cell's instance of the viral replication model according to conservation of species. We assumed that epithelial cells continue uptaking viral particles from the environment after infection until cell receptors are depleted.

Viral Replication
Our very simple proof-of-concept viral replication model was inspired by discussions with Paul Macklin and has a form similar to that published by Macklin but differs in equations and parameters [72] . It represents the replication of a generic virus and does not include several aspects of viral replication specific to coronaviruses and their timescales. The system of ordinary differential equations modeling the viral replication process is assigned as an independent copy to each uninfected and infected cell. The model contains four variables representing different states of the viral replication process: unpacking , replicating , U R packing , and assembly of new virion capsids . P A Here is the unpacking rate, is the maximum replication rate, is the translation rate and r u r max r t is the packing rate. The regulation of replication is represented by a Michaelis-Menten r p function of the amount of replicating viral material , where is the amount of at which r half R +r half r half R the replicating rate is . The viral replication model is specified as a readily sharable Antimony 2 rmax string that can be implemented as a standalone using the Tellurium package. The number of newly assembled virion capsids is passed to the cell's instance of the viral release model. where is the number of assembled virions, is a Hill coefficient and is the (cell) A h apo V apo amount of assembled virions at which the apoptosis probability is 0.5. Immune Submodels

Immune Cell Recruitment
The total immune cell population is governed by an ordinary differential equation of a state variable that represents immune response due to local conditions and long-distance S signaling. Our convention is that when , immune cells are recruited to the simulation S > 0 domain; likewise, immune cells are removed from the simulation domain when . We S < 0 accomplish this by imposing probability functions describing the likelihood of immune cell seeding and removal, Here the coefficient controls the sensitivity of immune cell addition and removal to the α immune state variable . The dynamics of are cast such that, in a homeostatic condition, a typical S S number of immune cells can be found in the simulation domain, and production of cytokine in the simulation domain results in additional recruitment via long-distance signaling ( i.e. , with some delay). We accomplish this by using the feedback mechanisms of the total number of immune cells in the simulation domain and a fraction of the total amount of decayed N immune cytokine .
Here is the total amount of decayed cytokine in the simulation domain and δ α sig δ models signaling by transmission of cytokine to some far-away source of immune 0 < α sig < 1 cells. With these mechanisms, we write the rate of as such, S Here and control the number of immune cells in the simulation domain under β add β sub homeostatic conditions, controls the delay between transmission of the cytokine and β delay immune response, and controls the return of to an unperturbed state ( i.e. , ). β decay S S = 0 At each simulation step the seeding probability is evaluated against a uniformly distributed random variable. To determine the seeding location, the simulation space is randomly sampled, and immune cells are seeded at the unoccupied location with the highest amount of the viral field. If no location is unoccupied, then the immune cell is not seeded. The removal probability is evaluated against a uniformly distributed random variable for each immune cell at each simulation step. Immune cells are removed by setting their volume constraint to zero.

Immune Cell Chemotaxis
Activated immune cells experience a motile force as a response to a signaling field. The immune cells chemotax on the chemical field representing cytokine signaling molecules. The chemotactic function measures the local gradient of the cytokine field and computes the effective energy associated with the gradient according to a prescribed chemotactic ℋ chemotaxis sensitivity parameter and calculated chemotactic force . The contribution of λ chemotaxis F chemotaxis to the change in the system total effective energy is calculated using when ℋ chemotaxis F chemotaxis considering spin flips. The chemotactic force is saturated by normalizing the chemotactic sensitivity parameter by the local concentration , (cell) c vir (15) ∇c . F chemotaxis = λ chemotaxis 1+c (cell) vir cyt Immune Cell Activation Immune cells have an associated boolean active state. All cells are initialized as inactive (Active= False). The activated state becomes true with a probability according to a Hill equation of the total cytokine bound to the cell , B cyt (immune cell, t) After one hour, an activated immune cell is deactivated, in which case evaluations of activation recommence. The immune cells "forget" a percentage the cytokine they have bound each time step,

Immune Cell Direct Cytotoxicity and Bystander Effect
Immune cells kill infected cells by direct contact. At each simulation step, neighbors of infected cells are evaluated. Apoptosis is triggered in an infected cell if it has an immune cell as one of its neighbors, in which case the cell changes type of dead. When an infected cell is killed by direct cytotoxicity, each of its first order neighbors is evaluated for bystander effect cytotoxicity. The neighbors have a probability of dying from bystander effect: Where the is the probability of a neighbor cell dying from bystander effect as a result of k bystander contact direct killing of an infected cell.

Immune Cell Oxidative Cytotoxicity
Immune cells when detecting a high cytokine concentration will release a short-range, diffusive oxidative agent. The oxidative agent kills any epithelial cell when its concentration inside the cell reaches a minimum concentration for death, . τ oxi death Figure 9: State diagram and interactions of epithelial cells . Immune cells can adopt two different generalized states: inactive and active. Inactive immune cells are recruited by the cytokine levels according to the immune recruitment submodel. Transition from inactive to active immune cells is determined by the immune activation submodel when cells are exposed to cytokine. Active immune cells amplify the cytokine signal by secreting cytokines to the extracellular environment. Active immune cells induce death of epithelial cells by direct cytotoxicity when coming into contact with infected cells, bystander effect by killing neighbors of infected cells and by releasing cytotoxic oxidative agents into the extracellular environment.

Mass Transport Submodels Viral Transport
The change in concentration of the viral field is calculated at each location in the c vir simulation domain by solving a reaction-diffusion equation using the following partial differential equation, Transport parameters such as the diffusion constant and decay rate are estimated Here is the dissociation coefficient. ζ

Oxidative Agent Transport
The oxidative agent field secreted by immune cells with an activated state diffuses according to the transport equation, (22) Δc c . ∂t ∂c oxi = D oxi oxi − γ oxi oxi + s oxi Bursts of oxidative agent are implemented as a source term for one time step at a rate of , which is uniformly mapped onto the source term . An oxidative (immune activated oxi) σ oxi − s oxi burst occurs in immune cells with an activated state when the cytokine in the immune cell's domain exceeds a threshold . τ oxi sec Initial and Boundary Conditions All simulations consisted of a domain of dimension 90 x 90 x 2 lattice sites. The initial cell configuration consisted of a 90 x 90 sheet of uninfected epithelial cells. Epithelial cells were "frozen", in that they were not permitted to translocate, 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 oxidative 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 its unpacking state variable of the viral replication model was set to a value of U one. All epithelial cells were initialized with a number of unbounded surface receptors equal R S to the number of initially unbound receptors . All immune cells, when introduced to the R o simulation by recruitment, were initialized not in an activated state, and with a boundy cytokine value ( ) equal to zero. During transition of an epithelial cell to the infected type, all state B cyt variables of the viral replication model were initialized with a value of zero. Secretion of viral particles by epithelial cells was only permitted for Infected-secreting types.

Simulation Specifications
Model implementation and all simulations were performed in 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 submodels developed in this work (Table 1). Non-dimensionalization was performed on all available model parameters in the literature for a lattice dimension of 4 μm per pixel along each dimension, at 120 s per MCS. For remaining model parameters, parameter estimation was performed such that, for the baseline set of parameters, spread of infection occurred throughout the domain by approximately the end of the simulation time. All parameter sets were simulated for ten trials, each consisting of 1,000 MCSs. Simulation data was collected at a frequency of 10 MCSs for all simulations, including the total number of all cell types, the total number of activated immune cells, the total diffusive virus and cytokine, and the value of the immune response signal ( ), S Two parameter sweeps were performed for submodel parameters of interest. In the first set, the virus-receptors affinity association and immune response delay coefficient k on β delay were varied. In the second parameter sweep, the replication rate was varied. For all r max variations, the baseline coefficient was multiplied by a factor in the set 10 -2 , 10 -1 , 10 0 , 10 1 , 10 2 . As such, the first parameter sweep consisted of 25 parameter sets, and the second parameter sweep consisted of 5 parameter sets, for a total of 300 total simulations.  Infection threshold 1 Initial target volume 9 Lambda volume λ volume 9 Initial number of immune cells 0 Lambda chemotaxis λ chemotaxis 1 Intrinsic Random Motility ℋ * 10 Contact coefficients (all interfaces) J 10 Table 1. Parameter values of baseline parameter set.