key: cord-0822626-lql6qoby
authors: Sego, T.J.; Mochan, Ericka D.; Bard Ermentrout, G.; Glazier, James A.
title: A Multiscale Multicellular Spatiotemporal Model of Local Influenza Infection and Immune Response
date: 2021-09-27
journal: J Theor Biol
DOI: 10.1016/j.jtbi.2021.110918
sha: 6ae0338dd6f8c7bd72e3adda2427a6f534f5fcb8
doc_id: 822626
cord_uid: lql6qoby

Respiratory viral infections pose a serious public health concern, from mild seasonal influenza to pandemics like those of SARS-CoV-2. Spatiotemporal dynamics of viral infection impact nearly all aspects of the progression of a viral infection, like the dependence of viral replication rates on the type of cell and pathogen, the strength of the immune response and localization of infection. Mathematical modeling is often used to describe respiratory viral infections and the immune response to them using ordinary differential equation (ODE) models. However, ODE models neglect spatially-resolved biophysical mechanisms like lesion shape and the details of viral transport, and so cannot model spatial effects of a viral infection and immune response. In this work, we develop a multiscale, multicellular spatiotemporal model of influenza infection and immune response by combining non-spatial ODE modeling and spatial, cell-based modeling. We employ cellularization, a recently developed method for generating spatial, cell-based, stochastic models from non-spatial ODE models, to generate much of our model from a calibrated ODE model that describes infection, death and recovery of susceptible cells and innate and adaptive responses during influenza infection, and develop models of cell migration and other mechanisms not explicitly described by the ODE model. We determine new model parameters to generate agreement between the spatial and original ODE models under certain conditions, where simulation replicas using our model serve as microconfigurations of the ODE model, and compare results between the models to investigate the nature of viral exposure and impact of heterogeneous infection on the time-evolution of the viral infection. We found that using spatially homogeneous initial exposure conditions consistently with those employed during calibration of the ODE model generates far less severe infection, and that local exposure to virus must be multiple orders of magnitude greater than a uniformly applied exposure to all available susceptible cells. This strongly suggests a prominent role of localization of exposure in influenza A infection. We propose that the particularities of the microenvironment to which a virus is introduced plays a dominant role in disease onset and progression, and that spatially resolved models like ours may be important to better understand and more reliably predict future health states based on susceptibility of potential lesion sites using spatially resolved patient data of the state of an infection. We can readily integrate the immune response components of our model into other modeling and simulation frameworks of viral infection dynamics that do detailed modeling of other mechanisms like viral internalization and intracellular viral replication dynamics, which are not explicitly represented in the ODE model. We can also combine our model with available experimental data and modeling of exposure scenarios and spatiotemporal aspects of mechanisms like mucociliary clearance that are only implicitly described by the ODE model, which would significantly improve the ability of our model to present spatially resolved predictions about the progression of influenza infection and immune response.

38 role in disease onset and progression, and that spatially resolved models like ours may be important 39 to better understand and more reliably predict future health states based on susceptibility of 40 potential lesion sites using spatially resolved patient data of the state of an infection. We can 41 readily integrate the immune response components of our model into other modeling and 42 simulation frameworks of viral infection dynamics that do detailed modeling of other mechanisms 43 like viral internalization and intracellular viral replication dynamics, which are not explicitly 44 represented in the ODE model. We can also combine our model with available experimental data 45 and modeling of exposure scenarios and spatiotemporal aspects of mechanisms like mucociliary 46 clearance that are only implicitly described by the ODE model, which would significantly improve 47 the ability of our model to present spatially resolved predictions about the progression of influenza 48 infection and immune response.

49 Keywords 50 influenza, immune response, multiscale modeling, cell-based modeling, cellularization 51 1 Introduction 52 Respiratory viral infections continue to be a serious public health concern, from mild seasonal 53 influenza strains to the highly pathogenic SARS-CoV-2 pandemic. In recent influenza strains 54 associated with highly pathogenic outcomes, excess inflammation and cytokine storm tend to be 55 major causes of mortality (1). Similarly, the recent COVID-19 epidemic has shown this 56 coronavirus induces a similar cytokine storm in many of the lethal infections (2). A deeper 57 understanding of the mechanisms involved in the initiation, proliferation, and reduction of the 58 inflammatory response is key to understanding the reasons why some infections can become lethal.

Spatiotemporal dynamics impact nearly all components of the resolution of a viral infection 60 both in vitro and in vivo. Viral replication rates, for example, depend on the type of cell the virus 61 has invaded, the family and strain of the virus, and the strength of the immune response deployed 62 against the pathogen. Viruses have been theorized to differ in replication rates between mucosal 63 and bronchial epithelial cells (3). In addition, some viruses have been shown to localize to certain 64 areas of the lungs rather than spread homogeneously throughout the respiratory tract. For instance, 65 the 2009 pandemic H1N1 strain has been shown to replicate more extensively throughout the lower 66 respiratory tract than either seasonal H1N1 or H5N1 (4). Seasonal H1N1 and H3N2 tend to 67 replicate primarily in the bronchi, while H5N1 replicates largely in alveoli (5). In addition to these 68 spatial differences, temporal differences in viral replication rates also play a part in differing levels 69 of pathogenicity between viral strains. Multiple experimental studies have shown that strains 70 exhibit distinct rates and mechanisms for cell entry, replication, and evasion of immune responses, 71 allowing certain strains to be more virulent than others (5-8). The immune response to viral 72 infection also includes many spatially-resolved biological processes, many of which are poorly 73 understood, such as the search strategies of CD4 + and CD8 + T cells leading to antigen recognition, 74 memory and effector T cell differentiation, migration via chemokinesis, chemotaxis and 75 haptotaxis, and cytotoxic killing of infected cells (9). Heterogeneous spread of virus and infected 76 cells has been theorized to affect the spread of infection through the lung; clusters of dead cells 77 near productively-infected cells may prevent the virus from spreading (10). This effect has been 78 seen after lethal H5N1 infection in ferrets (5); excessive damage in the lower respiratory tract 79 prevents the virus from spreading further through the lung and limits the peak of the viral load. 80 Thus, characterizing the spatial spread of the virus through the lung is critical to understanding the 81 intrahost immune response to the infection. 82 Mathematical modeling has long been used to explore various details of the immune 83 response to respiratory viral infections using ordinary differential equation (ODE) models (11-84 14). However, spatial effects cannot be explored in a typical ODE model, as these models are 85 founded on a well-mixed assumption that neglects spatially resolved biophysical mechanisms 86 (e.g., lesion shape). Spatial models of the immune response have been developed in recent years 87 to explore the effects of the spatial distribution of immune components on the resolution of 88 infection (10,15-21). However, to our knowledge, no spatial model exists that describes host-89 pathogen interactions during influenza infection with cellular resolution while considering detailed 90 descriptions of local and systemic aspects of both the innate and adaptive immune responses.

In this work, we combine the approaches of non-spatial ODE modeling and spatial, cell-92 based modeling to develop a multiscale, multicellular model of influenza infection and immune 93 response. We generate much of our spatial model from a calibrated ODE model that describes 94 infection, death and recovery of susceptible cells and innate and adaptive responses during 95 influenza infection (14) using cellularization (22), a recently developed method for generating 96 spatial, cell-based models from non-spatial ODE models. We develop models of cell migration 97 and other mechanisms not explicitly described by the ODE model, and determine new model 98 parameters to generate agreement between the spatial and original ODE models under certain 99 conditions. We compare results between the models to investigate the nature of viral exposure and 100 impact of heterogeneous infection on the time-evolution of the viral infection.

101 2 Models and Methods 102 The ODE model of in-host response to influenza A virus from which the spatial model is generated 103 describes infection and death of susceptible epithelial cells, and inflammatory, innate, adaptive 104 and humoral responses. Population dynamics consist of explicit expressions for uninfected, 105 infected and dead epithelial cells ( , and ), macrophages ( ), neutrophils in the blood and 106 infected tissue ( and ), antigen presenting cells (APCs, ), natural killer (NK) cells ( ), B cells 107 ( ), and CD4 + and CD8 + T cells ( and , respectively). Soluble signals of the model consist of 108 tumor necrosis factor (TNF, ), interleukins 10 and 12 (IL-10 and IL-12, and , respectively), 109 types I and II interferon (IFN, and , respectively), and generic chemokines ( ), antibodies ( ) 110 and reactive oxygen species (ROS, ). We employ the method of cellularization (22) to generate 111 a multiscale, multicellular, spatiotemporal model of local influenza A infection and immune 112 response in an epithelial sheet. For details of the complete ODE model and cellularized spatial 113 model, see Appendix 1 in Supplementary Materials.

Cellularization describes the relationships of measurements of quantity at various scales of 115 a biological system under well-mixed conditions. For a scalar quantity of a species at one scale, 116 a scalar quantity of the same species at another scale, and a field distribution of which = ( , ) 117 measures, according to cellularization,

119 where and are global and local scaling coefficients, respectively, and is the diffusion 120 coefficient of for diffusive species. For diffusive species, is the volume integral of over a 121 spatial domain, while for discrete objects of a particular type is the number of instances of the 122 type of object (e.g., the number of neutrophils). Reaction-diffusion equations for locally 123 heterogeneous soluble signals are generated from non-spatial descriptions. For a rate equation 124

for chemical species and and number of a cell type , the analogous = ( , ) + ( , ) 125 reaction-diffusion equation for is is written as a function of time to support changes in type (e.g., cell death). ( , )

Cellularization formulates cell-based stochastic dynamics using the Poisson cumulative 133 distribution function from reaction kinetics that describe the inflow, outflow, and transitions by 134 type (e.g., from alive to dead) of cell populations. For a number of cells of cell type with mean 135 inflow rate , mean outflow rate , and mean transition rate to cell type (i.e., 136 and for cells of type ) over a period , 146 Here a process for cell with total available contact surface area is mediated by contact with a 147 cell of type with total available contact surface area , and is the contact area of cell , 148 with a type cell. (14) to comparable work on multiscale, spatial, cell-based 152 modeling of viral infection and immune response (21) using cellularization. We consider a quasi-153 two-dimensional spatial domain in which local infection occurs in a fixed planar sheet of epithelial 154 cells. Recruitment of various immune cell populations is governed by organismal-level dynamics 155 coupled with signaling from within the spatial domain, where motile, locally acting immune cells 156 are recruited from outside the spatial domain and placed on top of the epithelial sheet. Likewise 157 organismal-level soluble signals are coupled with locally heterogeneous distributions of diffusive 158 species in the spatial domain. We refer to model objects whose dynamics are explicitly modeled 159 in the spatial domain as local, and likewise to those modeled with an ODE as global.

To model the spatial effects of infection, we model extracellular virus and uninfected, 161 infected and dead cells as local (Figure 1 ). Type I IFN is modeled as local to model the spatial 162 effects of antiviral resistance. Macrophages, chemokines and IL-10 are modeled as local to model 163 the spatial effects of macrophage migration and diffusion in local inflammatory signaling. 164 Likewise, NK and CD8 + T cells are modeled as local to model the effects of contact-mediated 165 killing of infected cells in the immune response. B cells and blood neutrophils are not present at 166 the local site of infection, and so are modeled as global. APCs primarily recruit other immune cell 167 types according to the ODE model, and so we neglect the spatial aspects of their type I IFN release 168 and model them as global. It follows that type II IFN, IL-12 and CD4 + T cells are modeled as 169 global, since they immediately affect global objects. Antibodies, ROS and TNF could reasonably 170 be modeled as local, however we assume that their rates of diffusion are sufficiently fast to 171 approximate them as uniform throughout the simulation domain and model them as global. It 172 follows that tissue neutrophils are modeled as global, since they release global ROS. other model objects are treated as homogeneously acting due to their absence in the spatial domain (e.g., blood neutrophils) or 177 their spatial properties (e.g., highly diffusive ROS) and are shown in solid borders. Analogous spatial and cell-based models of 178 processes within, and across the boundary of, the spatial domain are derived from the ODE model using cellularization.

To model migration of local immune cells, we use the Cellular Potts model (CPM, or 180 Glazier-Graner-Hogeweg model). The CPM is a lattice-based hybrid kinetic Monte Carlo method 181 that represents generalized cells and medium as discrete, deformable, volume-excluding objects 182 (23). Cell motility in the CPM is modeled as the stochastic exchanging of lattice sites by cells and 183 medium according to minimization of a system effective energy , in this work written as, ℋ

.

The first term implements a volume constraint in each cell by cell type, the second term models 186 adhesion at intercellular and cell-medium interfaces by cell type according to contact coefficients 187 using a neighborhood of each site, and the third term models logarithmic chemotaxis by ( ) 188 cell type and field distribution according to a chemotaxis parameter , local field concentration 189 and center-of-mass measurement of . In this work we use a second-order Manhattan 190 neighborhood for adhesion calculations, while applications of adhesion and chemotaxis modeling 191 are described in the following section. In general, the CPM randomly selects a pair of neighboring 192 lattice sites and considers whether the identification at one of the sites copies itself to the other 193 site, called a copy attempt, which occurs with a probability according to a Boltzmann acceptance 194 function, 195

Here denotes the copy attempt where the identification at copies to , is

′ ℋ * 197 the intrinsic random motility that affects the stochasticity of copy attempts, and is the change ∆ℋ 198 in due to the copy. One simulation step, called a Monte Carlo step (MCS), consists of having ℋ 199 considered a number of copy attempts equal to the total number of lattice sites. 206 Here is a model parameter of the ODE model.

Diffusive transport is assumed to occur in a homogeneous medium, where for extracellular 208 virus, chemokines, IL-10 and type I IFN we define the diffusion coefficients , , and , 209 respectively. The spatial model describes spatial and cell-based analogues of all mechanisms 210 described by the ODE model for each heterogeneous species using partial differential equations 211 (PDEs) of diffusive transport. Diffusive transport modeling of the extracellular virus includes 212 general decay, decay by the action of antibodies, mucociliary clearance, uptake by uninfected cells 213 and release by infected cells; of chemokines includes general decay, and release by macrophages 214 regulated by TNF and the presence of dead cells; of IL-10 includes general decay, release by 215 macrophages regulated by TNF and the presence of dead cells, and release by uninfected cells; of 216 type I IFN includes general decay, release by APCs, and release and uptake by infected cells 217 (Table 1) . 218 219 Table 1 . PDEs generated from cellularization of the influenza ODE model for virus , type I IFN , chemokines and IL-10 220 according to a general form for reaction-diffusion transport with diffusion coefficient. All symbols with subscripts are parameters 221 from the ODE model. and are the global and local scaling coefficients, respectively, according to cellularization. 223 type (e.g., uninfected , infected , macrophage ). , , and are the total antibodies, dead cells, APCs and TNF in ( , ) * 224 the spatial domain. is a mean cellular measurement of .

Decay rate Source rate

The ODE model employs the Allee effect with a critical population of uninfected cells, 227 above which recovery of uninfected cells occurs, and below which additional death of uninfected 228 cells occur. We cellularize this mechanism by splitting it into individual stochastic events, of which 229 a mean rate of death due to the Allee effect is considered for each uninfected cell, and a ( ) 230 mean rate of recovering a dead cell is considered. The process of cell recovery is ( ) 231 implemented as the transitioning of a dead cell to an uninfected cells (21). By treating both 232 mechanisms of the Allee effect as contact-mediated and applying the well-mixed conditions (see 233 Appendix 1 in Supplementary Materials for derivations), both rates can be described in terms of 234 local conditions of each uninfected and dead cell, Table 1 . , and are the total chemokines, APCs and 253 ROS in the spatial domain, and is the resistance of cell .

Transition rate

256 Beyond the cell-based models that can be generated from the ODE model using cellularization, 257 the ODE model implicitly describes spatiotemporal aspects of influenza A infection and immune 258 response that we can infer, impose or propose using additional data, assumptions and hypotheses. 259 For the simplest case, employing the CPM requires imposing a volume constraint on each cell, the 260 quantities, but not geometries, of which the ODE model describes. As such, we impose an 261 approximate cell diameter of 10 µm on all cells according to the typical size of epithelial cells and 262 simplification of negligible differences in typical volume among cell types (Table 3) .

263 The ODE model describes the killing of infected cells proportionally to the number of NK 266 and CD8 + T cells. In the spatial model, we place macrophages and NK and CD8 + T cells at the site 267 of infection and explicitly model their shape and motility, which provides the opportunity to 268 generate an explicit description of the spatiotemporal mechanisms involved in local immune cells 269 locating and eliminating infection. We 

Based on phagocytosis and inflammatory signaling by macrophages, we model 285 macrophages as chemotaxing up gradients of extracellular virus, and NK and CD8 + T cells as 286 chemotaxing up gradients of chemokines. We also model the specialization of CD8 + T cells as 287 their chemotactic sensitivity being twice that of NK cells. Based on contact-mediated killing of 288 infected cells, we model stronger adhesion of immune cells to infected cells compared to 289 uninfected and dead cells. We determined in early computational experiments that generating an 290 effective local immune response also required preferential attachment that prevents excessive 291 homotypic aggregation of immune cells but allows both heterotypic aggregates of immune cells 292 and dispersion of immune cell aggregates according to chemoattractant distributions. As such, we 293 model adhesion of immune cells to other immune cell types and the medium the same as to infected 294 cells, and to immune cells of the same type the same as to uninfected and dead cells. 295 We approximated the diffusive characteristics of local soluble signals by diffusion length 296 (i.e., for diffusion coefficient and decay rate ) in units of cell diameters, using the decay √( / ) 297 parameters of the ODE model. Diffusion of extracellular virus and chemokines were approximated 298 with diffusion lengths of five and ten cell diameters, respectively, based on previous, comparable 299 modeling work on local infection and immune response (21). The diffusion length of IL-10 was 300 assumed to be the same as that of chemokines, while type I IFN was modeled with a diffusion 301 length of two cell diameters to model local anti-viral signaling. 302 2.4 Implementation Details 303 Simulations were performed with comparable configurations to those in similar modeling work on 304 local infection and immune response (21). All simulations were executed in CompuCell3D (24) 305 with either a lattice planar dimension of 0.3 mm or 1.0 mm. Every lattice was discretized with a 306 discretization length of 2 µm for cells that, on average, occupied 25 sites (Table 4 ). The local 307 scaling coefficient 4×10 -8 µm -2 was calculated from the total number of epithelial cells = 308 according to the ODE model (250k) and cell volume constraint . The local scaling coefficient 309 was calculated as the ratio of the number of epithelial cells in the simulation domain to those in 310 the ODE model parameters, and was 0.0049 and 0.04 for lattices with planar dimensions of 0.3 311 mm and 1.0 mm, respectively. Neumann and periodic conditions were applied to boundaries 312 parallel and orthogonal, respectively, to the epithelial sheet. Epithelial cells were arranged in a 313 uniform grid of 5x5 squares. All simulations used a time step of one minute per step, which was 314 determined to be sufficiently small for numerically stability, particularly of type I IFN signaling. 315 All ODE model parameters were taken from (14). Scaling was performed by epithelial cell 316 population for an ODE model epithelial cell population of 250k cells. Using the prescribed volume 317 constraint, simulations using only the ODE model showed the potential for local immune cells to 318 exceed the available space in the immune cell layer, an event called overcrowding in 319 cellularization. To mitigate overcrowding, we employed the cellularization strategy of partially 320 homogenizing local immune cell populations, where all macrophages are explicitly represented in 321 the spatial domain (since they provide directional signaling), while 25% of NK and CD8 + T cells 322 act homogeneously as scalar-valued functions. All local immune cells were seeded into the 323 immune cell layer with a seeding fraction of 1% according to field values of their chemoattactrants 324 (i.e., by virus for macrophages, by chemokines for NK and CD8 + T cells 369 3.2 Disagreement in Large Epithelial Patches 370 Having shown acceptable agreement between the spatial and ODE models under marginally 371 stochastic initial infection conditions, we generated a spatial equivalent of the lethal scenario to 372 which the ODE model was calibrated by exposing large, uninfected epithelial patches to an initial 373 viral load. We simulated 50 replicas of 1.0 mm x 1.0 mm epithelial patches, which, for the chosen 374 spatial model parameters and 250k epithelial cells in the original ODE model simulations, 375 collectively total two model organisms of the ODE model.

In all simulation replicas, at most, marginal infection occurred, in strong disagreement with 377 ODE model predictions of a lethal outcome at around ten days (Figure 4) . Figure 4A shows spatial 378 distributions of one representative replica where any notable infection occurred, which consisted 379 of an isolated lesion of infected and dead cells that recovered within two weeks. During early 380 progression of such lesions, inflammatory signaling recruited significant numbers of macrophages, 381 which localized at the lesion and subsequently recruited local NK and CD8 + T cells. Some new, 382 later infection sites were also observed but well mitigated by present antibodies and quickly 383 eliminated by the already stimulated immune response. In these cases, present local immune cells 384 migrated with the general pattern of macrophages chemotaxing towards infected cells, followed 385 by present, and reinforced by newly recruited, local NK and CD8 + T cells. In many other 386 simulation replicas, no infection occurred, and the initial viral load decayed with no indication in 387 the epithelial patch of exposure to virus ( Figure 4B ). 

395 Because the initial viral load in the calibrated lethal scenario of the ODE model did not generate a 396 lethal outcome in the spatial model, we tested varying initial viral loads to determine at what order 397 of magnitude of initial viral load the spatial model generates a lethal outcome. Since the ODE 398 model was calibrated to both non-lethal and lethal scenarios, where the lethal scenario differed 399 from the non-lethal scenario only by a 10-fold increase in initial viral load, we performed a 400 logarithmic parameter sweep of initial viral load by beginning with the spatial model equivalent 401 to the non-lethal scenario, and increasing the initial viral load by a factor of 10 until the spatial 402 model produced lethal outcomes in twenty simulation replicas. 

We found that the spatial model begins to generate lethal outcomes when the initial viral 408 load is at least greater than the initial viral load of the calibrated lethal scenario by a factor of 100 409 ( Figure 5 ). Increasing the initial viral load of the lethal scenario ( Figure 5 , initial viral load of 10) 410 by factors of 10 and 100 did not produce a lethal outcome in any simulation replica over two weeks 411 of simulation time, though, besides the difference in outcome, the latter produced comparable 412 predictions to those of the ODE model. A 1k increase in initial viral load from the lethal scenario 413 produced at least nearly lethal outcomes in all simulation replicas, with the number of uninfected 414 cells reaching minima at least as low as 10 cells (i.e., 0.1% uninfected). Many simulation replicas 415 produced no uninfected cells at times as early as three and a half days, compared to about two and 416 a half days according to the ODE model (i.e., when the number of uninfected cells is less than 1 417 according to the ODE model). In some replicas, marginal numbers of uninfected cells persisted as 418 late as twelve and a half days, while two replicas demonstrated recovery of the epithelial patch and 419 a corresponding non-lethal outcome. For this initial viral load, spatial model results disagreed 420 otherwise only in amount of extracellular virus for replicas that produced a non-lethal outcome. 421 We found these differences to be due to the difference in treatment of cell populations (i.e., as 422 continuous quantities in the ODE model and as discrete quantities in the spatial model), where cell 423 populations being less than one exhibits no notable effects in the ODE model, whereas discrete 424 cell populations in the spatial model cease to exhibit any effects by having a number of cells equal 425 to zero. For the calibrated non-lethal scenario ( Figure 5 , initial viral load of 1), no infection 426 occurred in 19 out of 20 replicas, and in the one replica that did experience any infection, the 427 maximum number of infected was two orders less, and occurred about two days earlier, than that 428 of the ODE model.

For simulation replicas that did experience significant amounts of infection (e.g., for those 430 with lethal outcomes), we observed multiple sites of significant infection within the first day after 431 exposure ( Figure 6 ). These sites were locations of significant recruitment of local macrophages 432 and subsequent recruitment of local NK and CD8 + T cells, as well as localized type I IFN, which 433 later became more homogeneous due to production by nonlocal APCs. Spatial distributions of 434 chemokines and IL-10 showed gradients most apparently at around two days of simulation time, 435 with IL-10 being greater in regions with significant accumulation of local immune cells, and 436 became mostly homogeneous by around one week of simulation time when immune cells mostly 437 covered the epithelial patch. For simulation replicas with the highest initial viral load that 438 recovered (e.g., as in Figure Since much higher initial viral loads were required to generate significant infection in large 449 epithelial sheets using the spatial model compared to the ODE model, we then tested agreement 450 between the ODE and spatial models while varying initial infection fraction in 1 mm x 1 mm 451 epithelial patches. We varied the initial infection fraction in a logarithmic sweep at intervals of 452 0.1%, 0.5%, 1% and 5% and simulated twenty simulation replicas for each initial infection 453 fraction. As in Section 3.1, 5% initial infection fraction can produce fatal outcomes but in large 454 epithelial sheets, though at times no earlier than about five days, but can also produce non-fatal 455 outcomes (Figure 7) . For all simulation replicas subjected to 5% initial infection fraction, the 456 epithelial patch experienced infection comparably to that predicted by the ODE model, however 457 in some cases a few uninfected cells survived and initiated recovery. As initial infection fraction 458 decreased, peak extracellular virus and infected cells in the spatial model occurred later and with 459 lesser magnitude, and all simulation replicas produced a non-fatal outcome for initial infection 460 fraction less than or equal to 1%, all of which are fatal according to the ODE model. 465 466 510 spatial model epithelial patch. The cellularized Allee effect, which was recast to make associated 511 death and recovery mechanisms dependent on the state and local conditions of individual cells, 512 also produced differences in ODE and spatial model results by allowing recovery with very few 513 total uninfected cells in the spatial model. While we found differences in associated cell deaths to 514 be marginal (and hence, not shown), the spatial model can produce recovery of the epithelial patch 515 in scenarios where associated cell death and a corresponding fatal outcome occur in the ODE 516 model (e.g., Figure 5 ), depending on the state and local conditions of uninfected epithelial cells.

Differences in emergent dynamics due to considering local conditions highlight a 518 potentially important role of cell-based, spatial models to do model validation when developing 519 mathematical models to describe observations of multicellular systems. In previous work, 520 cellularization of simpler viral infection models has shown that spatial mechanisms like diffusion-521 limited transport of extracellular virus can result in well fitted, but significantly reduced, ODE 522 model parameters describing the infectiousness of a virus (22). In the present work, spatial 523 particularities of initial conditions demonstrated results of a comparable category in the sense that 524 initial exposure could only generate comparable dynamics between the spatial and ODE models 525 for initial exposures that were several orders of magnitude greater than previously considered when 526 using on the ODE model. In both the previous case of diffusion-limited viral transport and present 527 case of localized initial exposure, a cellularized model demonstrated significant biases of a 528 particular model structure (here, continuous non-spatial). Such biases could lead to misguided 529 parameter selection or even inappropriate acceptance of a model. In this sense, not only do cell-530 based, spatial models generated using cellularization like the one developed in the present work 531 present new opportunities for biological inquiry using mathematical modeling, they also present 532 opportunities for validation of existing and future non-spatial models by subjecting non-spatial 533 models to spatial conditions and effects. As shown in this work, such activities like considering 534 the spatial details of initial conditions can provide new insights into both a non-spatial model, and 535 the biological system for which the model was developed to describe. 

An approach to effectively modeling local features of exposure would significantly 550 improve the ability of our cellularized model to present spatially resolved predictions about the 551 progression of influenza infection and immune response, though will likely require considering 552 mechanisms that are only very implicitly described by the ODE model such as mucociliary 553 clearance. As such, future work should combine the cellularized model presented here with 554 available experimental data and modeling of exposure scenarios. Likewise, future work should 555 further explore and develop a cellular basis for the mechanisms represented by the Allee effect, in 556 particular what all is represented when imposing an organismal-level property like the number of 557 a particular cell type onto the fate of individual cells (e.g., levels of growth factors, hormones, 558 blood pressure). Such work would be particularly relevant for cellularized ODE models with fewer 559 mechanisms of death than the ODE model of this work, in which case the role of the Allee effect 560 would be clearer, and likely more significant. Furthermore, it may be possible to identify some 561 significant sources of differences to due discrete, stochastic dynamics independently of spatial 562 effects (as suggested by a reviewer of this manuscript) by employing modified forms of the ODE 563 model that impose important discontinuities onto the otherwise continuous ODE model variables 564 (e.g., a minimum population of infected cells to produce virus). Such investigations could inform 565 about which aspects of a spatial model are necessary to interrogate the mechanisms responsible 566 for observed host-pathogen interactions, and about which aspects of a non-spatial ODE model are 567 sensitive to small quantities.

The current type I IFN model constitutes the overwhelming majority of computational cost 569 of the spatial model. In particular, calculating a cellular property like resistance from the mean 570 value of a local diffusion field requires sufficiently small time steps, since its downstream effects 571 include regulation of future type I IFN production. We plan to make improvements to cellularized 572 mechanisms associated with type I IFN production, as well as to CompuCell3D, to permit larger 573 time steps and better facilitate computational performance. Such improvements are particularly 574 critical to modeling bigger tissue patches and more complicated tissue geometries, and simulating 575 longer scenarios.

Lastly, we are particularly interested in further developing other spatial aspects of the 577 cellularized model to further elucidate associated cellular and spatial aspects of influenza infection 578 and immune response, as well as the immune response, in general. We can use in vitro data of cell 579 migration to refine the proposed locomotion model of local immune cells at a site of infection, 580 both to better understand how adhesion and chemotaxis affect the effectiveness of the immune 581 response, and to isolate necessary additional model mechanisms to better represent local aspects 582 of immune response. At a broader level, we can perform similar activities to this project but for 583 other sites of interest associated with the ODE model. For example, the ODE model presents 584 systemic response data that can guide development of spatiotemporal, cell-based models of B cell 585 maturation, antigen presentation, and antibody production and circulation. We envision a 586 computational framework consisting of multiple compartments simulating spatiotemporal models 587 of various sites of interest throughout an organism (e.g., multiple sites of infection, lymph nodes, 588 thymus), which could be interconnected using the techniques of cellularization in similar fashion 589 to what was employed in this work.

Our simulation scenario (i.e., the periodic boundary conditions) implies that our simulation 470 replicas are constituent elements of a patterned system, the collection of which the ODE model infection and immune response under the same boundary and initial infection conditions

This is particularly important in therapeutics and modeling, alike, in that spatially 484 resolved patient data of the state of infection may elucidate future health states based on 485 susceptibility of potential lesion sites, which could be better understood and more reliably 486 predicted with spatially-resolved models of the type presented in this work. 487 Differential adhesion and chemotaxis parameters of the introduced spatial model 488 mechanisms were formulated qualitatively, and only roughly calibrated to recapitulate ODE model 489 results in Section 3.1. Interestingly, the employed differential adhesion was necessary to 490 recapitulate ODE model results, the role of which is currently not well-defined. These roles are 491 fairly intuitive when considering the observed temporary aggregation of NK and CD8 + T cells in 492 simulations. Were the differential adhesion scheme employed such that NK and CD8 + T cells show 493 a preferential attachment to each other, then the observed aggregations at sites of infection due to 494 recruitment would result in ineffective subsequent elimination of infected cells due to the 495 continued aggregation of NK and CD8 + T cells. As such, the model predicts that aggregation of 496 local cytotoxic immune cells is due to chemotactic signaling and preferential attachment to 497 infected cells

The most prominent differences between the spatial and ODE models all resolve to 501 localization of type I IFN and recovery. The ODE model, and correspondingly the cellularized 502 spatial model, describe saturated release of type I IFN, the saturation of which is diffusion-limited 503 in the spatial model. This leads to differences not only in total over-production of type I IFN in the 504 spatial models, but also in downstream over-production of virus (i.e., diffusion-limited anti-viral 505 resistance)

However, such differences between the ODE and spatial models were shown to be marginal under 507 certain exposure conditions (e.g., very high initial viral load or infection fraction) and, as 508 previously described, to be significant when the lack of representing localization of virus in the 509 spatial model significantly inhibits the progression of infection

We used the developed spatial model to show how exposure to virus 597 should be locally concentrated to generate significant infection in an epithelium, while uniform 598 exposure to virus is likely ineffective. We also used our developed spatial models to elucidate the 599 necessary roles of differential adhesion and local chemotaxis for an effective local immune 600 response, both concerning the locating by macrophages, and the eliminating by NK and CD8 + T 601 cells, of infected cells. 602 Funding 603 TJS and JAG acknowledge funding from National Institutes of Health grants U24 EB028887 and 604 R01 GM122424 and National Science Foundation grant NSF 1720625. GBE acknowledges 605 funding from National Science Foundation grant NSF 1951099. This research was supported in 606 part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive 607 Technology Institute

Declaration of Interest 610 JAG is the owner/operator of Virtual Tissues for Health, LLC, which develops applications of 611 multiscale tissue models in medical applications, and is a shareholder in Gilead Life Sciences

Appendix 2: Complete list of all behaviors, roles and properties of the cell types and fields of the 615 model

Source Code: Simulation package to run the model in CompuCell3D

Fatal outcome 619 of human influenza A (H5N1) is associated with high viral load and hypercytokinemia. 620 Nature Medicine

Longitudinal analyses 622 reveal immunological misfiring in severe COVID-19

Higher Level of 624 Replication Efficiency of 2009 (H1N1) Pandemic Influenza Virus than Those of Seasonal 625 and Avian Strains: Kinetics from Epithelial Cell Culture and Computational Modeling

Severity of pneumonia due to new H1N1 influenza virus in ferrets is intermediate 629 between that due to seasonal H1N1 virus and highly pathogenic avian influenza H5N1 virus. 630 J Infect Dis

Comparison of temporal and spatial dynamics of seasonal H3N2, pandemic 633 H1N1 and highly pathogenic avian influenza H5N1 virus infections in ferrets

Pathogenic Avian Influenza H5N1 Viruses Elicit an Attenuated Type I Interferon Response 637 in Polarized Human Bronchial Epithelial Cells

An 639 Ultrasensitive Mechanism Regulates Influenza Virus-Induced Inflammation

Avian 643 influenza viruses that cause highly virulent infections in humans exhibit distinct replicative 644 properties in contrast to human H1N1 viruses

T cell migration, search strategies and mechanisms. 647

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

Kinetics of Influenza A 651 Virus Infection in Humans

A dynamical model of human immune response to 653 influenza A virus infection

Modeling Within-655 Host Dynamics of Influenza Virus Infection Including Immune Responses

The 659 inflammatory response to influenza A virus (H1N1): an experimental and mathematical 660 study

Modelling dynamics of the type I 662 interferon response to in vitro viral infection

Linking Influenza 664 Virus Tissue Tropism to Population-Level Reproductive Fitness

A simple cellular automaton model for influenza A 668 viral infections

Hybrid approach to model the spatial 670 regulation of T cell responses

How is the effectiveness of immune surveillance impacted by the 673 spatial distribution of spreading infections?

A computational 676 model tracks whole-lung Mycobacterium tuberculosis infection and predicts factors that 677 inhibit dissemination

A 680 modular framework for multiscale, multicellular, spatiotemporal modeling of acute primary 681 viral infection and immune response in epithelial tissues and its application to drug therapy 682 timing and effectiveness

Generation of multicellular 684 spatiotemporal models of population dynamics from ordinary differential equations, with 685 applications in viral infection

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

Multi-Scale 689 Modeling of Tissues Using CompuCell3D

691 Comparing the Ex Vivo Fitness of CCR5-Tropic Human Immunodeficiency Virus Type 1 692 Isolates of Subtypes B and C

A 694 Dual Infection/Competition Assay Shows a Correlation between Ex Vivo Human 695

Immunodeficiency Virus Type 1 Fitness and Disease Progression

BioNetGen: software for rule-based 698 modeling of signal transduction based on the interactions of molecular domains