key: cord-0857535-whs3uret
authors: Rowlatt, C. F.; Chaplain, M. A. J.; Hughes, D. J.; Gillespie, S. H.; Dockrell, D. H.; Johannessen, I.; Bowness, R.
title: Modelling the within-host spread of SARS-CoV-2 infection, and the subsequent immune response, using a hybrid, multiscale, individual-based model. Part I: Macrophages
date: 2022-05-06
journal: bioRxiv
DOI: 10.1101/2022.05.06.490883
sha: f2e987c229ecbfcc7c88b75d13f10bb1cea110bc
doc_id: 857535
cord_uid: whs3uret

Individual responses to SARS-CoV-2 infection vary significantly, ranging from mild courses of infection that do not require hospitalisation to the development of disease which not only requires hospitalisation but can be fatal. Whilst many immunological studies have revealed fundamental insights into SARS-CoV-2 infection and COVID-19, mathematical and computational modelling can offer an additional perspective and enhance understanding. The majority of mathematical models for the within-host spread of SARS-CoV-2 infection are ordinary differential equations, which neglect spatial variation. In this article, we present a hybrid, multiscale, individual-based model to study the within-host spread of SARS-CoV-2 infection. The model incorporates epithelial cells (each containing a dynamical model for viral entry and replication), macrophages and a subset of cytokines. We investigate the role of increasing initial viral deposition, increasing delay in type I interferon secretion from epithelial cells (as well as the magnitude of secretion), increasing macrophage virus internalisation rate and macrophage activation, on the spread of infection.

Individual responses to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, the causative agent of the COVID-19 pandemic, are heterogeneous; for the original (Wuhan) strain, and subsequent Alpha to Delta variants, the majority of cases present a mild course of infection (with a significant proportion being asymptomatic) and do not require hospitalisation, whilst some infections develop into coronavirus disease (COVID- 19) , which not only require hospitalisation but can be fatal.

A dysregulated immune response is associated with severe COVID-19, characterised by an aberrant immune cell distribution, as well as elevated inflammatory profiles, [1, 2, 3, 4, 5] . Indeed, early studies demonstrated a loss of resident alveolar macrophages, an increase in inflammatory monocyte-derived macrophages and neutrophils, loss of regulatory T cell responses, and an overall decrease in lymphocytes, in patients with severe disease, as well as an increase in proinflammatory cytokines (such as interleukin 1 (IL-1), interleukin 6 (IL-6) and tumor necrosis factor (TNF-α)) [6, 1, 7, 8, 9] . Recent studies have

The presented model extends, and adapts, an established model for the within-host spread of Mycobacterium tuberculosis (M. Tb) infection [33] to simulate the spread of SARS-CoV-2 infection. Section 2.1 introduces the model design, whilst Sections 2.2-2.6 introduce the model components. The model supports a uniform monolayer of lung epithelial cells (see §2.5); extracellular virus diffusion (see §2.2), as well as virus-specific entry and replication pathways (see §2.4); subsequent innate immune response, such as macrophages, (see §2.6); and cytokine signalling from epithelial and immune cells (see §2.3). A schematic is given in Fig. 1 (a).

Let the time interval be denoted by I = (0, T ] ⊂ R such that the closure is denoted by I = [0, T ], where T ∈ R denotes the simulation time. Let Ω ⊂ R 2 denote a two-dimensional domain containing a uniform grid of lung epithelial cells and randomly distributed cross-sections of blood vessels (see Fig. 1(b) ), which is a reasonable assumption provided that the blood vessels are perpendicular to, and that there are no branching points through, the plane of interest (see Bowness et al. [33] and the references therein). We ignore any temporal dynamics or spatial changes of these vessels. Let Ω e ⊆ Ω and Ω b ⊂ Ω denote the sets of lung epithelial cells and blood vessel cross-sections, respectively, such that

where Ω e k ⊂ R 2 , k = 1, . . . , N epi , represents a single lung epithelial cell; Ω b k ⊂ R 2 , k = 1, . . . , N bves , represents the cross-section of a blood vessel; and N epi , N bves denote the number of epithelial cells and blood vessel cross-sections, respectively. Clearly, we have Ω = Ω e ∪ Ω b and N grid = N epi + N bves , where N grid denotes the size of the grid (see Fig. 1(b) ). Furthermore, we define Ω T = Ω × I, and ∂Ω T = ∂Ω × I, where ∂Ω denotes the boundary of Ω.

Initially, the monolayer of epithelial cells is randomly seeded with a user-defined number of individual virus particles, known as virions (see §2.2), which can be provided via a multiplicity of infection (MOI). Furthermore, the grid is initially populated with a user-defined number of tissue-resident immune cells (specifically, macrophages), which are assumed to be in a resting state (see §2. 6 ). 

Due to the high viral loads which can be observed [24, 21, 40] , modelling extracellular virus particles (virions) as individual discrete agents is computationally challenging. Therefore, we choose to model the extracellular virus as a continuous field using a reaction-diffusion equation. Let V : Ω T → R denote the number concentration of extracellular virus (units: virions/mm 2 ). Then the spatiotemporal evolution of the extracellular virus is given by

V (x, 0) = V 0 (x), x ∈ Ω, (2.2.1c) where D v denotes the spatially-dependent viral diffusion coefficient; s v denotes the source of virus from exocytosed virions following replication; u v denotes the uptake of virus following engagement with ACE2 receptors; d v denotes the decay of extracellular virus; n Ω denotes the outward unit normal to Ω from ∂Ω; and V 0 : Ω → R denotes the spatially-dependent initial condition, which is determined from either a user-defined number of virions or a user-defined initial multiplicity of infection (MOI). A similar PDE is used for the cytokine modelling (see §2.3). The total number of extracellular virions at time t ∈ I is given by

Therefore, the number of extracellular virions resident on an epithelial cell Ω e k , k = 1, . . . , N epi , is given by

The numerical approximation to (2.2.1) is obtained using a first-order forward Euler scheme in time and a second-order central finite difference scheme in space, where each finite difference node corresponds to either a blood vessel or an epithelial cell.

Following infection, host-cells utilise pattern recognition receptors (PRRs) to detect pathogen associated molecular patterns (PAMPs) and damage associated molecular patterns (DAMPs), leading to the secretion of a variety of cytokines and chemokines which regulate the host response to the infection. The balance of these cytokines and chemokines is crucial in ensuring an effective and efficient management, as well as clearance, of the infection with minimal damage to the host. Indeed, elevated levels of many pro-and anti-inflammatory markers have been observed to correlate with severe COVID-19 [5, 41, 2, 1] . In particular, recent evidence suggests IL-6 and GM-CSF play a central role in severe COVID-19 [5] . Due to limits on computational resources, it is impractical for our model to consider all of the implicated cytokines and chemokines. Therefore, we choose to consider the following subset: type I interferon (IFN-I), because it induces the expression of IFN-stimulated genes (ISGs) that can inhibit viral entry and replication processes [42] , has been linked to genetic susceptibility [17] , and a delay in production has been suggested to increase disease severity [43, 13] as has been observed for other coronaviruses [44, 45] ; interleukins (IL-6, IL-10), because of their importance in regulating inflammation [46, 47] ; and a generic mononuclear phagocyte chemokine, so that macrophages may be directed towards an apoptotic cell, which needs to be removed via phagocytosis. Similar to the model of the extracellular virus (see §2.

2), we model the cytokines and chemokines as continuous fields using reaction-diffusion equations. Let C i : Ω T → R, i = 1, . . . , N cyt , denote the molar concentration of the i-th cytokine (units: nanomolar (nM)), where N cyt denotes the number of cytokines and chemokines considered in the model. Then the spatiotemporal evolution of the i-th cytokine or chemokine is given by

where D c i denotes the spatially-dependent diffusion coefficient; s c i denotes the source of cytokine from host cells; u c i denotes the uptake of cytokine by host cells; d c i denotes the extracellular decay; n Ω denotes the outward unit normal to Ω from ∂Ω; and C i,0 : Ω → R denotes the initial condition, which is set to zero. Cytokine sources from all cells satisfy:

where b c i denotes the basal source; φ c i denotes the maximal production rate; K i denotes the maximal production half-max; h c i denotes the maximal production Hill coefficient; µ k denotes the inhibitory halfmax; h µ k denotes the inhibitory Hill coefficient; w c k denotes the inhibitory weights; and S c i denotes the signal. Therefore, cytokines can only be produced in the presence of a sufficiently strong signal; the basal source is required to initiate the secretion, but reduces at a sufficient high level; a maximal rate is reached at sufficiently high levels; and cytokine production may be inhibited by other cytokines, but not removed entirely because of the basal production. The signal is given by

where w j denotes the signal weights; S 50,j denotes the signal half max; h s j denotes the signal Hill coefficient; and S j denotes the component (such as a specific cytokine or extracellular/intracellular Grid spacing (∆x) 0.02 mm [33] Cytokine diffusion coefficient (D c i ) 0.036 mm 2 hr −1 [33] Virus diffusion coefficient (Dv) 0.0054 mm 2 hr −1 Calculated from [36] Cytokine decay coefficient (d c i ) 0.612 hr −1 [33] Virus decay coefficient (dv) 0.02 hr −1 Chosen to be approximately an order of magnitude lower than the cytokine decay coefficient virus) that produces the signal. Note that each host cell will define individual sources, uptakes and signals depending on cell class and state. Reaction-diffusion parameters are given in Table 1 , whilst cell-specific source, uptake and signal parameters can be found in Tables 3 and 4 , respectively.

SARS-CoV-2 virions consist of four main structural proteins: the spike (S) protein, which is exposed at the surface and facilitates entry into target cells; the membrane (M) and envelope (E) proteins, which are important for morphogenesis and budding; and the nucleocapsid (N) protein, which encapsidates the viral genome [48] . The S protein facilitates entry into target-host cells by binding with host-cell surface receptors; specifically, the angiotensin-converting enzyme 2 (ACE2) receptor [48] . ACE2 is expressed in a variety of human tissues, but is abundantly expressed in upper and lower respiratory tracts [49] . Following ACE2 receptor binding, the S protein is cleaved (either at the cell surface or endosome) by cellular proteases (such as TMPRSS2 [50, 51] or cathepsins [52] ), enabling fusion to take place which leads to the uncoating of the viral envelope and the subsequent deposition of the viral genome to the host-cell cytoplasm.

The SARS-CoV-2 genome is a single-stranded, positive-sense RNA genome that is capped and polyadenylated, and possesses cis-acting elements necessary for controlling viral RNA replication. The replication of the full length, positive-sense, viral genomic RNA produces a full length intermediate negative-sense RNA genome, which serves as a template for the production of new full length, positivesense, genomic RNA [48] . The newly synthesised positive-sense RNA genomes are either packaged into new virions or are used for further replication. A feature of coronavirus replication is the discontinuous transcription process that produces sub-genomic RNA (sgRNA) [53] . The sgRNA is important for the production of the structural and accessory proteins that are required for the assembly process [48] . In general, newly assembled virions are removed from the cell by exocytosis [54] .

In the model presented, the viral entry and replication procedures are defined to be surface and bulk components of a pathway, respectively. Each epithelial cell contains its own pathway model, allowing individual variation of parameters between epithelial cells. The viral entry and replication procedures discussed above are complex and therefore, as a starting point, we modify a pre-existing simplified model, presented in [36] , to include cytokine inhibition of both viral entry and replication. More complex models are available in the literature (see e.g. [55, 56] ) and their inclusion into the presented model is a subject of future work.

For each epithelial cell Ω e k , k = 1, . . . , N epi , the viral entry model, presented in [36] , considers external unbound ACE2 receptors ([R eu ] k ), external virus-bound receptors ([R eb ] k ), internalised virus-bound receptors ([R ib ] k ) and internalised unbound receptors ([R iu ] k ), and is given by

where the functions F : R 2 → R and G : R 2 → R describe the rate reduction caused by an inhibiting cytokine (such as type I interferon) and whether the virus should be internalised, and are given by

where r denotes the rate to be reduced; C denotes the inhibiting cytokine; h c denotes the Hill coefficient; C 50 denotes the half-max value of the inhibiting cytokine; n v denotes the number of virions; and τ denotes the time at which n v first exceeds one. Therefore, the function G is a blocking function which prevents continuous internalisation of virions until a single whole virion is detected. Note that the blocking function is required to prevent spurious spread of infection through the epithelial monolayer at low levels of extracellular virus. Similarly, for each epithelial cell Ω e k , k = 1, . . . , N epi , the viral replication model, presented in [36] , considers fully endocytosed virions ([V ] k ), uncoated viral RNA ([U ] k ), synthesised viral RNA ([R] k ), synthesised viral proteins ([P ] k ) and fully assembled virions ([A] k ), and is given by

where exocytosis of newly assembled virions occurs once a single newly assembled virion is detected. Therefore, the source and uptake of extracellular virus is given by

where |Ω e k | denotes the area (or volume) of Ω e k . Figure 2 illustrates the viral entry (2.4.1) and replication (2.4.3) dynamics of a single epithelial cell infected by a single virion in the absence of viral diffusion (D v = 0), and compares against the solution of the ODE systems 1 given by (2.4.1) and (2.4.3) . Clearly, excellent agreement is seen over both short and long timescales for both entry and replication dynamics. Viral entry and replication parameters are given in Table 2 . The numerical solution of (2.4.1) and (2.4.3) is obtained using a fourth-order Runge-Kutta scheme, which is available as part of the GNU Scientific Library [57] . Note that we have not observed any numerical instabilities due to discontinuous export.

As discussed in §2.1, the grid of epithelial cells Ω e is assumed to be uniform and therefore, each epithelial cell is chosen to be a square with width and height of 20 µm. Each epithelial cell can exist in five possible states: healthy, infected, infectious, apoptotic and removed. A healthy epithelial cell is placed in an infected state as soon as the number of internal bound receptors are non-zero ([R ib ] k > 0, k = 1, . . . , N epi ). Therefore, an infected epithelial cell is an epithelial cell which has been infected by the virus but is not yet virus-producing (infectious). Due to the blocking function (G in (2.4.1)), an infected epithelial cell becomes infectious when a single whole virion is detected on the inside of the cell. As discussed in §2.4, each epithelial cell has a copy of the viral entry (2.4.1) and replication (2.4.3) models. An infectious epithelial cell may be placed into an apoptotic state due to either high levels of type 1 interferon (IFN-I) or high levels of intracellular, newly assembled, virions ([A] k ). In both cases, a Hill function H : R + → [0, 1] is used to determine the likelihood of triggering apoptosis in an epithelial cell, and is given by [U ]

[R]

[P ]

[A]

(d) Viral replication T = 48 hrs Table 2 .

where r apop max,S denotes the maximum apoptosis rate for signal S; S apop k denotes the apoptosis signal (such as IFN-I or [A] k ) for the k-th epithelial cell; h apop S denotes the Hill coefficient for signal S; and S apop k,50 denotes the half-max value. An epithelial cell will then become apoptotic provided r < H(S apop k ), where r ∈ [0, 1] is a uniform random number. Apoptotic epithelial cells are targeted for phagocytosis by active macrophages (see §2.6.4), and are therefore placed into a removed state upon successful phagocytosis.

We do not attempt to model homeostasis and therefore, healthy epithelial cells do not secrete cytokines. However, once a cell is infected, we assume that it produces cytokines. For computational reasons, we assume that all cytokines considered in the model can be produced by infected epithelial cells, according to the procedure detailed in §2.3, but in varying quantities (see Table 3 ). As discussed in §2.3, each cytokine produced requires a specific signal with specific parameters for that cytokine. For each cytokine produced by infected epithelial cells, the signal is defined as the amount of intracellular, newly assembled, virions, with the exception of IL-10 whose signal is the level of IL-6. Epithelial cells in the apoptotic state can only produce the generic mononuclear phagocyte chemokine, which chemotactically directs macrophages to the apoptotic epithelial cell so that phagocytosis may take place. Epithelial cells in a removed state cannot produce any cytokines. Additionally, the diffusion coefficients of cytokines and extracellular virus at the location of a removed epithelial cell are set to zero, to simulate an empty space on the grid. Binding rate (r bind ) 0.06 hr −1 Calculated from [36] Recycle rate (r recycle ) 0.6 hr −1 Calculated from [36] Endocytosis rate (r endo ) 0.6 hr −1 Calculated from [36] Release rate (r release ) 0.06 hr −1 Calculated from [36] Uncoating rate (r U ) 0.6 hr −1 Calculated from [36] RNA synthesis rate (r P ) 0.6 hr −1 Calculated from [36, 35] RNA decay rate (λ R ) 0.021 hr −1 Calculated from [36] Protein synthesis rate (r S ) 0.6 hr −1 Calculated from [36] Assembly rate (r A ) 0.6 hr −1 Calculated from [36, 35] Protein decay rate (λ P ) 0.021 hr −1 Calculated from [36] Export rate (r E ) 0.6 hr −1 Calculated from [36, 35] Maximum RNA synthesis (rmax) 0.3 hr −1 Calculated from [35] RNA synthesis half-max (r 1/2 ) 2000 [35] Half Table 3 . Table 3 . 

Although the model presented supports macrophages, monocytes, neutrophils, and natural killer cells, we choose to focus only on macrophages for the purposes of this article. The reason for this choice, is primarily due to model complexity, but is also informed by mounting evidence suggesting an important role for macrophages [58, 59] . Evidence suggests that there is an early loss of resident alveolar macrophages, as well as a recruitment of inflammatory interstitial and monocyte-derived macrophages (see e.g. [60] and the references therein). However, for our initial model, we assume that for each class of immune cell (such as a macrophage) there is only a single type (such as alveolar macrophage) and we do not distinguish between phenotypes. We recognise the clear limitation to this assumption, and extending the model to incorporate different types and phenotypes is a subject of future work. Immune cells may, therefore, exist in three possible states: resting, active and apoptotic. The transition from resting to active state, and vice-versa, is discussed in §2.6.3. Similar to an epithelial cell, we assume an immune cell may be placed in an apoptotic state by high levels of IFN-I, or by high levels of intracellular virions. Apoptotic immune cells may be removed from the grid following successful phagocytosis by phagocytic immune cells (such as macrophages). Again, as we do not attempt to model homeostasis, resting immune cells do not secrete cytokines. However, once an immune cell becomes active, then cytokines are secreted. Similar to epithelial cells, for computational reasons we assume that all cytokines considered in the model can be produced by immune cells, according to the procedure detailed in §2.3, but in varying quantities (see Table 4 ). Therefore, each immune cell is considered a major source of cytokines, but not the only source. As discussed in §2.3, each cytokine produced requires a specific signal with specific parameters for that cytokine and this is discussed in more detail in §2.6.4. Furthermore, we assume that immune cells in the apoptotic state cannot produce any cytokines, nor can they migrate.

Each immune cell can migrate from epithelial cell to epithelial cell and is only allowed to reside on a single epithelial cell at any one time (lattice-based movement across finite difference grid). The migration may be pseudo-random or may be directed in response to cytokine levels. Each immune cell migrates according to a characteristic movement rate, which is allowed to vary for each immune cell (see Table  4 ). Note that we assume that apoptotic immune cells cannot migrate.

Let Ω nb k ⊂ Ω denote the set of immediate neighbours of epithelial cell Ω e k , k = 1, . . . , N epi , such that

where Ω nb k,j ⊂ R 2 , j = 1, . . . , 8, can be either an epithelial cell or a blood vessel cross-section. Note that we do not allow immune cells to migrate onto blood vessels. If an immune cell is resident on the epithelial cell Ω e k , k = 1, . . . , N epi , then for each neighbour Ω nb k,j , j = 1, . . . , 8, we calculate a migrating signal S M k,j , which is defined as the weighted sum of chemokines on the j-th neighbour of the k-th epithelial cell. Note that the chemokines involved are potentially specific for each immune cell, and the weights enable a migration bias towards certain chemokines for a particular immune cell. For each immune cell, we then define the total migrating signal from the immediate neighbours of the k-th epithelial cell as

where the superscript (·) m b denotes a movement bias, which is allowed to vary for each immune cell. If the total migrating signal from the immediate neighbours is zero (S M k = 0), or if the immune cell state is resting, then the immune cell migrates randomly. Otherwise, the immune cell will migrate chemotactically according to the following procedure:

denote the sequence of partial sums given by

and let r ∈ [0, 1] denote a uniform random number; then for each i = 1, . . . , 8, an immune cell which is resident on the k-th epithelial cell will migrate to the i-th neighbour provided S M k,i > r · S M k and there is sufficient space available on the i-th neighbour. Note that the available space is determined by the difference between the area of the i-th neighbour and the sum of the area of any immune cells that are resident on the i-th neighbour.

The pseudo-random migration of the immune cells is illustrated in Fig. 4 (a), where the immune cell movement is largely clustered around its initial location before moving away. Chemotactic migration is illustrated in Fig. 4 (b), where directed migration of the immune cell towards the signal ( Fig. 4 (d)) can be seen.

The immune cells are recruited through the blood vessels, where, for each t ∈ I, each blood vessel is checked to see whether a single immune cell of each available class (for example, a macrophage) needs to be recruited from that vessel. For each neighbour of the k-th blood vessel, we calculate a recruitment signal (S R k,j ), which is defined as the difference between the weighted sum of pro-inflammatory cytokines and the weighted sum of anti-inflammatory cytokines on the j-th neighbour of the k-th blood vessel. Similar to the migration procedure above, the cytokines involved are potentially specific for each immune cell, and the weights enable a recruitment bias towards certain cytokines for a particular immune cell. The total recruitment signal from the immediate neighbours of the k-th blood vessel is given by

The Hill function H : R + → [0, 1] given in (2.5.1) is then used to determine the likelihood that an immune cell will be recruited from the k-th blood vessel. An immune cell will be recruited if Note that a grid size of 120 × 120 has been used; a red dot denotes any immune cell; a cyan dot denotes a blood vessel; and a black line denotes an immune cell path.

r ∈ [0, 1] is a uniform random number, and if there is sufficient space available on a randomly chosen neighbour of the k-th blood vessel. Note that if the recruitment signal is negative, then an immune cell is not recruited. All immune cells are recruited in a resting state and cannot be recruited onto a blood vessel.

Remark 1. We note that the recruitment procedure is sensitive to the choice of Hill function parameters. If the recruitment half-max, or the Hill coefficient, are too low then immune cell recruitment may be seen early in the simulation, which is undesirable as we expect the resident immune cells to be the only respondants during the early stages of infection. However, if the recruitment half-max is too high, then we may not observe any recruitment. Furthermore, if the Hill coefficient is too large, then switch-like recruitment may take place, where no recruitment is seen until the half-max, followed by saturating recruitment soon after, leading to a flooding of immune cells which can significantly affect computational time. Therefore, the Hill function parameters require manual tuning for each set of simulation parameters. Reducing this parameter sensitivity is a subject of future work. Figure 4 (c) illustrates the recruitment of immune cells, where following a sufficiently strong chemotactic signal (Fig. 4(d) ) immune cells enter the domain from a nearby vessel, before chemotactically migrating towards the signal, as expected.

The activation and deactivation of the immune cells is carried out in two ways: according to local cytokine levels; and probabilistically if the extracellular/intracellular viral load is sufficient. Specifically for macrophages, if the internalised viral content exceeds one, then the macrophage will activate provided r < p activate , where r ∈ [0, 1] is a uniform random number and p activate is the probability of activation. In practice, many cytokines may be involved in the activation of the immune cells. Consequently, an immune cell which is resident on the k-th epithelial cell has an activating signal (S A k ), which is defined as the difference between the weighted sum of pro-inflammatory cytokines and the weighted sum of anti-inflammatory cytokines, where the cytokines involved are potentially specific for a given immune cell, similar to the migration and recruitment procedures above. Once again, the weights enable an activation and deactivation bias towards certain cytokines for a particular immune cell. The Hill function H : R + → [0, 1] given in (2.5.1) is then used to determine the likelihood that an immune cell resident on the k-th epithelial cell will activate or deactivate. Therefore, an immune cell will activate if r < H(S A k ), where r ∈ [0, 1] is a uniform random number. Note that a macrophage will deactivate provided r > H(S A k ), r 2 > p activate and the internalised viral content is less than one, where r, r 2 ∈ [0, 1] are uniform random numbers.

Each macrophage is assumed to be a circle, with a diameter of 20 µm. Initially, a user-defined number of resident macrophages are randomly placed in the tissue environment. Aside from the secretion of cytokines, the main role of macrophages is to patrol tissue environments and clear away any dead tissue and cells, as well as extracellular debris and pathogens [61] . In our model, both resting and active macrophages may internalise extracellular virus, whilst only active macrophages may phagocytose apoptotic epithelial and immune cells. Additionally, to aid viral clearance, we assume that an active macrophage may phagocytose an infected epithelial cell by implicitly triggering apoptosis, however this occurs infrequently in our model. As mentioned at the beginning of §2.6, we do not attempt to model homeostasis and therefore, only active macrophages produce cytokines. We assume that active macrophages produce all the cytokines considered in the model, but at varying quantities. As discussed in §2.3, cytokine production requires a specific signal with specific parameters for that cytokine. For the cytokines IFN-I and IL-6, the signal is defined as the internalised viral content, whereas for the cytokine IL-10, the signal is defined as the pro-inflammatory cytokine IL-6.

The internalisation of extracellular virus is governed by Table 4 ), then the macrophage may rupture, releasing the viral contents into the extracellular environment. Phagocytosis of an infected or apoptotic epithelial cell, as well as apoptotic immune cells, is assumed to take place over a time scale which is set by the target cell (for example, an apoptotic epithelial cell). Once the phagocytosis process has begun, the target cell is placed into a temporary phagocytosed state, where all cellular processes stop. Additionally, the macrophage performing the phagocytosis enters a temporary phagocytosing state, where cytokines may still be produced but the cell cannot change state or migrate until the phagocytosis process is complete. Once phagocytosis is complete, the target cell may be placed into a removed state (in the case of epithelial cells) or it is removed from the grid (in the case of immune cells).

Resting macrophages may only migrate pseudo-randomly, whilst active macrophages may migrate chemotactically, as well as pseudo-randomly, according to the procedure detailed in §2.6.1. The chemotactic migration of macrophages is controlled solely by a generic mononuclear phagocyte chemokine, which is released by apoptotic cells (both epithelial and immune). Recruitment of resting macrophages from the blood vessels is carried out according to the procedure detailed in §2.6.2. Recruitment of the macrophages is encouraged by: IL-6; and the generic mononuclear phagocyte chemokine, whilst the anti-inflammatory cytokine IL-10 discourages recruitment.

Activation and deactivation of a macrophage is carried out according to the procedure detailed in §2.6.3. Activation of resting macrophages is encouraged by: internalised viral content [V ]; IFN-I; IL-6; and the generic mononuclear phagocyte chemokine, whilst the anti-inflammatory cytokine IL-10 discourages the pro-inflammatory activation. Therefore, we are only considering classical activation of the macrophages in our model. All cytokines involved in the activation and deactivation of macrophages are obtained from the location where the macrophage resides. Similar to phagocytosis, apoptosis of a macrophage is assumed to take place over a time scale which is set by the macrophage. Once the apoptosis process has begun, the macrophage is placed into a temporary apoptosed state, where we assume that cytokines may still be released but the cell cannot change state, migrate or carry out phagocytosis; this is a reasonable approximation, particularly during the early stages of the apoptotic programme. Once apoptosis has been completed, the macrophage is placed into a apoptotic state, where the generic mononuclear phagocyte chemokine is released to encourage efferocytosis, but all other cytokine secretions stop, and the macrophage may not migrate nor carry out phagocytosis. Additionally, once apoptosis is completed, we assume that a random portion of the macrophage's internalised viral content is scattered to the neighbouring epithelial cells.

Some of the parameter values, which are given in Tables 1-4 , have been obtained from the modelling literature: namely, Bowness et al. [33] , Sego et al. [35] , and Getz et al. [36] , which in turn obtained their parameters either heuristically, from other modelling work, or by estimation, approximation and calibration from experimental studies. Therefore, we do not claim all of the parameters are specific to SARS-CoV-2, but can be from related viruses such as SARS-CoV-1. Furthermore, not all parameter values could be directly related. Therefore, the reader will notice that specific language is used when justifying a parameter value; specifically, assumed, chosen, heuristically chosen, heuristically estimated from and calculated from. Parameter values which are chosen have not been heuristically varied, instead they are chosen for a particular reason or to simulate a specific behaviour; heuristically chosen parameters are those which have been chosen from computational experiments; heuristically estimated from denotes parameters which were obtained from the literature but then varied heuristically for our model; and calculated from denotes parameters which were obtained from the literature but altered, due to a change in units, for example.

The code used in this article can be found on GitHub:

https://github.com/Ruth-Bowness-Group/CAModel. The data used in this article has been archived within the Zenodo repository: https://doi.org/10.5281/zenodo.6514656.

To investigate the influence of the initial viral deposition, we consider three values for the initial multiplicity of infection (MOI = 0.001 (= 11 virions), 0.01 (= 102 virions), 0.1 (= 1015 virions) ). Due to the possibility of the immune response being influenced by the grid size, each initial viral deposition is randomly distributed over a subdomain located in the centre of the tissue environment; the full grid size is 200 × 200, whilst the central subdomain has a grid size 100 × 100, which leaves a buffer around the initial subdomain for the infection to spread into. Therefore, the number of virions to be deposited is calculated using the initial MOI and the number of epithelial cells in the subdomain. Rigorous investigation into the grid size dependency is a subject of future work. Figure 5 shows the spread of infection at the end of the simulation (T = 120 hours) as the MOI is increased (baseline parameter values are given in Tables 1-4 ). In Fig. 5(a) , we can clearly see nonuniform spread of infection which is caused by each epithelial cell possessing a different number of ACE2 receptors. Infected (not virus producing) epithelial cells form a thin boundary around infectious (virus producing) epithelial cells. We see a high number of macrophages in either a resting, active or apoptotic state. Furthermore, we observe that the majority of macrophages at the site of infection are either in an active or apoptotic state, whilst the majority of resting macrophages are either local to the boundary of infection or away from it. As the MOI is increased (Figs. 5(b)-(c) ), we see that the spread of infection increases, as well as increased number of active macrophages and, in particular, apoptotic macrophages. For MOIs = 0.01, 0.1 (Figs. 5(b) ,5(c)), we see an almost square spread of infection due to the high viral deposition for the number of epithelial cells in the initial subdomain. Fig. 6(a) , we see an increasing number of infectious epithelial cells by the end of the simulation, for increasing MOI. Additionally, there is an initial delay (≈ 12 hours) before infectious epithelial cells are first found; this initial delay is seen for all MOI values but is most easily observed for MOI = 0.1, and is a consequence of the timescale of viral entry and replication, as well as the blocking function (see §2.4), which allows a whole virion to be produced before being exocytosed. Interestingly, following the initial delay, we observe a second, longer delay before a continuous increase in infectious epithelial cells is observed, which is likely to be caused by both the blocking function and IFN-I inhibition (see §2.4 and Fig. 8(a) ). Epithelial cells may be placed into an apoptotic state due to high intracellular viral loads or due to high IFN-I; in Fig. 6(b) , we see an increase in apoptotic epithelial cells after t = 24, where, for each MOI, the start of the increase correlates with both an increase in total (across all epithelial cells) intracellular viral load (Fig. 7(a) ) and extracellular IFN-I ( Fig. 8(a) ). Furthermore, we observe an increase in the number of apoptotic epithelial cells with increasing MOI. Epithelial cells may be placed into a removed state after successful macrophage phagocytosis; in Fig.  6 (c) we see an increase in removed epithelial cells after t = 48 hours, which correlates with an increase in active macrophage numbers (Fig. 6(e) ) for each MOI considered. Interestingly, for the largest MOI considered, we see a spike in resting macrophages at around t = 48 hours, which is not observed for the other MOIs (Fig. 6(d) ), nor for active and apoptotic macrophage numbers (Figs. 6(e),(f)), and may be caused by the sharper increase in apoptotic epithelial cells at around t = 48 hours (Fig. 6(b) ), as well as extracellular IL-6 levels ( Fig. 8(b) ). Alongside the increase in apoptotic epithelial cells, we see an increase in apoptotic macrophages (Fig. 6(f) ), which in conjunction with increases in IL-6 levels, likely contributes to the recruitment of resting macrophages observed in Fig. 6(d) . Moreover, we observe a decline in resting macrophage numbers towards the end of the simulation (Fig. 6(d) ), which is likely caused by increased activation as well as increasing IL-10 levels (Fig. 8(c) ), which discourages both recruitment and activation. Generally, we see an increase in resting, active and apoptotic macrophages with increasing MOI. illustrates the intracellular and extracellular viral load per grid, as the simulation progresses, for each MOI considered. In Fig. 7(a) a local maximum can be observed for all MOI values; as the MOI is increased, we see the local maximum appear earlier in the simulation, but the magnitude of the local maximum does not monotonically increase as the MOI is increased. Indeed, for MOI = 0.01, we see a greater local maximum than when MOI = 0.1. Our hypothesis is that the local maximums appear due to a sudden increase in IFN-I and are therefore, affected by both the delay in epithelial IFN-I secretion as well as the epithelial IFN-I secretion rate, and is discussed in more detail in §3.2. It is interesting, therefore, that the local maximum may be overcome to produce a rebound in intracellular virions per grid, where the higher the MOI, the steeper the rebound, and may be caused by viral replication inside epithelial cells where the local extracellular IFN-I levels are lower, or because the IFN-I inhibition of viral entry and repliation has saturated whilst the local extracellular viral load is increasing. In Fig.  7(b) , it appears that for t < 48 hours, the extracellular viral loads are approximately similar for each MOI considered. However, for t > 48 hours, the viral loads begin to differ (although not substantially until t = 96 hours) and by the end of the simulation we observe an increasing number of extracellular virions as the MOI is increased, as expected. Around t = 96 hours, for MOI = 0.001, we see a temporary flattening of the number of extracellular virions which correlates with the local maximum in intracellular virions per grid (Fig. 7(a) ). Figure 8 illustrates the levels of each cytokine considered as the simulation progresses, for each MOI value. In all cases, we see an increasing amount of cytokine by the end of the simulation as the MOI is increased; similarly, the initial increase in cytokine level occurs earlier in the simulation, and generally has a steeper increase, as the MOI is increased. For MOI = 0.1, we see a sharp increase IFN-I levels around t = 48 hours, which correlates with the observed local maximum in the intracellular viral load per grid ( Fig. 7(a) ). However, we do not observe any changes in cytokine levels that correspond to the rebound in intracellular viral load per grid. IL-10 levels (Fig. 8(c) ) seem to flatten off towards the end of the simulation for the largest MOI considered, due to a short flattening in IL-6 levels ( Fig. 8(b) ), whilst are still increasing for the other MOI values. [44, 45] , which has been suggested to increase disease severity [42, 21] . Therefore, we investigate the impact of a delayed IFN-I secretion from epithelial cells on the spread of infection, epithelial and macrophage numbers, as well as viral load and cytokine levels. We consider four values of the epithelial IFN-I secretion delay (S 50,V = 1, 10, 100, 1000 virions) so that IFN-I is only secreted from infected and infectious epithelial cells when 1, 10, 100 or 1000 virions are detected in the cell cytoplasm. Note that we fix the MOI = 0.01 in this section, and do not consider delays in IFN-I inhibition of viral entry and replication pathways. Figure 9 illustrates the spread of infection at the end of the simulation (T = 120 hours), for fixed MOI = 0.01, as the epithelial IFN-I secretion delay is increased. Clearly, we observe an increase in the spread of infection for longer epithelial IFN-I secretion delays, as the viral entry and replication processes are able to proceed for longer without being inhibited by IFN-I. Similar to Fig. 5 , we observe mostly active and apoptotic macrophages at the site of infection, with resting macrophages mostly localised to the boundary of the infection and away from it. As the epithelial IFN-I secretion delay is increased, we see an increasing number of resting and active macrophages but a decrease in apoptotic macrophages Tables 1-4 .

due to lower levels of IFN-I. Figures 10(a) -(c) illustrate the number of epithelial cells in an infectious, apoptotic or removed state, respectively, for each value of the epithelial IFN-I secretion delay. As expected, in Fig. 10(a) , we see an increasing number of infectious epithelial cells by the end of the simulation, for increasing IFN-I secretion delay due to the enhanced spread illustrated in Fig. 9 . Additionally, for t < 48 hours, we do not observe significant differences between the number of infectious, apoptotic or removed epithelial cells, for each IFN-I secretion delay. Epithelial cells may be placed into an apoptotic state due to high intracellular viral loads or due to high IFN-I; in Fig. 10(b) we see an increase in apoptotic epithelial cells at the end of the simulation, as the IFN-I secretion delay is increased, due to an increase in intracellular viral load (see Fig.  11(a) ). Epithelial cells may be placed into a removed state after successful macrophage phagocytosis; in Fig. 10(c) we see an increase in removed epithelial cells at the end of the simulation, as the IFN-I secretion delay increases, which correlates with an increase in active macrophage numbers (see Fig.  10 (e)) for each IFN-I secretion delay considered. Interestingly, the number of resting macrophages seems to approximately flatten after t > 72 hours (Fig. 10(d) ). Furthermore, reduced apoptotic macrophages are seen for the longest epithelial IFN-I secretion delay (Fig. 10(f) ), as illustrated in Fig. 9 . Figure 11 illustrates the intracellular and extracellular viral load per grid as the simulation progresses, for each IFN-I secretion delay considered. For both the intracellular (Fig. 11(a) ) and extracellular ( Fig.  11(b) ) viral loads per grid, we see a very significant increase by the end of the simulation, for the longest IFN-I secretion delay (= 1000 virions). Additionally, the local maximum that can be observed in the intracellular viral load per grid for epithelial IFN-I secretion delays of 1, 10, 100 virions, cannot be observed for an IFN-I secretion delay of 1000 virions; instead, we observe a decrease in the growth rate of intracellular viral load per grid, further suggesting that IFN-I (relative to the amount of intracellular virus) is responsible for the local maximum in the intracellular viral load per grid observed earlier ( Fig.  7(a) ). Therefore, the rebound in intracellular viral load per grid observed earlier (Fig. 7(a) ), is still Tables 1-4 .

observed for IFN-I secretion delays of 1, 10, 100 virions, but not for a IFN-I secretion delay of 1000 virions ( Fig. 11(a) ). Figure 12 illustrates extracellular IFN-I, IL-6, IL-10 and generic mononuclear phagocyte chemokine levels as the simulation progresses, for each epithelial IFN-I secretion delay considered. As expected, IFN-I levels decrease as the secretion delay increases (Fig. 12(a) ), but IL-6 levels increase as IFN-I secretion delay increases (Fig. 12(b) ) due to the increasing spread of infection and thus, infectious epithelial cells (Figs. 9 and 10(a) ). Correspondingly, we see an increase in IL-10 levels as IFN-I secretion delay increases (Fig. 12(c) ) due to the increase in IL-6 levels; however, for the longest IFN-I secretion delay, IL-10 levels seem to be declining towards the end of the simulation. Additionally, we observe an increase in generic mononuclear phagocyte chemokine as IFN-I secretion delay increases (Fig. 12(d) ) due to an increase in the total number of apoptotic epithelial cells and macrophages (Figs. 10(b),(f) ). Tables 1-4 .

To further investigate the influence of IFN-I secretion on the spread of infection, epithelial and macrophage numbers, as well as viral load per grid and cytokine levels, we consider three values of the epithelial IFN-I secretion rate (φ c i = 0.00126, 0.0126, 0.126 nM hr −1 ), for fixed MOI = 0.01 and epithelial IFN-I secretion delay = 1 virion. Figure 13 illustrates the spread of infection at the end of the simulation (T = 120 hours), as the epithelial IFN-I secretion rate increases. Clearly, as the epithelial IFN-I secretion rate is increased, we see a decreased number of infectious epithelial cells, as well as a decrease in resting, active and apoptotic macrophages at the site of infection; however, there is an increase in macrophage apoptosis away from the site of infection (Fig. 13(c) ) due to larger amounts of extracellular IFN-I. Tables 1-4 . Figures 14(a),(b) ,(c) illustrate the infectious, apoptotic and removed epithelial cell numbers, as the simulation progresses, for each IFN-I secretion rate considered. As shown in Fig. 13 , we see a decreasing number of infectious epithelial cells, at the end of the simulation, as IFN-I secretion rate increases ( Fig.  14(a) ). However, for t < 48 hours we do not see any significant differences in the number of infectious, apoptotic or removed epithelial cells, for each IFN-I secretion rate considered. As epithelial cells may be placed into an apoptotic state because of high IFN-I or high intracellular virus, we observe higher numbers of apoptotic epithelial cells for the smallest and largest IFN-I secretion rates considered, when compared to the median IFN-I secretion rate; however, for the largest IFN-I secretion, we see a significantly higher number of apoptotic epithelial cells (Fig. 14(b) ), by the end of the simulation. In Fig. 14(c) , we see a decreasing number of removed epithelial cells at the end of the simulation, as the IFN-I secretion rate is increased, which correlates with the decreasing levels of active macrophages (Fig. 14(e) ). For the largest IFN-I secretion rate, we do not observe an expansion of resting macrophages (Fig. 14(d) ) due to the lower levels of IL-6 ( Fig. 16(b) ); therefore, for t > 48 hours, we see a decrease in the number of resting macrophages as the simulation progresses due to activation and apoptosis. For the other IFN-I secretion rates, we observe an expansion of resting macrophages with the greatest expansion observed for the smallest IFN-I secretion rate. Due to the lower levels of macrophages observed for the highest IFN-I secretion rate, we observe lower levels of apoptotic macrophages (Fig. 14(f) ); however, for the other IFN-I secretion rates considered, we do not see any difference in apoptotic macrophage numbers over the course of the simulation. Figure 15 illustrates the intracellular and extracellular viral loads per grid, as the simulation progresses, for each epithelial IFN-I secretion rate considered. In Fig. 15(a) , we see an increasing intracellular viral load per grid with increasing IFN-I secretion rate, by the end of the simulation. Interestingly, for the smallest IFN-I secretion rate (0.00126 nM hr −1 ) we do not observe a local maximum in intracellular viral load, which can be seen for the other IFN-I secretion rate values; instead, similar to the longest epithelial IFN-I secretion delay (Fig. 11(a) ), we observe a decrease in the growth rate of intracellular viral load per grid. Additionally, following the local maximum, the rebound in the intracellular viral load which has been observed previously (e.g. Fig. 7(a) ), is not clearly observed for the largest IFN-I secretion value. Figure 16 illustrates the IFN-I, IL-6, IL-10 and generic mononuclear phagocyte chemokine levels, as the simulation progresses, for each epithelial IFN-I secretion rate considered. As expected, increasing the IFN-I secretion rate increases the IFN-I levels ( Fig. 16(a) ); higher IFN-I levels reduces the number of infectious epithelial cells (Fig. 14(a) ), resulting in decreasing IL-6 and IL-10 levels as the IFN-I secretion rate is increased (Figs. 16(b),(c) ). We observe slightly higher levels of generic mononuclear phagocyte chemokine for the largest IFN-I secretion rate (Fig. 16(d) ), which correlates with a higher total number of apoptotic cells (Fig. 14(b) ,(f)); however, we do not correspondingly observe an increase in macrophage numbers for increasing IFN-I secretion rates (Fig. 14(d) ), even though the generic mononuclear phagocyte chemokine is involved in the recruitment of macrophages, due to insufficient IL-6 levels. 

We consider variations of two macrophage parameters in our simulation: macrophage virus internalisation rate and macrophage activation half-max. Figure 17 illustrates the spread of infection at the end of the simulation, as the macrophage virus internalisation rate increases. There does not appear to be any significant differences between the spread of infection at the end of the simulation. This is confirmed in Figs. 18(a),(b) ,(c), where we see similar amounts infectious, apoptotic and removed epithelial cells, by the end of the simulation, as the macrophage virus internalisation rate increases. However, in Fig. 17 , we do observe (slight) decreases in macrophage numbers as the virus internalisation rate increases. This is confirmed in Figs. 18(d) ,(e),(f), where we observe decreasing levels of resting, active and apoptotic macrophages as the virus internalisation rate increases. Additionally, we observe only slight differences in intracellular and extracellular viral loads per grid (Fig. 19) , as well as only slight differences in the cytokine levels (Fig. 20) . Figure 22 illustrates the number of infectious, apoptotic and removed epithelial cells, as well as resting, active and apoptotic macrophage numbers, as the simulation progresses, for increasing activation halfmax values. For the lowest activation half-max value (meaning macrophages can activate readily), we observe higher numbers of infectious and apoptotic epithelial cells (Figs. 22(a),(b) ); however, we see a lower level of removed epithelial cells (Fig. 22 (c)) due to lower numbers of active macrophages (Fig.  22(e) ). For the lowest activation half-max value, there does not appear to be a significant expansion of macrophages (Figs. 22(d) ,(e)). Interestingly, for the highest activation half-max value considered, there appears to be a flattening of infectious and apoptotic epithelial cells (Figs. 22(a),(b)); however, we clearly see a larger number of removed epithelials (Fig. 22 (c)) due to higher levels of active macrophages ( Fig. 22(e) ), which are most likely observed due to the higher levels of resting macrophages (Fig. 22(d) ). Figure 23 illustrates the intracellular and extracellular viral loads per grid, as the simulation progresses, for all activation half-max values considered. The intracellular viral loads per grid follow similar trajectories up to ≈ t = 72 hours, for each activation half-max value considered; however, for t > 72 hours, the viral load per grid for the largest activation half-max value reaches a lower minimum than for the other activation half-max values ( Fig. 23(a) ). However, there does not appear to be any significant changes in the IFN-I levels around t = 72 hours that would explain the lower minimum observed in the intracellular viral load for the largest activation half-max value. Similarly, we observe a lower extracellular viral load per grid for the largest activation half-max value ( Fig. 23(b) ), due to the larger number of active macrophages (Fig. 22(e) ). Figure 24 illustrates the IFN-I, IL-6, IL-10 and generic mononuclear phagocyte chemokine levels, as the simulation progresses, for each activation half-max value. The trajectories for IFN-I and IL-6 are identical (Figs. 24(a),(b) ), where the highest levels of IFN-I and IL-6 are observed for the lowest activation half-max value. The higher levels of IL-10 that are observed for the largest activation halfmax value (Fig. 24(c) ), are most likely due to the higher numbers of active macrophages observed for the largest activation half-max value. Additionally, higher levels of the generic mononuclear phagocyte chemokine are observed for the largest activation half-max value (Fig. 24(d) ), due to the higher total number of apoptotic cells (Fig. 22(b) ,(f)). 

The experimental study of host-pathogen dynamics generally involve utilising a combination of cell lines (see e.g. [62] ), organoids (see e.g. [63, 64, 65] ), animal models (see e.g. [66, 44, 67] ), as well as patient and post-mortem data (see e.g. [40, 1, 6, 4] ). Mathematical and computational models offer an additional perspective which can supplement experimental studies to enhance the understanding of host-pathogen dynamics (see e.g. [21, 22, 23, 24] Tables 1-4. focus on macrophages only, as well as a limited set of cytokines, in order to keep the model complexity to a minimum. The model is then used to study the influence of initial viral deposition, by increasing the initial MOI (Figs. [5] [6] [7] [8] ; the influence of delayed IFN-I secretion from epithelial cells, as well as the Understanding how the initial viral deposition relates to more severe infection and disease outcomes is complicated by the highly nonlinear interactions between the virus, and the innate and adaptive immune systems. Never-the-less, many studies have observed correlations between a high viral load and more severe disease outcomes (see e.g. [1, 68, 69] ). Therefore, although our model is limited at this stage (see §4.1), we investigated the role of increasing initial viral deposition (through increasing the initial MOI) on the spread of infection. An increase in infection was observed (Fig. 5) , with a higher number of infectious, apoptotic and removed epithelial cells being seen for the highest MOI value considered (Figs. 6(a)-(c)). The increase in infection correlated with an increase in interferon, cytokine and chemokine levels (Fig. 8) . A similar correlation has been observed from longitudinal studies (see e.g. [1] ). However, care should taken when interpreting the observed correlation; our model is by necessity a simplification, therefore we do not consider other macrophage phenotypes, other immune cell types, or other anti-inflammatory actions, which one would expect to further regulate the interferon, cytokine and chemokine levels. Unexpectedly, we observed a local maximum in the intracellular viral load per grid (Fig. 7) , which moved earlier in the simulation as the MOI increased. Within the confines of our model, we hypothesise that a local maximum can only be observed if the export of the intracellular virions exceeds the production (which may happen due to an increase in IFN-I levels relative to the amount of intracellular virus), or if a cluster of infectious epithelial cells have been removed by apoptosis or macrophage phagocytosis, or a combination of both. At the time of the observed local maximum in intracellular viral load per grid, we observed an increase in IFN-I levels but not significant increases in removed or apoptotic epithelial cells, suggesting that the local maximum arises due to IFN-I inhibition of viral entry and replication.

Type I interferons (IFN-I) are potent antiviral cytokines that are produced from a variety cell types in response to a viral infection [42] . Given their potent antiviral action, viruses have evolved to not only evade host-cell detectors so that IFN-I induction is suppressed, but also interfere with interferonstimulated gene (ISG) induction [14] . Consequently, a delay in IFN-I activity has been suggested to increase disease severity [43, 13] , as has been observed for other coronaviruses [44, 45] . Therefore, we investigated the role of increasing delay of IFN-I production from epithelial cells, as well as the magnitude of IFN-I secretion from epithelial cells, on the spread of infection. Increasing the IFN-I secretion delay from epithelial cells increases the spread of infection (Fig. 9) , with higher infectious, apoptotic and removed epithelial cells (Figs. 10(a)-(c)) being observed for the longest IFN-I secretion delay. Furthermore, for the longest IFN-I secretion delay, reduced levels of IFN-I were observed (as expected) accompanied by higher levels of IL-6 ( Figs. 12(a),(b) ). Although increased levels of IL-6 have been consistently implicated in severe COVID-19 (see e.g. [6, 1, 7, 8, 9] ), for our model the increase in IL-6 level is merely due to a larger infection. The inclusion of more cytokines and immune cell types, as well as more detailed viral entry and replication pathways, to our model would enable us to ). Note that the MOI = 0.01 and the other baseline parameters are given in Tables 1-4. of infection (Fig. 13) . Furthermore, for the highest epithelial IFN-I secretion rate considered, we did not observe a significant rebound in level of intracellular viral load per grid ( Fig. 15(a) ), enhancing the suggestion that viral replication in the presence of low extracellular IFN-I levels is responsible for the observed rebound in intracellular viral load per grid.

In our model, recruitment of resting macrophages to the site of infection is primarily encouraged by IL-6, with help from the generic mononuclear phagocyte chemokine, and discouraged by IL-10. Therefore, generally, when there is a larger amount of infection, we see higher levels of IL-6 and consequently, greater recruitment of resting macrophages (Figs. 6(d), 10(d), 14(d)). The increase in resting macrophage numbers results in increasing numbers of active and apoptotic macrophages, as expected. To further investigate the role of macrophages on the spread of infection, we varied two macrophage parameters: the macrophage virus internalisation rate, and the macrophage activation half-max. We did not observe increasing infection with increasing macrophage virus internalisation rate (Fig. 17) . Furthermore, we did not observe any significant changes in the levels of infectious, apoptotic and removed epithelial cells (Fig.  18) ; in the levels of intracellular and extracellular viral loads per grid (Fig. 19) ; or in the cytokine levels (Fig. 20) . We argue that this is unsurprising, due to the high viral loads that are observed (e.g. Fig.  19(b) ); however, it does potentially suggest that the macrophage virus internalisation rate is either too low or too high to have an observable impact, alternatively it could suggest that the simulation is not sensitive to this parameter, but this cannot be determined without sensitivity analysis, which is beyond the scope of this article. Interestingly, we do observe a difference in resting macrophage recruitment (Fig. 18(d) ), even though there are no observable differences in cytokine levels; we attribute this to the stochasticity of the recruitment procedure. On the other hand, for the largest macrophage activation half-max value, we observe a significant influx of resting macrophages. An influx of macrophages/monocytes has been frequently observed in COVID-19 (see e.g. [6, 1, 40, 4] ). In our model, we expect the large recruitment of resting macrophages to be a consequence of larger infection and higher IL-6 levels; however, for the largest macrophage activation half-max value, we observe reduced infection ( Fig. 22(a) ), as well as lower IL-6 levels ( Fig. 24(b) ), when compared to the smallest macrophage activation half-max value. We do observe a slight increase in generic mononuclear phagocyte chemokine, which could (in combination with IL-6) be enhancing the recruitment of resting macrophages, suggesting a possible feedback between apoptotic cells (epithelial and macrophage) and the recruitment of macrophages, driving tissue injury and inflammation. However, this feedback is likely to be enhanced in our model because resting macrophages do not produce cytokines.

The presented model has several limitations: (1) we assume only a single immune cell population, with only a single phenotype; (2) we assume that macrophages may only activate via a single (classical) route; (3) healthy epithelial cells and resting macrophages do not produce cytokines; (4) only a limited set of cytokines is considered; and finally (5) there is a fair amount of uncertainty surrounding the parameter values. Mathematical and computational models require initial simplification, as well as subsequent iterative development, but can never-the-less, offer a unique perspective that can help enhance understanding. Therefore, the presented model serves as an initial starting point from which we can build, by including further complexity and more detailed study which can inform or enhance experimental studies.

When responding to a viral infection, no component of the immune system is expected to work by itself. However, mathematical and computational studies, as well as experimental studies, necessarily require reductions and simplifications, as well as iterative development. Indeed, in this article, we chose to focus only a single type and phenotype of macrophages, as well as a limited set of cytokines and chemokines. Therefore, extending the model to include multiple immune cell populations, as well as subpopulations, and a larger set of cytokines and chemokines is a subject of future work. Furthermore, we chose to focus on a single (classical) activation route, whereas in reality macrophages can alternatively activate; therefore, extending the model to allow for multiple activation routes is also a subject of future work. Multiscale, individual-based models are inherently complex and contain many parameter values which, in general, cannot be expected to directly fit to other modelling work, let alone experimental data. Therefore, improving the parameter estimation and fitting to experimental data, either directly or indirectly via ODE models (see e.g. [70] ) is a subject of future work. The presented model has the capability to include more complex intracellular pathways. Therefore, in the future, we would like to study, using more detailed ODE models, the crosstalk between virus replication, IFN-I induction and signalling, as well as pro-inflammatory cytokine induction and signalling, and how the crosstalk impacts the spread of infection. Additionally, the recently discovered Omicron variant seems to display an altered tissue tropism [71, 72] , and therefore we would like to consider additional models for virus entry and replication that are specific to different tissue environments. Finally, we would like to take the model to larger scales, either through computational improvements (such as massive parallelisation) or through mathematical models which can bridge-the-gap between the scales. 

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19

Dysregulation of Immune Response in Patients With Coronavirus

Tissue-Specific Immunopathology in Fatal COVID-19

Inflammatory profiles across the spectrum of disease reveal a distinct role for GM-CSF in severe COVID-19

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

Complex Immune Dysregulation in COVID-19 Patients with Severe Respiratory Failure

Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications

Heightened Innate Immune Responses in the Respiratory Tract of COVID-19 Patients

Tissue-based SARS-CoV-2 detection in fatal COVID-19 infections: Sustained direct viral-induced damage is not necessary to drive disease progression

Association Between Administration of Systemic Corticosteroids and Mortality Among Critically Ill Patients With COVID-19: A Meta-analysis

Interleukin-6 Receptor Antagonists in Critically Ill Patients with Covid-19 -Preliminary report. preprint, medRxiv

Type I and III interferon responses in SARS-CoV-2 infection. Experimental & Molecular Medicine

Type I interferons and SARS-CoV-2: from cells to organisms

The type I interferonopathies: 10 years on

Autoantibodies against type I IFNs in patients with life-threatening COVID-19

Genetic mechanisms of critical illness in COVID-19

Inborn errors of type I IFN immunity in patients with life-threatening COVID-19

A single-cell atlas of the peripheral immune response in patients with severe COVID-19

HIV-1 Dynamics in Vivo: Virion Clearance Rate, Infected Cell Life-Span, and Viral Generation Time

COVID-19 virtual patient cohort suggests immune mechanisms driving disease outcomes

Mechanistic Modeling of SARS-CoV-2 and Other Infectious Diseases and the Effects of Therapeutics

Modeling Within-Host Dynamics of SARS-CoV-2 Infection: A Case Study in Ferrets

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

Kinetics of Influenza A Virus Infection in Humans

An Age-based Multiscale Mathematical Model of the Hepatitis C Virus Lifecycle During Infection and Therapy: Including Translation and Replication

Validated models of immune response to virus infection

Host-pathogen kinetics during influenza infection and coinfection: insights from predictive modeling

Translational approaches to treating dynamical diseases through in silico clinical trials

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

Spatiotemporal Dynamics of Virus Infection Spreading in Tissues

A simple cellular automaton model for influenza A viral infections

Modelling the effects of bacterial cell state and spatial location on tuberculosis treatment: Insights from a hybrid multiscale cellular automaton model

NF-kB Signaling Dynamics Play a Key Role in Infection Control in Tuberculosis

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

Iterative community-driven development of a SARS-CoV-2 tissue simulator. preprint, bioRxiv

Comparative Computational Modeling of the Bat and Human Immune Response to Viral Infection with the Comparative Biology Immune Agent Based Model

Spatially distributed infection increases viral load in a computational model of SARS-CoV-2 lung infection. preprint, bioRxiv

A hybrid PDE-ABM model for viral dynamics with application to SARS-CoV-2 and influenza

Virological assessment of hospitalized patients with COVID-2019

Cytokine storm induced by SARS-CoV-2

Type I and Type III Interferons -Induction, Signaling, Evasion, and Application to Combat COVID-19

Should we stimulate or suppress immune responses in COVID-19? Cytokine and anti-cytokine interventions

Dysregulated Type I Interferon and Inflammatory Monocyte-Macrophage Responses Cause Lethal Pneumonia in SARS-CoV-Infected Mice

IFN-I response timing relative to virus replication determines MERS coronavirus infection outcomes

IL-6 in Inflammation

IL-10: The Master Regulator of Immunity to Infection

Coronavirus biology and replication: implications for SARS-CoV-2

Tissue distribution of ACE2 protein, the functional receptor for SARS coronavirus. A first step in understanding SARS pathogenesis

A Multibasic Cleavage Site in the Spike Protein of SARS-CoV-2 Is Essential for Infection of Human Lung Cells

SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor

Cell entry mechanisms of SARS-CoV-2

Coronaviruses use Discontinuous Extension for Synthesis of Subgenome-Length Negative Strands

SARS-CoV-2 structure and replication characterized by in situ cryo-electron tomography. preprint, bioRxiv

Comparing antiviral strategies against COVID-19 via multi-scale within host modelling

Intracellular Life Cycle Kinetics of SARS-CoV-2 Predicted Using Mathematical Modelling

GNU scientific library reference manual: for GSL version 1.12. Network Theory, Bristol

Post-mortem dissection of COVID-19: a pathogenic role for macrophages? Intensive Care Medicine

Monocytes and macrophages in severe COVID-19 -friend, foe or both?

Immunology of COVID-19: Current State of the Science

Alveolar macrophages: plasticity in a tissue-specific context

Differential immune activation profile of SARS-CoV-2 and SARS-CoV infection in human lung and intestinal cells: Implications for treatment with IFN-beta and IFN inducer

Modeling COVID-19 with Human Pluripotent Stem Cell-Derived Cells Reveals Synergistic Effects of Anti-inflammatory Macrophages with ACE2 Inhibition Against SARS-CoV-2. preprint, Research Square

Human Lung Stem Cell-Based Alveolospheres Provide Insights into SARS-CoV-2-Mediated Interferon Responses and Pneumocyte Dysfunction

Three-Dimensional Human Alveolar Stem Cell Culture Models Reveal Infection Response to SARS-CoV-2

Susceptibility of ferrets, cats, dogs, and other domesticated animals to SARS-coronavirus 2

Evasion by Stealth: Inefficient Immune Activation Underlies Poor T Cell Response and Severe Disease in SARS-CoV-Infected Mice

Virology, transmission, and pathogenesis of SARS-CoV-2. BMJ :m3862

SARS-CoV-2, SARS-CoV, and MERS-CoV viral load dynamics, duration of viral shedding, and infectiousness: a systematic review and meta-analysis

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

SARS-CoV-2 Omicron spike mediated immune escape and tropism shift. preprint, bioRxiv

SARS-CoV-2 Omicron efficiently infects human airway, but not alveolar epithelium. preprint, bioRxiv

). Note that the MOI = 0.01 and the other baseline parameters are given in Tables 1-4. investigate the IFN-I and IL-6 relationship more closely, but this is beyond the scope of the current article. Interestingly, for the longest IFN-I secretion delay, we did not observe a local maximum in the intracellular viral load per grid ( Fig. 11(a) ), further suggesting that IFN-I is responsible for the observed maximum, emphasising its key role in antiviral responses. Moreover, following the local maximum in the intracellular viral load per grid, we observed a rebound from a local minimum (Figs. 7(a), 11(a) ), which is commonly observed in low MOI studies and, in our model, is believed to be caused by viral replication inside epithelial cells with low extracellular IFN-I levels. To investigate this further, we studied the role of increasing IFN-I secretion rate on the spread of infection and observed a clear decrease in the spread