A multiscale model of viral infection in epithelial tissues

18 Simulations of tissue-specific effects of primary acute viral infections like COVID19 19 are essential for understanding disease outcomes and optimizing therapies. Such 20 simulations need to support continuous updating in response to rapid advances in 21 understanding of infection mechanisms, and parallel development of components by 22 multiple groups. We present an open-source platform for multiscale spatiotemporal 23 simulation of an epithelial tissue, viral infection, cellular immune response and tissue 24 damage, specifically designed to be modular and extensible to support continuous 25 updating and parallel development. The base simulation of a simplified patch of epithelial 26 tissue and immune response exhibits distinct patterns of infection dynamics from 27 widespread infection, to recurrence, to clearance. Slower viral internalization and faster 28 immune-cell recruitment slow infection and promote containment. Because antiviral drugs 29 can have side effects and show reduced clinical effectiveness when given later during 30


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 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]. cells can die through bystander-effect mechanisms. In acute infections, adaptive immune 158 response leads to pathogen neutralization and clearance [48]. Viral loads usually 159 decrease rapidly as adaptive immune cells like CD8+ T-cells enter the tissue and 160 eliminate infected cells. Cells also begin to send out anti-inflammatory signals, shutting 161 down the immune response as viral clearance proceeds. Antigen presentation within the 162 immune system also induces activation of naive B-cell lymphocytes into antibody-163 producing memory B-cells and plasma cells, which leads to the production of antibodies. 164 The adaptive immune response remembers its exposure to previous pathogens and 165 8 provides the body with pathogen-specific effector cells and antibodies which neutralize 166 and clear them, providing long term immunity [49]. Tissue damage results from virus and 167 cytokine-induced cell death (which is first noticeable after 2 or 3 days) and from killing of 168 infected and uninfected cells by immune cells, which increases steadily until the end of 169 viral clearance (7-10 days). Tissue recovery and healing start around the time of viral 170 clearance and may last for several weeks. 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 238 cytotoxicity (I4). Activated immune cells participate in oxidative cytotoxicity (I4) and secrete oxidative agents into the 239 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

253
Immune cells are initially inactive and do not participate in the oxidative cytotoxicity (I4) or chemotax towards higher 254 concentrations of the cytokine field (I2). They become activated when they experience activation (I1). In all panels, 255 dashed arrows with barbed heads represent transformations, solid arrows with barbed heads represent transport and 256 solid arrows with lollipop heads represent regulation.

258
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). 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 (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)  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 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 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 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 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%. 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)  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 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. 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). ′. 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 ( 3 ) 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  Contact coefficients (all interfaces) 10 Selected comparably to [68] for low adhesion immune cell-immune cell and immune cellmedium interfaces 1 The diffusivity in water for a virus of radius 0.1 microns like SARS-CoV-2 according to Stokes-Einstein is about 3 microns 2 /s. The 1567 average steady-shear viscosity for lung mucus varies significantly and is shear thinning, but in the more viscous regions is found to 1568 vary for frequencies between 10 -4 and 102 Hz, spanning viscosity values as high as 103 Pa-s and as low as 10 -2 Pa-s. In general, at 1569 low shear rates, the viscosity of human mucus is as high as 104−106 times that of water [74]. Thus the minimal diffusion constant 1570 possible would be 3 x 10 -6 microns 2 /s and the maximal rate in water would be 3 microns 2 /s. 0.01 microns 2 /s is a reasonable geometric 1571 interpolation.