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

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.

2 infection, we studied the effects on progression of treatment potency and time-of-first 31 treatment after infection. In simulations, even a low potency therapy with a drug which 32 reduces the replication rate of viral RNA greatly decreases the total tissue damage and 33 virus burden when given near the beginning of infection. Many combinations of dosage 34 and treatment time lead to stochastic outcomes, with some simulation replicas showing 35 clearance or control (treatment success), while others show rapid infection of all epithelial 36 cells (treatment failure). Thus, while a high potency therapy usually is less effective when 37 given later, treatments at late times are occasionally effective. We illustrate how to extend 38 the platform to model specific virus types (e.g., hepatitis C) and add additional cellular 39 mechanisms (tissue recovery and variable cell susceptibility to infection), using our 40 software modules and publicly-available software repository. The model is used to investigate how potential treatments influence disease progression. 47

53
The current global pandemic of COVID-19, caused by the novel coronavirus 54 Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), has motivated the 55 study of beta coronavirus diseases at multiple spatial and temporal computational 56 modeling scales [1]. The time course, severity of symptoms and complications from 57 SARS-CoV-2 infection are highly variable from patient to patient [2]. Mathematical 58 modeling methods integrate the available host-and pathogen-level data on disease 59 dynamics that are required to understand the complex biology of infection and immune 60 response to optimize therapeutic interventions [3][4][5]. Mathematical models and computer 61 simulations built on spatial and ODE frameworks have been extensively used to study in-62 host progression of viral infection [6], with a recent acceleration in the development of 63 spatial COVID-19 viral infection models in response to the global pandemic [7,8]. 64 Building multiscale models of acute primary viral infection requires integration of 65 submodels of multiple biological components across scales (e.g., viral replication and 66 internalization, immune system responses). Non-spatial, coupled ordinary differential 67 equation (ODE) models can represent many aspects of pathogen-host interaction. In the 68 context of viral infection dynamics, specialized ODE models can describe both the entire 69 virus-host response at the tissue and organ levels and different stages of the viral 70 7 The temporal dynamics of the concentration of extracellular virus varies greatly 143 among virus families, tissues and host species [44]. However, for many viruses, including 144 influenza and coronaviruses, once infected cells begin to release virus, the amount of 145 extracellular virus increases exponentially over a period of a few days, reaching a peak 146 during an early phase of infection [45]. As the viral load increases, immune signaling 147 increases rapidly (this increase is associated with the onset of fever and other symptoms) 148 recruiting more circulating cells of the innate immune system to the infection site [46]. The model represents the extracellular environment as a space above the 217 epithelium which provides the space in which immune cells are recruited and move, and 218 into which cells release viruses and chemicals. We include a single type of immune cell 219 that exhibits many key immune-cell behaviors associated with macrophage, neutrophil, 220 NK cell and T-cell roles, including activation, chemotaxis, relaying and amplification of 221 cytokine signals, contact killing, secretion of diffusible killing factors, and bystander killing, 222 to represent the main types of immune-cell mediated cell killing, rather than any particular 223 immune cell phenotype. We omit macrophages, neutrophils and their phagocytosis and 224 signaling, which can be represented using simple extensions of the immune-cell type. We 225 do not model contact interactions between immune cells (e.g., by T-cells and APCs). We behaviors. Epithelial-cell processes include viral internalization (E1), viral replication (E2), viral release (E3) and cell 237 death (E4). Immune cell processes include activation (I1), chemotaxis (I2), contact cytotoxicity (I3) and oxidative 12 oxidizing-agent field (T3). Activated cells become inactive after 1 hour. The virus field (T1), cytokine field (T2) and 240 oxidizing-agent field (T3) describe spatially-varying densities of extracellular components. Field processes describe 241 diffusive transport and clearance of material in the extracellular environment and activated transport to the lymph nodes.

242
The lymph node is a single-compartment model whose pro-or anti-inflammatory state specifies the recruitment or 243 removal rate (L1) of immune cells in the epithelial tissue. The transport of cytokines to the lymph node promotes its We simulate extracellular-virus particle density as a continuum field and particle 259 transport and clearance as continuous diffusion and decay. We approximate the discrete 260 process of a cell's internalization of a virus particle by a stochastic virus internalization 261 event (E1) determined by time elapsed, the local concentration of the virus field, and the 262 number of available cell-surface receptors on the cell. We simplify the complexity of viral 263 replication into four steps: unpacking, viral genome replication, protein synthesis and 264 packaging/assembly (E2, Fig 2B) [7,[52][53][54]. The subcellular kinetics of viral replication 265 determine the rate of release of new virus particles into the extracellular environment (E3). 266 To represent the combined effect of the many types of virus-induced cell death, each 267 13 infected epithelial cell has a probability of dying that depends on the number of assembled 268 viral particles inside the cell per unit time (E4). 269 We simplify the complex biochemistry of the many molecular signals active in 270 epithelial tissues, which include chemokines, interferons and viral fragments, into a single 271 generic extracellular cytokine field that acts both as a tissue-local and systemic pro-272 inflammatory signal. Infected epithelial cells and immune cells both secrete cytokine (T2). 273 The cytokine field produces local immune effects such as activation (I1) and chemotaxis 274 (I2) of immune cells. Activated immune cells can revert back to inactive immune cells 275 when the cytokine signal ceases. The cytokine field also recruits immune cells to the 276 tissue through long-distance signaling via the lymphatic system (L1). We model 277 recruitment of immune cells to the simulation domain using an ordinary differential 278 equation for the immune signal ( ), which represents the balance between pro-and anti-279 inflammatory signaling and the delay due to antigen-presenting cell transport from the 280 tissue through the lymphatic system to the lymph node and due to the time required for 281 T-cell amplification. In the absence of infection, the lymph node maintains a small resident 282 immune cell population in the tissue. Immune cells can cause epithelial cell death (E4) by 283 three mechanisms: 1) contact cytotoxicity; 2) the bystander effect; and 3) through the 284 release of an oxidative agent. Immune cells kill infected epithelial cells by contact 285 cytotoxicity, in which case neighboring uninfected, infected and virus-releasing epithelial 286 cells can also die through a bystander effect (I3). In regions of the tissue with high cytokine 287 levels, immune cells secrete a diffusive oxidative agent to the environment (T3) that kills 288 uninfected, infected and virus-releasing epithelial cells (I4). 289 290 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020

291
We begin by presenting our base multicellular model of viral infection in an 292 epithelial tissue, along with a simulation for a baseline set of parameters and basic 293 analyses. We then explore the simulation's dependence on critical parameters that are 294 important to the dynamics of acute primary viral infection in airway epithelial cells. All 295 simulations and spatial, population and system-level metrics presented in this section 296 follow the specifications in Simulation Specifications. We performed simulations using the 297 open-source modeling environment CompuCell3D [55]. Downloading and running the 298 simulation provides instructions on how to run these simulations. 299 We initialize the simulations with periodic boundary conditions parallel to the plane 300 of the sheet and Neumann boundary conditions perpendicular to the plane of the sheet. 301 Initially, the extracellular environment does not contain any extracellular virus, cytokines, 302 oxidative agents or immune cells. We introduce infection by creating a single infected 303 epithelial cell at the center of the epithelial sheet, comparably to (but less than) initial 304 conditions of similar works that model discrete cellular distributions [18,31] . To illustrate 305 the full range of dynamics of viral infection in the presence of an immune response, we 306 established a baseline set of parameters (Table 1) for which the immune response is 307 strong enough to slow the spread of the infection, but insufficient to prevent widespread 308 infection and death of all epithelial cells (Fig 3  The infected cell immediately starts releasing cytokines into the extracellular 339 environment. After an incubation period (~150 minutes, 2 ½ hours), the first infected 340 epithelial cell (green) transitions from infected to virus releasing (red), and starts releasing 341 viruses into the extracellular environment. Initial release of extracellular virus causes 342 additional epithelial cells to become infected. Release of cytokines leads to delayed 343 addition of immune cells to the simulation domain (Fig 3D). By 4000 minutes (67 hours, 344 2 ¾ days), the number of infected cells increases 10-fold and the epithelial cells start 345 undergoing virally-induced death as the infection spreads radially outward from the initial 346 site. The increase in the number of infected cells and the local cytokine concentration is 347 accompanied by a similar increase in the local immune cell population. By 8000 minutes 348 18 (133 hours, 5 ½ days), the number of dead epithelial cells around the initial infection site 349 increases sharply. Following this phase of rapid cell death, the number of infected, virus-350 releasing and dead epithelial cells continues to increase exponentially but at a slower 351 rate. This transition results in the formation of an annular region of infected cells spreading 352 radially outwards from the initial infection site (Fig 3A) (of 1/900 per cell in our analysis), rendering the frequency of recurrence negligible 418 (Fig 4F). A special case of clearance (Failure to infect) occurs when the initially 419 infected epithelial cells fail to infect any other epithelial cells (Fig S2B).   Table 1  We ran ten simulation replicas for each parameter set, increasing and decreasing 471 the baseline parameter values 10-fold and 100-fold (Figs 5-7). For each simulation 472 replica, we examined the number of uninfected epithelial cells (Fig 5), the number of 473 infected epithelial cells (Fig 6), the total extracellular virus (Fig 7), the number of dead 474 epithelial cells (Fig S3), the number of virus-releasing cells (Fig S4)

512
. CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. The resulting cytokine secretion elicits a moderate to high response of the immune 546 recruitment signal (Fig S6) and number of immune cells (Fig S5). The early recruitment 547 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint 29 of immune cells leads to many epithelial cells dying before 4000 minutes (66 hours, 2 ⅔ 548 days) (Fig S3). For a high probability of virus internalization (high ), even low amounts 549 of extracellular virus are sufficient to cause recurrence. For moderate to low and 550 moderate to high (upper left unshaded subplots in Figs 5-7), rate of new infection 551 of epithelial cells is slightly faster than the immune system's response, resulting in a 552 combination of widespread infection, slowed infection and containment, and a few cases 553 of clearance. The immune system is only moderately responsive to the cytokine signal, 554 resulting in a slow to moderate increase in the immune recruitment signal (Fig S6) and in 555 the number of immune cells (Fig S5). Cases of clearance and containment occur for a 556 smaller final number of dead epithelial cells (Fig S3) compared to previously discussed 557 cases. 558 559 Even moderate inhibition of genomic replication by antiviral therapies significantly 560

reduces the spread of infection, but only when initiated soon after infection 561
Optimal therapeutic use of antiviral drugs requires considering the relationship 562 between molecular efficacy (how effectively the drug blocks a particular aspect of the viral 563 life cycle at saturation concentration), potency of therapy (the effect of the drug at a 564 molecular level at a given dose) and clinical effectiveness (how well the drug reduces the 565 severity or duration of the infection), as well as the tradeoff between side effects and 566 bioavailability. One drug might have moderate efficacy but have few side effects. Another 567 drug might have high efficacy, but have severe side effects at high doses that limit the 568 maximum tolerated dose or use of even moderate doses in prophylaxis. A drug might 569 also have a combination of beneficial and adverse effects (e.g., it might reduce viral 570 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint 30 replication early in infection, but also be immunosuppressive) [13,37]. Antiviral drugs like 571 Tamiflu retain their ability to block aspects of the viral life cycle (efficacy), but become 572 much less clinically effective as the time before treatment increases: (in adults Tamiflu is 573 most effective when given within 48 hours after exposure and thus is often used 574 prophylactically) [38]. 575 In this section we use our model to show the trade-offs between time-of-use and 576 potency of a drug that targets viral genome replication in a host cell. Several antiviral 577 medications for RNA viruses reduce the net viral replication rate by inhibiting synthesis of 578 viral RNA by the viral RNA polymerase. We focus on RNA-synthesis blockers in this paper 579 because viral genome synthesis exponentially increases the production rate of viruses 580 per cell, while the other stages of viral replication have linear amplification effects (see 581

Equations (6)-(9) in Models and methods). 582
To simulate the effects of treatment that targets RNA synthesis using different drug 583 efficacies and times of administration, we generated a series of simulations in which we 584 reduced , the replication rate of genomic material in the viral replication model 585 (Equation (7) in Models and methods), by different amounts and at different times in the 586 simulation. The "viral replication multiplier" represents the potency of the treatment, the 587 factor by which is reduced (either a low dose with high efficacy, or a high dose with 588 a less efficacy). The "time delay of application" is the simulation time at which is 589 reduced, which corresponds to the time after infection at which the treatment is 590 administered. To characterize therapeutic effectiveness, we distinguished three classes 591 of simulation outcomes: 592 593 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint Positive outcomes: effective treatment, where at least 50% of the epithelial cells 594 remain uninfected at the end of the simulation (green-shaded subplots). 595 Negative outcomes: ineffective treatment, where less than 10% of the epithelial 596 cells remain uninfected at the end of the simulation (orange-shaded subplots). 597 of the number of uninfected epithelial cells (Fig 8), virus-releasing epithelial cells (Fig 9), 604 the total amount of extracellular virus (Fig 10), the number of dead epithelial cells (Fig S7)  green-shaded subplots). If the drug is administered prophylactically or very soon after 633 infection (at 0 minutes) the treatment potency needs to be only 25% to achieve mostly 634 positive outcomes (effective treatment). Increasing the time to treatment increases the 635 potency required to achieve similar numbers of positive outcomes: the treatment is 636 effective for a potency of at least 37.5% if administered by 4000 minutes (67 hours, 2 ¾ 637 days), and at least 87.5% if administered by 6000 minutes (100 hours, 4 days). For all 638 potencies greater than 12.5%, early intervention prevents significant increase in the 639 number of virus-releasing cells (Fig 9,, and a small number of 640 immune cells suffices to stop the spread of infection (Fig S9,. In 641 this region, delaying treatment results both in a higher level of extracellular virus (Fig 10,642 green-shaded subplots) and more dead epithelial cells at the end of simulation (Fig 8,643 green-shaded subplots). With inhibited viral replication in the infected epithelial cells, the 644 extracellular virus decays until it is mostly cleared by the end of simulation (Fig 10). of simulation is relatively high (around 50%) and the treatment is usually moderately 689 effective (Fig 8). For potencies below 50%, the number of virus-releasing epithelial cells 690 remains approximately constant or continues to increase after treatment (Fig 9) and 691 significant levels of extracellular virus remain at the end of the simulation, and so in most 692 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint cases the treatment is ineffective (Fig 10). In this regime, variability between outcomes 693 for the same parameter values is higher than for potencies above 50%. 694 695 696 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint potency) produced simulation replicas that had instances of all three outcomes (Fig 11). 710 In a simulation replica with a positive outcome (Run 2, Fig 11A), the first uninfected 711 epithelial cell dies (as well as a few uninfected epithelial cells) before 4000 minutes (67  712 hours, 2 ¾ days), after which the total extracellular virus gradually decreases. At around 713 4000 minutes (67 hours, 2 ¾ days), an epithelial cell near the initially-infected cell 714 becomes infected, causing a recurrence of infection, whose spread was stopped by the 715 treatment. In contrast, simulation replicas with intermediate and negative outcomes (Runs 716 8 and 4, respectively, Fig 11A)

Model extensions 729
In this section we demonstrate the deployment of extensions to the framework 730 described in Models and methods, which we will refer to as the "main framework" when 731 discussing extensions in this and subsequent sections, as well as particularization of the 732 framework to specific biological problems like a different virus. We accomplish this parameter set). Implementation of randomly selected heterogeneous susceptibility is 755 available in the add-on modules library and is hosted on our repository for public use with 756 the module name "RandomSusceptibility". For each set of replicas we simulated ten 757 replicas using the RandomSusceptibility module while varying the fraction of 758 unsusceptible cells (like different locations in the lungs), from 10% to 50% unsusceptible 759 in intervals of 10% (Fig 12).   (Table S2). The integrated HCV model is implemented in the add-on 822 modules library using the aforementioned software architecture and is hosted on our 823 repository for public use with the module name "HCVIntegrated" (Fig 13A).  radially advanced outward from the initial site of infection ( Fig 13B). However, one notable 843 difference was that the prominent band of virus-releasing cells during spread of infection 844 only occurred in some simulation replicas (Fig 13B, replica "a"), while in other replicas 845 only small, isolated groups of cells became virus releasing (Fig 13B,replicas "b" and "c"). 846 Variability of outcomes was particularly notable among the ten simulation replicas. Some 847 simulation replicas produced widespread infection at a comparable timescale to that of 848 the baseline parameter set, between one and two weeks (Fig 13C). Such simulation 849 replicas were those that produced comparable spatial distributions of infected and virus-850 releasing cells (specifically, with a prominent band of virus-releasing cells) to those of the 851 main framework using the baseline parameter set. In other simulation replicas (e.g.,  the same potency to be more effective than later treatment, our model suggests that, for 968 antivirals, time of treatment is a more significant factor than potency in determining the 969 effectiveness of the therapy. Our model thus suggests that drugs that interfere with virus 970 replication are significantly more effective if used even at very low doses prophylactically 971 or very soon after infection, than they would be if used even at a high dose as a treatment 972 given later after exposure. Specifically, a prophylactic treatment in simulation which 973 reduces the rate of viral RNA synthesis by only 35% (35% potency) is more effective than 974 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint a treatment with 100% potency given two and a half days after infection, and has about 975 the same efficacy as a treatment with 50% potency given one day after infection. Our 976 model also showed that because of stochasticity in viral spread, later treatment at 977 moderate to high potency may still be effective in a subset of individuals. 978 Both parameter sweeps had regions with little variation in outcome between 979 replicas (e.g., the upper-right and lower-left corners of Fig 5). In regions of the parameter 980 space between these extremes (e.g., the unshaded areas in Figs 5 and 9), different 981 replicas showed dramatically different outcomes. One such parameter set in our drug 982 therapy simulations produced three distinct qualitative outcomes (i.e., positive, 983 intermediate and negative outcomes, Fig 11). For these parameters, replica outcomes 984 were particularly sensitive to stochasticity early in infection when only a few cells were 985 infected (Fig 11A), with delayed spread of infection from the first infected cell producing 986 more positive outcomes. Simulation replicas with negative outcomes (Fig 11A, Run 8) 987 had higher extracellular virus levels at earlier times than those with intermediate 988 outcomes (Fig 11A, Run 4), even though the fraction of each cell type was similar. Since 989 the viral replication module is deterministic, the primary cause of this difference is the 990 spatial distribution of cells. Spatial structure (e.g., infection of neighboring cells), 991 stochastic events (e.g., early cell death of infected cells before significant virus release) 992 and cell-to-cell variation (e.g., difference in viral release between cells) all affect the 993 variation between replicas. We are working to implement validated non-spatial models of viral infection and 1076 immune response as agent-based spatial models (e.g., viral production, cytokine 1077 secretion, tissue damage). By starting with a validated model that uses ordinary 1078 differential equations and adding spatial components gradually, we can calibrate our 1079 spatial models and validate our results. In ongoing work, we have developed a formal 1080 method for spatializing ODE models and employing their parameters such that these 1081 analogous spatial models reproduce the ODE results in the limit of large diffusion 1082 constants. Using this method, we can combine the ability to do rapid formal parameter 1083 identification of ODE-based models and to leverage published ODE model parameters 1084 with the flexibility of spatial modeling. In these cases, any differences between the ODE 1085 results and the spatial model can be definitively attributed to spatialization (e.g. The COVID-19 crisis has shown that drug discovery and therapy development both 1103 require new predictive capabilities that improve their effectiveness and efficiency. We 1104 have developed our framework to explore the relationship between molecular, cellular-1105 level and systemic mechanisms and outcomes of acute viral infections like SARS-CoV-1106 2, and to support development of optimal, patient-specific treatments to combat existing 1107 and new viruses. 1108 1109

1110
In this section we first present our model as a high-level conceptual model where 1111 we list each process included in an implementation-independent manner. We then detail 1112 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. H2O2 and nitric oxide. The mechanism of action of such agents varies but we assume that 1235 we can generalize such effects by modeling a single diffusive and decaying oxidizing 1236 agent field in the extracellular environment. The oxidizing agent is secreted by activated 1237 immune cells after persistent exposure to cytokine signals (I4→T3). We assume that the 1238 range of action of the oxidizing agent is short. Cell death is induced in uninfected, infected 1239 and virus-releasing epithelial cells when sufficiently exposed to the oxidizing agent 1240 (T3→E4). 1241 1242 I1 -Immune cell activation. Module I1 models immune cell maturation due to cytokine 1243 signaling. Immune cells mature at the recruitment site before being transported to the 1244 infection site as inactive immune cells (L1→Immune Cells). After infiltration, immune 1245 cells need to be exposed to local cytokine signals before activating (T2→I1). Once 1246 activated, immune cells chemotax along the cytokine field (I2) and amplify immune 1247 signaling by releasing cytokine molecules into the extracellular environment (I1→T2). 1248 Immune cells can also deactivate after a period of activation (I1 and Fig 2C). 1249 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a  ′. The final term, ℋ ℎ , models chemotaxis-directed cell motility, and is prescribed 1320 by module I2. The cell configuration evolves through asynchronous lattice-site copy 1321 attempts. A lattice-site copy attempt starts by random selection of a site in the lattice as 1322 a target, and a site ′ in its neighborhood as a source. A configuration update is then 1323 proposed in which the value ′ from the source site overwrites the value of in the target 1324 site. The change in total effective energy ℋ due to the copy attempt is calculated, and 1325 the update is executed with a probability given by a Boltzmann acceptance function, 1326 Here the intrinsic random motility ℋ * controls the stochasticity of accepted copy attempts. 1327 Updates that reduce the system's effective energy are always accepted. The unit of  ( Here ℎ is a Hill coefficient, is the cell's initial number of unbound receptors, is 1356 the virus-receptor association affinity, is the virus-receptor dissociation affinity, 1357 is a characteristic time constant of uptake and is the time represented by one MCS. At 1358 each simulation time step, the uptake probability is evaluated against a uniformly-1359 distributed random variable. When uptake occurs, the uptake rate is proportional to the 1360 local amount in the virus field (T1), and the probability of uptake is used to define the 1361 amount ( ) of virus taken up during the MCS, 1362 The amount absorbed by each cell ( ) is uniformly subtracted from the virus field 1363 over the cell's domain and the cell's number of cell surface receptors is reduced by the 1364 same amount. The amount of virus taken up ( ) is also passed to the cell's instance 1365 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. .

71
of the viral replication model (E2). Infected epithelial cells continue taking up viral particles 1366 from the environment until their cell surface receptors are depleted. 1367 1368 E2 -Viral replication. Our viral replication model combines equations and parameters 1369 from several sources to represent the replication of a generic virus [7,9,52,53]. The model 1370 contains four variables representing viral quantities in different states of the viral 1371 replication process: internalized virus (Equation (6), the process of unpacking), viral 1372 genome taking part in genomic replication (Equation (7), the process of viral genome 1373 replication), synthesized protein (Equation (8), the process of protein synthesis), and 1374 assembled and packaged virions (Equation (9), the process of assembly and 1375 packaging). Biologically, the only process which can exponentially increase the rate of 1376 virus production by a single cell is viral genome replication, so the equations include the 1377 positive feedback by in Equation (7). Biologically, factors like the cell's metabolism, 1378 limited number of ribosomes, maximum rate of endoplasmic reticulum function and 1379 activation of intracellular viral suppression pathways all limit production of viral 1380 components, so we include a Michaelis-type saturation term for the rate of replication in 1381 Equation (7). Each uninfected, infected and virus-releasing cell in the simulation contains 1382 an independent copy of the system of ordinary differential equations modeling the viral 1383 replication process, 1384 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint Here is the unpacking rate, is the viral replication rate, is the translating rate 1385 (rate at which viral genomes turn into RNA templates for protein production) and is the 1386 packaging rate.
is defined in E1 and is defined in E3. The saturation of 1387 the rate of viral genome replication is represented by a Michaelis-Menten function, where ℎ is the amount of at which the viral genome replication rate is reduced to 1389 2 (and the flux is reduced to  Table 1), with the additional complication that in our model, 1393 an epithelial cell does not release virus until reaches a threshold value. The timescale 1394 for tenfold increase of virus release when viral replication is maximal is 10 ≈ (10) (7.7 1395 hours for the reference parameter set in Table 1). The number of newly assembled virions 1396 is passed to the cell's instance of the viral release module (E3). See Fig 2B for   is governed by an ordinary differential equation of a dimensionless state variable that 1434 represents immune response due to local conditions and long-distance signaling. Our 1435 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020 75 convention is that when > 0, immune cells are recruited to the simulation domain; 1436 likewise, immune cells are removed from the simulation domain when < 0. Probability 1437 functions of describe the likelihood of immune cell seeding and removal, 1438 Pr(remove immune cell) = erf(− ) , < 0.
Here the coefficient is the sensitivity of immune cell addition and removal to the 1439 state variable . The dynamics of are cast such that, in a homeostatic condition, a typical 1440 number of immune cells can be found in the simulation domain, and production of cytokine 1441 (T2) results in additional recruitment via long-distance signaling (i.e., with some delay). 1442 We model this homeostatic feature using the feedback mechanism of the total number of 1443 immune cells in the simulation domain. Cytokine signaling is modeled as 1444 perturbing the homeostatic state using the term . Here is the total amount of 1445 decayed cytokine in the simulation domain and > 0 models signaling by transmission 1446 of cytokine to some far-away source of immune cells. We write the rate of change of as 1447 Here and control the number of immune cells in the simulation domain under 1448 homeostatic conditions. controls the delay between transmission of the cytokine to 1449 the lymph node and corresponding immune response by adjusting the rate of recruitment 1450 due to total cytokine (i.e., increasing increases the resulting delay). regulates 1451 the return of to an unperturbed state (i.e., = 0, increasing increases the rate of 1452 return to = 0). To determine the seeding location, the simulation space is randomly 1453 sampled times, and an immune cell is seeded at the unoccupied location with the 1454 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020 After ten hours, an activated immune cell becomes inactive, in which case evaluations of 1476 activation (Equation (15) 1498 I4 -Immune cell oxidizing agent cytotoxicity. Immune cells release a short-range, diffusive 1499 oxidizing agent when exposed to high cytokine concentration (T3). The oxidizing agent 1500 kills an epithelial cell of any type when the total amount of oxidizing agent in the domain 1501 of the cell ( ) exceeds a threshold for death ℎ , 1502 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020 Here is the diffusion constant of the extracellular virus and is the decay rate is 1507 the decay rate. Uptake and release by a cell at each location are determined using the 1508 viral internalization (E1) and the viral release (E3) modules, and are uniformly applied 1509 over all sites of the domain of the cell. 1510 1511 T2 -Cytokine transport. The change in concentration of the cytokine field is obtained 1512 by solving a reaction-diffusion equation of the following general form, 1513 The decay term represents cytokine leaving the simulation domain (e.g., in 1514 immune recruitment). To model immune signaling, the rate of cytokine secretion is 1515 described by an increasing Hill function of ( ( )) with Hill exponent ℎ = 2. The 1516 meaning of ( ( )) depends on the cell type and the Hill exponent can differ for other 1517 cell types and states. The rate of cytokine secretion is written as, 1518 Here ( ( ( ), )) is the maximum cytokine secretion rate for the cell type at , 1519 ( ( )) is a quantity that affects the cells cytokine secretion, ( ( ( ), )) is the 1520 cytokine uptake rate of the cell type at and ( ( ( ), )) is a dissociation coefficient 1521 of cytokine secretion for the cell type at .

1982
The number of uninfected epithelial cells for each simulation replica for each parameter set, plotted on a logarithmic 1983 scale, vs time displayed in minutes.

1989
The number of uninfected epithelial cells for each simulation replica for each parameter set, plotted on a logarithmic  Table 1).

2010
The number of uninfected epithelial cells for each simulation replica for each parameter set, plotted on a logarithmic 2011 scale, vs time displayed in minutes.

2012
2013 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint a proxy for the cytoplasmic plus-strand RNA molecules of the HCV model. Both 2029 modifications are described here. 2030 According to the HCV model, in the cytoplasm, 2031 internalization events, subgenomic replication within a particular cell occurs due to 2044 internalization of virus by the cell (6). As such, the final term of (S1) was added during 2045 integration such that internalized virus acts as a source for . 2046 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020.  Having selected and as the shared biological object of the two models, mass 2057 action → of the viral replication model of the main framework requires modification. 2058 We assume that decay of described in the HCV model leads to production of 2059 through intermediate processes. The viral replication model of the main framework ((6)-2060 (10)) then takes the modified form, 2061 = − , (S10) 2062 = , (S11) 2063 . CC-BY 4.0 International license (which was not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint this version posted September 26, 2020. . https://doi.org/10.1101/2020.04.27.064139 doi: bioRxiv preprint