A mathematical model describing the localization and spread of influenza A virus infection within the human respiratory tract

Within the human respiratory tract (HRT), virus diffuses through the periciliary fluid (PCF) bathing the epithelium. But virus also undergoes advection: as the mucus layer sitting atop the PCF is pushed along by the ciliated cell’s beating cilia, the PCF and its virus content are also pushed along, upwards towards the nose and mouth. While many mathematical models (MMs) have described the course of influenza A virus (IAV) infections in vivo, none have considered the impact of both diffusion and advection on the kinetics and localization of the infection. The MM herein represents the HRT as a one-dimensional track extending from the nose down towards the lower HRT, wherein stationary cells interact with IAV which moves within (diffusion) and along with (advection) the PCF. Diffusion was found to be negligible in the presence of advection which effectively sweeps away IAV, preventing infection from disseminating below the depth at which virus first deposits. Higher virus production rates (10-fold) are required at higher advection speeds (40 μm/s) to maintain equivalent infection severity and timing. Because virus is entrained upwards, upper parts of the HRT see more virus than lower parts. As such, infection peaks and resolves faster in the upper than in the lower HRT, making it appear as though infection progresses from the upper towards the lower HRT, as reported in mice. When the spatial MM is expanded to include cellular regeneration and an immune response, it reproduces tissue damage levels reported in patients. It also captures the kinetics of seasonal and avian IAV infections, via parameter changes consistent with reported differences between these strains, enabling comparison of their treatment with antivirals. This new MM offers a convenient and unique platform from which to study the localization and spread of respiratory viral infections within the HRT.


Introduction
Mathematical models (MMs) of influenza A virus (IAV) infection kinetics are mainly based on ordinary differential equations (ODEs) that describe the time evolution of infection (virus as a function of time) with the implicit assumption that all cells see all virus and vice-versa, as discussed in several reviews [10,23,29]. A different approach from ODE MMs is to use agent-based MMs which treat the dynamics of each cell, and even each virus, individually and track local cell-virus interactions [5,7,71]. Such MMs are computationally intensive, restricting the number of cells, and therefore the area, which can be represented, and often lack support or validation in the form of experimental localization data [4,7,50]. Infection localization can also be modelled by dividing the human respiratory tract (HRT) into compartments corresponding to different regions of the HRT, where each compartment has different parameters based upon physiological differences [62].
IAV infection spread faces a tremendous physical, spatial barrier in the form of the thick mucus layer which lines the HRT, trapping and carrying free virions rapidly upwards towards the nose and throat and out of the HRT [26]. This physical clearing process is typically taken somewhat into account in non-spatial (1) represents the HRT as a one-dimensional track, as illustrated. The MM considers only the virus that is located in, and diffusing within the PCF, which also moves along with the PCF upwards at a fixed advection speed. What happens to the virus when it reaches the bottom (reflective) or top (absorptive) of the MM-represented HRT is indicated. MM (1) describes the interactions between the moving virus concentration in the PCF, and the stationary cells that carpet the HRT.
ODE MMs via a term for the exponential clearance of virions. This simplification, however, does not account for the fact that not all cells will have equal exposure to the virus: cells that are downstream of the mucus flow (higher in the HRT) will have greater exposure to virus than those upstream, located deeper within the HRT. To our knowledge, the effect of virus advection due to entrainment by the mucus layer on IAV infection localization and spread within the HRT has never been evaluated.
In this work, a spatiotemporal MM for the spread of IAV infection within the HRT using partial differential equations (PDEs) is constructed. Through its one-dimensional representation of the HRT, the MM is used to study the effect of viral transport modes on the course of an IAV infection in vivo. Via the addition of only two additional spatial parameters whose value is relatively well-established, namely the rate of diffusion and advection of virions, the MM able to produce a richer range of IAV kinetics. The effect of cellular regeneration and a simplified immune response are also considered. The PDE MM is also used to explore differences in the kinetics of IAV infection, in the absence and presence of antiviral therapy, for IAV infection with either a seasonal or an avian-adapted strain, with the latter being typically more severe [56].

A spatial MM for IAV in the HRT
In the proposed spatial MM, the HRT is represented as a one-dimensional tract running along the x-axis indicating the depth within the HRT, with x = 0 cm located at the top of the HRT (nose), and x = 30 cm terminating somewhere within the bronchi [60,70], as illustrated in Figure 1. It is an extension of the 'standard' MM for IAV in vitro [55,59,68] which adds: (1) the diffusion of virions through the periciliary fluid (PCF) which lies between the cells' apical surface and the thick mucus blanket which lines the airways; (2) their advection due to upwards entrainment of the PCF by the mucus layer; and (3) their effect on the one-dimensional, depth-dependent fraction of non-motile (stationary) cells in various stages of infection. The spatial MM is formulated as The fraction of uninfected, target cells, T (x, t), located at a depth x at time t are infected at rate βV (x, t), proportional to the concentration of virions, V (x, t), at that depth and time. Newly infected cells are in the eclipse phase, E i (x, t), i.e. they are infected but are not yet producing virus. After an average time τ E , infected cells leave the eclipse phase to become infectious cells, I j (x, t), producing virions at constant rate p for an average time τ I until they cease virus production and undergo apoptosis. As in [55,68], the eclipse and infectious phases are each divided into n E = n I = 60 age classes so that the time spent by cells in each phase follows a normal-like distribution, consistent with biological observations [36].
The MM assumes virions are released from stationary infected cells into the PCF, so V (x, t) is the concentration of virus in the PCF. The produced virions diffuse through the PCF at rate D PCF , and are transported upwards (towards the nose) within the PCF at speed v a , as the PCF is entrained by the mucus blanket pushed along by the beating cilia of the ciliated cells that line the HRT [47,69]. When the virus in the PCF moves beyond the bottom or top of the HRT, it either 'bounces back' (reflective) or is lost (absorptive), respectively. Absorption of virions into the mucus blanket (out of the PCF), their loss of viral infectivity over time (in the PCF), and other modes of non-specific virion clearance are all taken into account via a single exponential viral clearance rate term, cV (x, t). Virions contributing to the infection in this MM are those remaining in the PCF, in direct contact with infectable cells.
In the spatial MM simulator, infection is initiated by a spatially localized, spray-like virus inoculum deposited at depth x d within the HRT. The spray-like inoculum is represented by a Gaussian centred at the site of deposition, with a standard deviation of 0.5 mm, about 10× the size of a large cough droplet [80]. At sites far from x = x d , V (x, t = 0) ≈ 0. The baseline values of the spatial MM's parameters were all set to values estimated in Baccam et al. [1], obtained by fitting a non-spatial MM to patient data from experimental primary infections with the influenza A/Hong Kong/123/77 (H1N1) virus. Table 1 lists the initial conditions and parameters used. A complete description of the spatial MM is provided in the Methods.

IAV kinetics in the presence of virus diffusion and advection
Infection of cells by the IAV, and the subsequent release of virions, occurs almost exclusively at the apical surface of the epithelium [15,51]. When advective motion of the PCF is neglected, virions can be assumed to undergo Brownian motion and the virus concentration distribution in the PCF will be governed by diffusion. Figure 2(a-c) shows IAV infection kinetics in the presence of diffusion alone, where quantities are shown averaged over all spatial sites as a function of time (see Methods). The infection takes place at a much slower pace in the spatial (diffusion only) MM than in the non-spatial MM which implicitly assumes an infinite diffusion coefficient. Figure 2(d-f) shows the spatial extent of the infection dissemination (with diffusion only) through the HRT at certain times post-infection. Virus becomes available to target cells gradually as it diffuses away from the site of initial infection, x d = 15 cm, and the infection wavefront moves outwards initial fraction of uninfected target cells * These values are used in all simulations unless otherwise stated. The spatially averaged initial inoculum, V (x, t = 0) , equals the value in [1] (see Methods). Values for τ E and τ I were chosen so as to lie near the middle of the ranges of [6,10] h and [10,40] h, respectively, obtained from MMs of IAV infections in vitro [55,68]. The value for p is discussed in Section 2.2, and that for p Baccam is explained in Methods.
from that site symmetrically towards the two ends of the HRT at x = 0 and x max . S1 Video depicts the spatiotemporal course of the infection, in the presence of diffusion only, in the form of a video. Along the length of the HRT, a layer of mucus about 0.5 µm-5 µm thick covers the PCF [26]. The collective motion of the underlying epithelial cells' beating cilia, dubbed the mucociliary escalator, drives this mucus layer upwards. It also leads to an upward advection of the PCF, at a speed similar to that of the mucus layer [47], entraining any virus in the PCF upwards at that speed. Given the advection speed of the mucus and PCF (v a ≈ 40 µm/s), any newly produced or deposited virion would be cleared from the HRT in less than ∼12 min (30 cm/v a ). Thus, it is not surprising that adding advection to the spatial MM, while still using the same parameter values, results in a subdued, low viral titer, slow-growing infection which still has not peaked by 7 dpi. In contrast, viral titer in patients infected with influenza A/Hong Kong/123/77 (H1N1) virus in Baccam et al. [1] peaks 2-3 dpi, which the non-spatial MM with these same parameters reproduces well. Figure 3(a-c) shows the IAV infection kinetics in the presence of both diffusion and upward advection as the virus production rate, p, is increased from the base value estimated by Baccam et al. [1] for their nonspatial MM, p Baccam (see Methods). Increasing the virus production rate to 11 × p Baccam yields an infection that peaks at ∼ 3 dpi. This adjusted value for the virus production rate (p = 11 × p Baccam , see Table 1) is used in the remainder of this work so that the spatial MM qualitatively reproduces the approximate timing of viral titer peak in an IAV infection.
Using the new value for the virus production rate, Figure 3(d-f) shows the spatial extent of the infection dissemination in the presence of both diffusion and advection in the spatial MM. It illustrates the protective effect of advection: preventing the infection from travelling much beyond its initial deposition depth (x d = 15 cm), as seen from the target cell depletion shown in Figure 3(d). Figure 3(g-i) explores the effect of varying this depth of deposition of the initial virus inoculum, x d . When the inoculum deposits lower in the HRT (as x d increases), the fraction of HRT consumed by the infection increases, resulting in higher viral titer peak and total virus yield. However, the virus titer curve is largely insensitive to the depth of deposition for   In the presence of both diffusion and advection in the spatial MM, the initial virus inoculum travels quickly from its deposition site at a fixed speed up the HRT, leaving in its wake a small, equal fraction of infected cells at all sites above the deposition point (∀x < x d ). As the newly infected cells begin to release virus which is entrained up the HRT, uninfected cells are only exposed to virus produced by cells lower than (upstream of) their own depth (greater x value), such that cells higher in the HRT (smaller x value) are exposed to the most virus. Consequently, infection in the upper HRT proceeds much faster than that in the lower HRT. Figure 3(j-l) shows the localized fraction of target and infectious cells, and infectious virus concentration over time at specific sites or depths along the HRT. Despite infection kinetic parameters being the same at all sites, infection progresses faster at sites located higher in the HRT, downstream of the virus flow.
Interestingly, the more rapid progression of the infection higher in the HRT makes it appear as though the infection is moving downwards, from the top (nose, x = 0) towards the bottom (x = 30 cm) of the HRT. Both the target cell depletion wavefront in Figure 3(d) and the infected cell "pulse" in Figure 3(e) appear to be moving to the right (towards the lower HRT) as time advances. This can be seen more clearly in the video (S2 Video) which depicts the spatiotemporal course of the infection, in the presence of both diffusion and advection. Cells in the upper HRT are consumed first, after which the pulse of infected cells appears to slowly move backwards, down the HRT, back to the initial deposition site, where it then slowly dissipates. Figure 4(a) provides a more "zoomed-in" view of what happens to the pulse in Figure 3(e) at times beyond 7 dpi. The pulse is shrinking and slowing down as it approaches the initial inoculum deposition depth of x d = 15 cm. Figure 4(b) quantitatively displays the (nearly exponential) decay of both the amplitude and speed of the pulse, as its peak's position (µ p ) approaches its final, asymptotic position of ∼ 14.1 cm, a little short of its initial deposition depth. Figure 4(c) shows the difference between the pulses obtained in the presence of both diffusion and advection minus that in the presence of advection alone (no diffusion). While advection dominates over diffusion, the latter still plays a role in allowing a small amount of virus to travel upstream, against the advection, which results in slightly more cells infected upstream and slightly less downstream. The effect of diffusion in the presence of advection is most likely negligible (one part in 10 5 in the fraction of cells infected) in light of typical, in vivo experimental variability (∼10-fold variability).

Target cell regeneration during IAV infection
During an uncomplicated IAV infection, viral loads typically fall below detectable levels between 6 dpi-8 dpi [2,13]. An experimental study in which mice were sacrificed at various times over the course of an IAV infection [61] found that by 3 dpi, ciliated cells were rarely observed on the tracheal surface, and instead the latter was mostly composed of a layer of undifferentiated cells. By 5 dpi, signs of repair were apparent, by 10 dpi the tracheal surface was covered with ciliated and nonciliated cells, and by 14 dpi the surface was indistinguishable from the uninfected surface. Another experimental study in hamsters suggests that 6 d-7 d after mechanical injury, most of the epithelium is comparable to that of control hamsters [42]. A similar study in guinea-pigs suggests that 5 d after mechanical injury, the epithelium was ciliated and differentiated, and by 15 d it was indistinguishable from that of control guinea-pigs [25]. Together, these studies suggest that the process of cellular regeneration following damage caused by an IAV infection would at least begin, if not be well underway, before the IAV infection has fully resolved.
Following airway injury, numerous factors are thought to trigger the start of the repair process [18]. Once triggered, it is believed neighbouring epithelial cells will try to stretch to cover the denuded area, and divide so as to regain their normal shape while maintaining coverage. As damage becomes more significant, progenitor cells, mainly basal cells exposed to the PCF as the epithelial cells above them detach or are removed by apoptosis or damage, also undergo proliferation and differentiation until the epithelium is restored to its pre-injury state [64]. With these processes in mind, cellular regeneration was implemented in the spatial MM by replacing the equation for uninfected, target cells in MM Eqn. (1) with where is the fraction of dead cells, i.e. the extent of the damage or injury. As such, cellular regeneration in Eqn.
(2) proceeds at a rate that is greater in the presence of greater damage (higher D) and of greater fraction of cells available to regenerate (higher T ). Parameter r D sets the scale of the regeneration rate (which also depends on T and D), and τ D is the regeneration delay such that the current regeneration rate depends on the fraction of target cells (T ) currently available to repopulate the area, and on the amount of damage (D) that was perceived some time τ D ago. The delay between damage and regeneration accounts for the time required for the damage to activate appropriate signalling pathways, and for both division and differentiation to take place so that newly regenerated cells are susceptible to infection. This same equation has been used by others to represent target cell regeneration during an IAV infection [12,28], but the authors therein did not include a delay (τ D = 0). An experimental study of cell regeneration in hamsters [42] suggests that in many cases, between 18 h-24 h following a mechanical injury, cells have already stretched or migrated to cover the complete denuded area and cell division has begun. A similar study in guinea-pigs [25] also suggests that at 15 h after mechanical injury, cells have migrated to cover the denuded area and cell proliferation is underway. This suggests that the delay between damage and the start of cellular division is around 15 h-24 h. Herein, this delaynamely the time elapsed between damage and the start of both cellular division and differentiation to an extent that newly regenerated cells are susceptible to infection -was chosen to be τ D = 1 d, based on  these experimental studies. To select an appropriate value for the regeneration rate (r D ), HRT epithelial cell regeneration following mechanical injury was simulated using Eqn. (2). Results are shown in Figure 5 for different regeneration rates and delays: the regeneration rate determines the steepness of the regeneration, while the regeneration delay sets its timing. These parameters should be easily identifiable if experimental data was available. An intermediate value of r D = 0.75 d −1 was chosen to ensure that with a delay of 1 d, regeneration is well underway by 5 d-8 d, and completely resolved by 12 d-14 d [25,42,61]. Figure 6 shows the infection kinetics in the presence of cellular regeneration. At this stage, and over the range of values considered here, the spatial MM is insensitive to the choice of regeneration rate or its delay. This is largely due to the fact that, in the absence of an immune response, the HRT above the inoculum deposition point is decimated by the infection at a rate much faster than the density-dependent regeneration can counter. As damage increases, the number of target cells available to replenish the lost cells also decreases dramatically, preventing any significant regeneration. To get a more realistic view of the role and impact of cellular regeneration on infection kinetics, the protective role of the immune response must be considered.

Immune response to an IAV infection
Past MM for IAV infection have incorporated one or more immune response components with varying degrees of success [11,12,23]. Difficulties in incorporating an immune response when modelling IAV infections include the complex nature of the interactions (large networks of cells and signals), and the lack of appropriate data (quantity and quality) to inform the MMs [29]. The simplified immune response considered herein comprises an innate response based on interferon (IFN), a humoral response represented by antibodies (Abs), and a cellular response embodied by cytotoxic T lymphocytes (CTLs). Since our aim is to display the MM's range of kinetics -rather than identify parameters, analyze data, or challenge hypotheses -a parsimonious approach was preferred wherein key immune components are represented using empirical curves that broadly reproduce the scale and timing of these experimentally measured quantities (see Methods). Figure 7 shows the empirical MM curves against their corresponding experimental time courses, for IFN, Abs, and CTLs, during experimental, in vivo IAV infections.
Although IFN is known to have many effects [40], herein its effect is to reduce the virus production rate p [1,30], via a resistance parameter, f 50 , analogous to the IC 50 used to describe antiviral resistance. The resistance is defined such that if f 50 = 0.8, the virus production rate is halved when IFN concentration is  [41,49]; and (c) cytotoxic T lymphocytes (CTLs) [48,49], taken over the course of in vivo IAV infections in mice, unless otherwise specified.  [38,66], (b,e) the Ab response [41] or (c,f) the CTL response [43,77,81] were disabled. In these experimental studies, the animal model is mice, except Seo et al. [66] which was conducted in pigs. Experimental viral load is measured in TCID 50 /mL for [66], in TCID 50 /mouse (of lung homogenate) for [38], in EID 50 /mL for [81], [43] and for [77] and in pfu/mL for [41]. The MM results are shown with f 50 = 0.5, k A = 500 h −1 and k C = 50 h −1 , or disabled with F (t) = 0, A(t) = 0 or C(t) = 0 in (d,e,f), respectively.
80% of its peak value (see Methods). Figure 8(a-c) shows the effect of IFN in the MM as resistance to IFN, f 50 , is lowered (increased sensitivity). The initial viral titer peak is reduced due to IFN presence, with lower resistance (smaller f 50 ) having a greater impact. However, once IFN starts to decay, its effect rapidly dissipates, leading to a rise, or even a rebound, in the viral load. On its own, IFN does not lead to infection resolution in this MM. The effect of Abs herein, like elsewhere [9,28,30,44,49], is to enhance infectious virus clearance such that c becomes c+k A A(t), where k A represents the neutralization rate of infectious virus by Abs, A(t). Figure 8(df) shows the effect of IFN+Abs in the MM as the neutralization rate, k A , is increased. Low binding affinity Abs (k A ≤ 200 h −1 ) cannot clear the infection. As k A increases, viral titer decay rates increase, leading to infection resolution within 11 dpi-16 dpi. However, the time to infection resolution remains longer than the 6 dpi-8 dpi typically observed for IAV infections in humans [2,13].
The effect of CTLs herein, like elsewhere [9,28,44,49], is to increase the rate of loss of infected cells expressing IAV peptides on their MHC-1. Infected cells begin expressing IAV peptides ∼ 4 h post-infection [9,78], or approximately halfway through the eclipse phase (8 h in the spatial MM, see Table 1), such that all infected cells past the mid-point of their eclipse phase will be removed at rate k C C(t), where k C represents the killing rate of infected cells by CTLs, C(t). Figure 8(g-i) shows the effect of IFN+Abs+CTLs in the MM as the killing rate, k C , is increased. Higher CTL killing rates (k C ) cause the infection to resolve earlier. Table 2: Parameter values used to reproduce the IAV infection kinetics in Figure 10.   [38,66] at early time as both suggest IFN acts early and helps to reduce the initial viral titer peak. The experimental Ab knockout experiment [41] suggests Abs help reduce viral load at a later time in the infection but does not affect the initial viral titer peak. The MM prediction also suggests Abs help reduce the viral load at a later time but also affects the initial viral titer peak. This could be better reproduced by the MM by having a lower initial amount of Abs (A 0 ) so that Abs appear later and thus act later. Finally, the MM prediction of the CTL knockout experiment is in good agreement with the experimental data of CTL knockout experiments [43,77,81] as both suggest that CTLs act only at a late stage in the infection, helping reduce infection duration. Generally, the experimental knockout appear to show a greater impact for the knockouts than predicted by the MM. This is likely due to the intricate interaction network between the different immune response components in vivo which makes it difficult to experimentally disrupt one component without affecting another.

Capturing the kinetics of in vivo infections with seasonal and avian IAV strains
The difference in infection severity between patients naturally infected with seasonal IAV strain or avianorigin H5N1 strains is thought to be the result of a number of different possible factors, including more rapid infection and replication kinetics, cell tropism, lack of pre-existing immunity, and/or an aberrant immune response [14, 16, 19, 39, 66-68, 72, 74, 75]. Figure 10(a) shows viral load measurements from pharyngeal swabs of patients (each point is a single measure from a single patient) naturally infected with either a seasonal (H3N2 or H1N1) IAV strain or a strain of the H5N1 subtype. This study wherein these measures were taken over the same time period and following the same methodology, such that viral loads for infection with seasonal and H5N1 strains can be readily compared, is the only one of its kind. Since patients were recruited into the study days after illness onset, and began receiving antiviral therapy immediately upon hospital admission, untreated infection kinetics leading up to and after admission are not available. The measured pharyngeal viral load (determined via qRT-PCR) is ∼100-fold higher for patients infected with H5N1 than seasonal IAV strains. This could also be true of the infectious viral titer (typically measured via TCID 50 or PFU) which was not measured as part of this study, although the ratio of infectious to total IAV is known to change over the course of an infection, and inconsistently so between experiments [58]. Another feature is that H5N1-infected patients reported to the hospital 3 d-4 d later after illness onset than those infected with seasonal IAV strains. It is unclear whether this is due to a slower progression of infection with an H5N1 strain, or because patients infected with H5N1 strains mostly came from remote provinces and took longer to reach the city hospital than those infected with seasonal strains which came from the city itself or neighbouring provinces. From the data in Figure 10(a), two hypothetical infection time courses or portraits were developed for infection with IAV H5N1 strains, and one for seasonal strains, shown in Figure 10(b). The portrait for infection with a seasonal strain peaks at ∼ 10 3 -10 4 TCID 50 /mL at ∼ 3 dpi, and resolves by ∼ 8 dpi [34,63]. For infection with a strain of the H5N1 subtype, two options are considered: the infection peaks early, at higher titers which are sustained over a longer period; or the infection grows at a rate similar to infection  [19]). The shaded polygons trace out the extent of the data points. (b) Proposed time courses for in vivo IAV infections with either a seasonal (grey) or H5N1 (black) IAV strain. The solid and dashed black lines represent two plausible time courses for infection with an IAV H5N1 strain: a more rapid rise to higher viral titers than for infections with seasonal strains (solid); or a more moderate rise, similar to that seen for infections with seasonal strains, that grows to higher titer due to lack of rapid and effective immune control (dashed). MM parameters for these curves are listed in Table 2. (c) The possible effect of a delayed CTL response, i.e. one which peaks at 12 dpi with a killing efficacy of k C = 50 h −1 (grey) in both the early (solid) and late (dashed) infection time courses which are otherwise shown without a CTL response (black). with a seasonal strain, but that growth is sustained over a longer period and thus peaks later, reaching higher titers than infection with a seasonal strain. For infection with an H5N1 strain, viral titer peak was chosen to be 100-fold higher, as observed for total viral loads, and the reported 3 dpi-4 dpi delay in hospital admission was captured either as a longer, sustained infection in the early time course, or as a longer delay to reach peak titer in the late time course. The IAV H5N1 portraits were constructed from the MM of the seasonal portrait by shifting parameters controlling cell-virus interactions and immunity in keeping with differences believed to be responsible for that shift. These parameters are presented in Table 2. The shaded areas, which depict the extent of the data points in Figure 10(a) where time is based on days since illness onset, are shown time-shifted in Figure 10(b) where time is days post-infection. For a seasonal infection, the (pale grey) area was shifted to one day later since there is generally a delay of one day between infection and symptom onsets [13]. For avian infections, the (dark grey) area was shifted to 4 days later since reports suggest ∼ 3-5 days elapse between symptom onset and known potential exposure to infected poultry [17,35,53].
For cell-virus interactions, the early time course relies on increased virus production (p), increased virus infectivity (β), and a shorter eclipse phase (τ E ), consistent with shifts in these parameters estimated from in vitro infections of A549 cells with seasonal (H1N1) versus H5N1 and H7N9 IAV strains [68]. In contrast, the late time course relies on lower virus infectivity and a longer eclipse phase, along with a higher virus production rate. The early and late portraits also depend on 8-and 10-fold longer periods of virus production (infectious cell lifespan, τ I ), respectively, compared to that for cells infected with a seasonal strain. In capturing differences in immunity, both the early and late portraits of infection with IAV H5N1 strains rely on 10-fold increased resistance to the effect of IFN (f 50 ), decrease in pre-existing immunity in the form of a 200-fold decrease in the neutralizing effect of Abs, and the absence of a CTL response. The resistance of H5N1 strains to the effect of IFN-α and -γ has been reported in vitro in porcine cell cultures and in vivo in pigs [66]. The decrease in the neutralizing efficacy of Abs is captured in the simplistic Ab MM utilized herein as a 200-fold decrease in their initial quantity (A 0 ) which means a longer time (∼8 d longer) to reach maximal Ab activity, but can equally be captured by a 200-fold decrease in the neutralization efficacy (k A ) of Abs. A recruitment study of patients naturally infected with avian H7N9 IAV strains reported that the patient group with the earliest recovery, had a prominent CTL response by 10 d post admission and a prominent Ab response, 2-3 d later [76]. In contrast, the patient group with the most delayed recovery had an Ab and CTL response that remained low even after 20 d post admission. This suggests that an infection with an H5N1 strain, characterized by a prolonged viral shedding period, would have delayed Ab and CTL responses. Because the portraits in Figure 10(b) only go up to 12 dpi, the effect of these delayed CTLs would not be apparent. Figure 10(c) illustrates what the time course of infection with an avian strain could look like in the presence of a delayed CTL response, but one which would be as effective (same k C ) as that for infection with a seasonal strain. It shows an important role for CTLs in controlling the persistent shedding in these infections, consistent with the findings that a delayed CTL response correlated with slower recovery from infection with H7N9 strains [76]. Figure 11 explores the impact of antiviral therapy with neuraminidase inhibitors (NAIs), administered at various times post-infection, on the time courses for infection with seasonal or avian IAV strains. Since NAIs block the release of newly produced virions, their effect is implemented in the MM, as elsewhere [8,11,[20][21][22]31,45,54], as a reduction in the virus production rate, namely (1−ε NAI )p, from the time treatment is administered, where ε NAI ∈ [0, 1] is the drug efficacy. The efficacy of NAIs was set to ε NAI = 0.98 to study the case of treatment with a high efficacy. Table 3 quantitatively compares the impact of NAI treatment for various endpoints: reduction in resolution time, peak viral titer, and area under the viral titer curve (AUC).
Human volunteer studies in patients experimentally infected with seasonal (H1N1) IAV strains and treated with NAIs report that early treatment (24-32 h post infection) is effective in reducing viral load [3,32,33], and one such study reports that delayed treatment (50 h post infection) is less effective [33], in line with the MM results. A study of H5N1-infected patients recruited days after illness onset, and beginning NAI treatment immediately after admission, reports that late treatment (4-8 d after illness onset) is still effective in some patients, reducing viral load and thus possibly also reducing time to infection resolution. This would be consistent with some avian strain-infected patients experiencing an infection time course similar to the early portrait with moderate to no benefit from delayed NAI treatment, and others experiencing infection  Table 2.  ND − 16 = ND 10 7.3 /10 6.8 = 10 0.5 10 6.8 /10 6.7 = 10 0.1 * For NAI therapy administered at various times post infection (t adm ), the change (∆) in 3 commonly reported endpoints are shown for the viral titer curves in Figure 11. The change in resolution time is computed as a difference (∆ = without − with), whereas that for the area under the viral titer curve (AUC) and viral titer peak corresponds to the fold-change (∆ = without/with). ND stands for not determined in cases where infection resolution (V (t) < 1 TCID 50 /mL) was not achieved without treatment. more consistent with the late portrait and benefiting from NAI therapy even when treatment initiation is delayed.

Discussion
Mathematical models (MMs) for the course of an influenza A virus (IAV) infection in vivo typically assume the infection is spatially homogeneous, i.e. that all cells see all virus, and vice-versa, instantly over all space. This simplification is reasonable in vitro in cell culture infections whose spatial extent is often no more than one or two square centimetres [37], but it is unclear whether it remains appropriate at the scale of the entire HRT. With simplicity and parsimony in mind, the spatial MM introduced herein represents the HRT as a one-dimensional track that extends from the nose down to a depth of 30 cm. It implements two modes of viral transport: advection of virus upwards towards the nose, and diffusion of the virus within the periciliary fluid.
When diffusion alone is considered, infection in the spatial MM proceeds at a slower pace than in the non-spatial MM, but virus production is sustained for longer. The diffusion initially acts to spatially restrict the number of cells available to the infection, and then releases these cells progressively as the diffusing infection wave reaches ever further in the HRT. When advection is added to diffusion, the former dominates the infection kinetics, requiring a 10-fold increase in the virus production rate to restore the timing and level of viral titer peak to that in the non-spatial MM. This is noteworthy, firstly, because it shows that advection constitutes an effective physiological mechanism to suppress infection. Secondly, it shows that use of a non-spatial MM to analyze infection data possibly underestimates the virus production rate, and consequently also the total amount of virus produced over the course of an infection. Since the more virus are produced, the more mutations accumulate [24,31,57], underestimating the virus production rate means underestimating the rate and likelihood of emergence of drug resistance.
Interestingly, the MM suggests that the depth at which the initial virus inoculum deposits also plays a major role in determining the extent of HRT involvement in the infection. Specifically, the MM predicts no target cell consumption below the depth at which the initial virus inoculum deposits. While deposition at depths lower than ∼ 15 cm did not affect the viral titer time course, initial depositions above that point both delayed and reduced peak titer. As such, higher virus production rates would be required to compensate for shallower deposition depths to obtain equivalent viral titer time courses. This is one more way a nonspatial MM might underestimate the rate of virus production, and thus the likelihood of antiviral resistance emergence.
A more unexpected finding was that although kinetically in the spatial MM the infection deposits at some depth and is dragged upwards by advection, the infection appears to be moving from the upper to the lower HRT. This is because the upper HRT is downstream of virus advection and sees most of the virus produced, while the lower HRT is upstream and sees little or none. Therefore, infection progresses, peaks and resolves faster in the upper HRT and slower in the lower HRT. It is this difference in time scale which makes it appear as though the infection is moving from the upper to the lower HRT. Experimental studies of mice infected with IAVs engineered to be fluorescent or bioluminescent observe infection spread from the upper towards the lower regions of the lung [46,73], consistent with the MM's predictions. In the absence of sufficient immune control, this could serve to maintain the infection with, wherein the slower, lower HRT infection could re-ignite more rounds of infection in the upper HRT, depending on the delicate timing between infection time course in the lower HRT and cell regeneration in the upper HRT.
The spatial MM was expanded to include density-dependent cell regeneration, and a simplistic immune response comprising IFN acting to down-regulate the rate of virus production by infected cells, Abs neutralizing infectious virions, and CTLs killing infected cells. In the presence of this immune response, the MM-simulated IAV infection was well controlled and resolved fully by 8 dpi. The spatial MM could also largely reproduce the key features of in vivo immune response knockout experiments [23]. In the MMsimulated infection (see Figure 8g), 10% of the HRT was involved in the infection by 3 dpi, around 6 dpi the fraction of cells involved in the infection peaked at 30%, and by 10 dpi 10% of the epithelium had yet to regenerate, although viral titers had fallen below the detection limit (see Figure 8i). These numbers align well with a 1989 report in Russian cited by Bocharov and Romanyukha in their 1994 IAV MM paper [9] of 10% damage at symptom onset, 30-50% of the upper airway destroyed at the peak of the disease, and resolution of disease at a time when as much as 10% of the normal epithelium is still damaged.
The full spatial MM was applied to simulate the kinetics of infection with either a mild, seasonal IAV strain or a severe infection with an avian strain such as H5N1. Since complete, untreated, infection kinetic time courses for infection with avian IAV strains are not available, two different portraits were constructed to represent plausible time courses: one peaking early with sustained, high titers over several days, and another rising slowly, similarly to seasonal infections, but continuing on to reach high titers 3-4 days later. Differences in the time courses for infection with a seasonal vs avian IAV strain could be captured by the MM by shifting the parameters in ways that are consistent with known or expected differences between infections with seasonal and avian strains. The portraits were used to evaluate the impact of antiviral therapy with neuraminidase inhibitors. Treatment was most effective for the seasonal and late peak avian IAV strain portraits when administered early during the infection at 2 dpi. Later treatment was of limited effectiveness for the seasonal IAV strain portrait, but was still effective in the late peaking portrait of infection with an avian IAV strain.
The spatial MM presented herein, despite its simplicity, offers a novel and interesting opportunity to study the spatial localization and dissemination of IAV infections within the HRT. In the future, the addition of an explicit, dynamical immune response -to replace the empirical equations used herein as stand-ins for immune response components -would enable the use of this MM to reproduce the kinetics of re-infection, similar to the work done by others [12,79].

Acknowledgements
This work was supported in part by Discovery Grant 355837-2013 (to CAAB) from the Natural Sciences and Engineering Research Council of Canada (www.nserc-crsng.gc.ca), Early Researcher Award ER13-09-040 (to CAAB) from the Ministry of Research and Innovation of the Government of Ontario (www.ontario.ca/page/early-researcher-awards), and by Interdisciplinary Theoretical and Mathematical Sciences (iTHEMS, ithems.riken.jp) at RIKEN (CAAB).

Financial disclosures
CAAB received financial support in the form of a research contract and a consultancy fee from F. Hoffmann-La Roche Ltd. in the early stages of this project. MBR was employed by F. Hoffmann-La Roche Ltd. when first participating in this work, and is currently employed by Array BioPharma Inc.. Neither company had any role in the study design, data collection and analysis, or decision to publish.

Numerical implementation of the mathematical model
The Euler method was used to numerically solve the cell population equations, namely I n+1 j,m = I n j,m + ∆t where T n m = T (x m , t n ), x m = m ∆x, t n = n ∆t, and ∆x and ∆t are the chosen spatial and temporal step sizes, respectively. We define N x = x max /(∆x) = 3000, the number of spatial boxes or sites that make up the simulated HRT such that ∆x = 0.3 m/(3000 sites) = 100 µm/site, the diameter of ∼4-5 epithelial cells. This N x was chosen by verifying that choosing a larger number of sites (smaller ∆x) did not affect the solution. When presenting the fraction of target cells or infectious cells, or the virus concentration, as a function of time only, averaged over space, those are computed as For the virus concentration, for simplicity, the diffusion and advection terms are each treated separately, and are applied before the production and clearance terms are considered. To solve the diffusion term, the Crank-Nicolson method was used, namely where α = D PCF (∆x) 2 /(∆t). The rate of viral diffusion in the PCF was estimated based on the Stokes-Einstein equation for IAV diffusing in plasma at body temperature as D PCF ≈ 10 −12 m 2 /s [6,37]. An absorbing boundary condition was used at the top of the HRT, V (0, t) = 0 (or V n 0 = 0), to capture virus flow out through the mouth and nose. A reflective boundary condition, V (x max + ∆x, t) = V (x max , t) (or V n Nx+1 = V n Nx ), was used at the bottom of the HRT. The bottom boundary condition becomes irrelevant once advection is introduced as the flow at x = x max becomes negligible. With these boundary conditions in place, the diffusion of virus over a step size ∆t can be expressed as The virus advection term can be written as follows This simplistic numerical scheme is known to lead to a dispersion (diffusion) of the solution [52]. However, the spatial MM simulator imposes v a ∆t/(∆x) = 1, such that Eqn. (4) simplifies to V n+1 m = V n m+1 , a trivial translation of the solution, which ensures the latter will remain stable and will not disperse. The speed of the PCF is believed to be approximately the same as that of the mucus blanket and is estimated to be v a ≈ 40 µm/s [47]. This requires the time step for the spatial MM simulator to be ∆t = ∆x/v a = (100 µm)/(40 µm/s) = 2.5 s. Here, the boundary condition at the top is irrelevant since the solution depends only on the downstream element (V n+1 1 = V n 2 ), and the boundary condition at the bottom is chosen such that V n Nx = 0, i.e. there is no virus beyond the end of the HRT. Once diffusion and advection have been applied, the Euler method is used to solve the remaining terms of the virus equation, namely The droplet-like, initial IAV distribution V (x, t = 0) is represented by a Gaussian centred at the site of deposition, where x d is the site of deposition, σ = 0.5 mm, and V 0 * is chosen so that V (x, t = 0) = 1

Nx
Nx m=1 V 0 m = 7.5 × 10 −2 TCID 50 /mL, i.e. the spatial average of the initial virus inoculum concentration in the spatial MM herein is equal to V (t = 0) in the non-spatial ODE MM by Baccam et al. [1].
While most parameters were taken directly from [1], some were adapted. Whereas in Baccam et al. [1] T + E + I = 4 × 10 8 corresponds to the number of cells, in the spatial MM herein, fraction of cells in each state are considered instead such that T + E + I = 1. As such, what we consider the Baccam et al. [1] virus production rate is converted from that reported in their paper as Once advection is introduced (see Results), a value of 11 × p Baccam is used so that the virus titer will peak at roughly 2-3 dpi, consistent with that observed in Baccam et al. [1].

Empirical, mathematical representation of the immune response
The time course for interferon (IFN) is represented as F (t) = 2 e −λg(t−tp) + e λ d (t−tp) , where λ g and λ d are the growth and decay rates of IFN respectively, and t p is the time of peak. Since F (t) ∈ [0, 1], F (t) represents the fractional amount of IFN, relative to peak IFN, at time t. In the MM, λ g = 2 d −1 , λ d = 1 d −1 and t p = 3 dpi to match experimental data. Its effect in the MM is to attenuate virus production, p, captured as where f 50 is the amount of IFN required to reduce the virus production rate to one half its normal value, i.e. p → p 2 when F = f 50 . In the MM, F (t) is normalized so as to have a maximum value of 1 (F max = 1). As such, if f 50 = 0.1, then p → p 2 when F = 0.1F max = 0.1. The time course of antibodies (Abs) is represented as where A 0 is the initial amount of Abs, and α is the growth rate of Abs. Since A(t) ∈ [0, 1], A(t) represents the fractional amount of Abs, relative to peak Abs, at time t. In the MM, α = 0.75 d −1 and A 0 = 2 × 10 −3 to match experimental data. The effect of Abs in the MM is to enhance virus clearance, c, captured as where k A (h −1 ) represents the binding affinity of Abs. The time course of cytotoxic T lymphocytes (CTLs) is represented as where λ g and λ d are the growth and decay rates of CTLs, respectively, and t p is the time of peak. Since C(t) ∈ [0, 1], C(t) represents the fractional amount of CTLs, relative to peak CTLs, at time t. In the MM, λ g = 2 d −1 , λ d = 0.4 d −1 and t p = 8 dpi to match experimental data. The effect of CTLs in the MM is to increase the rate of loss of IAV-infected cells which CTLs recognize as infected. Newly infected cells begin to present IAV peptides on their MHC-1 for CTL recognition ∼ 4 h after IAV infection [9,78]. In our MM, this corresponds approximately to half the duration of the eclipse phase (τ E /2 = 4 h): cells in E i=[1,n E /2] would not be presenting peptides, while those in E i=n E /2+1,n E would be recognizable and thus could be killed by CTLs. As such, killing of recognizably-infected cells by CTLs in the MM is captured as ∂E i (x, t) ∂t = · · · − 0 for i = 1, 2, ..., n E 2 ∂E i (x, t) ∂t = · · · − k C C(t) E i (x, t) for i = n E 2 + 1, ..., n E ∂I j (x, t) ∂t = · · · − k C C(t) I j (x, t) for j = 1, 2, ..., n I where k C (h −1 ) provides the scale for the rate of infected cell killing by CTLs.