Accounting for the geometry of the lung in respiratory viral infections (2024)

3.1 Diffusion of extracellular virions strongly influences viral dynamics on tubular tissues

In the literature on spatial viral dynamics, an assumption is sometimes made that the extracellular viral density is approximately uniform across the tissue, or, equivalently, that free viral transport is sufficiently fast as to be effectively instantaneous [5, 30, 17]. We tested whether, for a given model cell population size and diffusion coefficient, this assumption would be equally applicable to tube-shaped tissue geometries as it is to more typical square geometries. We began by constructing model tissues in a tubular geometry with a fixed number of cells but varying aspect ratio. Throughout this experiment, we use tissues of $4096$ cells with tube circumference of $C$ cells and a length of $L$ cells such that the aspect ratio of the tissue is $C/L$. Note that aspect ratio is defined based on cell numbers in each dimension and not the actual size of the tissue. Here, and throughout this work, we define tissues in tube geometry as rectangular grids of hexagonal cells with periodic boundary conditions on the upper and lower edges of the sheet and no–flux boundary conditions on the left and right ends of the tube. For this experiment, we initiate infection with four infected cells randomly placed on the left edge of the tissue. We illustrate this setup and present simulations in Figure 2. Figure 2(a) shows the topology of the tube, and Figure 2(b) shows the unrolled representation of the tube with the experimental setup. For full details of the implementation, refer to Section 2 and the Supplementary Information.

In Figure 2(c) we plot tissue snapshots from a representative simulation, in this case with an aspect ratio of 0.25 and diffusion coefficient $D=100~{}\text{CD}^{2}\text{h}^{-1}$ (where CD is a cell diameter, here taken to be 10$\mu$m as the size of a typical epithelial cell in the respiratory tract [10]). Figure 2(c) shows typical dynamics of the model. Infection begins on the left end of the tube. It is initially sluggish, but eventually accelerates and invades deeper into the tube. These snapshots show the formation of a wide infection front moving from left to right, followed by a region of dead cells. We observe the formation of distinct, tightly clustered infection foci beyond the infection front. These distal foci are caused by extracellular viral diffusion away from the infected front, which predominantly spreads via cell–to–cell infection in accordance with our parameter selections.

In Figure 2(d), we plot the cumulative proportion of infected cells in the tissue as a time series for a range of choices of the extracellular virus diffusion coefficient $D$, ranging from zero to infinity. Tissue aspect ratio is fixed here at 0.25, and we plot ten trajectories for each value of $D$. Figure 2(d) shows cumulative infected proportion growing approximately linearly with time, but with a dramatic dependence on the value of $D$. Simulations with zero diffusion (i.e., where infection can only progress through cell–to–cell infection, shown in dark green) and those where diffusion is infinite (i.e., where viral concentration is uniform, yellow) exhibit extremely different time scales of infection. Intermediate values of diffusion give rise to relatively evenly spaced trajectories between these extremes.

As a metric of sensitivity to diffusion for a given choice of tissue aspect ratio, we compute the ratio between the time taken to (almost) completely infect the tissue with no diffusion and the time taken with infinite diffusion. We use the notation $t_{0.95}^{D}$ for the mean time taken to infect 95% of the cell sheet when the diffusion coefficient is $D$. We report the time until 95% infection instead of total infection due to the high degree of noise associated with the latter. In Figure 2(e) we plot $t_{0.95}^{0}/t_{0.95}^{\infty}$ for varying tissue aspect ratios in orange. Figure 2(e) shows that sensitivity to the diffusion coefficient is dramatically increased when the tissue aspect ratio shrinks, that is, when the tube becomes increasingly narrow and elongated. Figure 2(e) demonstrates that, compared to standard tissue geometries (where the aspect ratio is approximately 1), the value of the diffusion coefficient ($D$) imposes a far greater difference on system dynamics when the tissue geometry is tubular, representing the respiratory airways. We found that if the circumference and the length of the tube were equal (1:1), the infection time was approximately five–fold faster in the infinite diffusion case compared to the zero diffusion case, whereas when the aspect ratio was much smaller (1:256), this infinite diffusion case was around 80 times faster. We also compared the ratio of the time to 95% infection for the $D=100~{}\text{CD}^{2}\text{h}^{-1}$ case — which we use in the remainder of this work — compared to the infinite diffusion case, and show these results on the same axes in purple. While this curve shows closer agreement with the infinite diffusion infection time compared to the zero–diffusion case, there is nonetheless a more than six–fold difference in infection time in the most elongated geometry, and even when the aspect ratio is 0.25, the $D=100~{}\text{CD}^{2}\text{h}^{-1}$ case is still around 25% slower than the infinite diffusion case. Hence, even with this fairly large diffusion coefficient, we still obtain significantly different dynamics compared to when diffusion is infinite. As such, the assumption that viral diffusion can be ignored and approximated as infinite — sometimes justified on square tissue geometries [5, 30, 17] — does not apply to infections on more realistic, tubular geometries.

3.2 Lung branching structure imposes directionality of infection

To investigate the role that the branching structure of the lung imposes on viral dynamics, we constructed a simple tissue geometry which serves as an approximation of a segment of the airway tree, and perform simulations which are presented in Figure 3. An illustration of this construction is shown in Figure 3(b). Figure 3(b) shows, from left to right, a single tube of tissue which splits into two offspring tubes at each of a sequence of branching points. The regions of tissue between branching points are the generations of the branching tree. We will make the simplifying assumption throughout this work that each tube division is even and binary — that is, if the circumference of tubes in Generation $g$ is $C_{g}$, then the circumference of tubes in Generation $g+1$ is $C_{g}/2$ — and moreover, that each generation is of equal length. We briefly note that, in reality, measurements of the anatomy of the human respiratory tree suggest that airways in subsequent generations tend to only shrink by around 20%, and as such, the sum of circumferences in subsequent generations increases substantially [22, 9]. In choosing to keep this quantity fixed, our model geometry should therefore be considered a caricature of the structure of the lung, instead of an exact representation of the biological reality. This approach, however, enables us to focus our results specifically on the presence of branching structure within the model tissue, and the presence of wider and narrower regions of tissue — both of which certainly are found in the lung — without the conflating effects of a domain which grows in overall width from left to right. Moreover, this geometry allows for the airway tree to be represented as a two-dimensional grid of cells with judiciously placed periodic boundaries, ensuring computational tractability and avoiding mathematical complications near generation splits. We show a schematic of the unrolled branching tissue in Figure 3(a), which shows the pairs of periodic boundaries lining either side of each branch of tissue within a given generation, defining the structure of the branching tree. Throughout this work, we model five tissue generations of length 100 cells each, where the circumference of the Generation 1 tube is 64 cells and the Generation 5 tubes therefore have circumference 4 cells. As with the tube, we place no–flux boundary conditions on the left and right edges of the branching tree. For full details of implementation, refer to Section 2 and the Supplementary Information.

We simulated viral infection of the branching tree by initiating infection with a single infected cell, where the seed cell was placed randomly on the left edge of a given tissue generation. We carried out ten simulations for each source generation and show the resulting infected proportion time series in Figure 3(c). Figure 3(c) shows the time series of each individual simulation along with the mean trajectory (excluding simulations where the infection quickly died out). This plot shows that, with the exception of the infection seeded in Generation 1, the dynamics form a clear peak of infection, and are generally delayed the deeper into the tree the infection is seeded. In the Generation 1 case, by contrast, the infection appears to be delayed compared to the Generation 2 case, and moreover forms a plateau in infected cell proportion with no clear peak of infection.

To probe these observations further, for each of the simulations in Figure 3(c) we plotted the time series for infected proportion in each tissue generation, along with the overall time series, which we show in Figure 3(d). To track how the infection peaks seen in Figure 3(c) is affected by source generation, we indicate the mean time of the infection peak at each source generation, for which we write $t_{\text{\tiny peak}}$. As an alternative metric of the time course of infection, we also show the mean time to 50% infection, $t_{0.5}$. Figure 3(d) confirms the observation made above that the peak of infection is delayed the deeper into the branching tree the infection is seeded. Furthermore, it shows that the simulations seeded in Generation 1 have the special property that the individual time series for each tissue generation are tightly peaked and non-overlapping, meaning the overall infected proportion time series is therefore fairly flat. In each other case, the dynamics in individual tissue generations form multiple peaks and overlap substantially. However, in each of these cases, the orange curve — representing the infected proportion dynamics in Generation 1 — always forms a single peak which correlates with the overall infected proportion peak. Taken together, these observations reveal that local availability of target cells governs the rate of infection spread. Since the spread of infection is spatially limited by the diffusion of extracellular virus and cell–to–cell infection, infection is faster in regions of the branching tree with the shortest path to the largest number of susceptible cells, which is Generation 1. Here, the infection spreads at its fastest and reaches its peak. Infection spread is likewise initially sluggish when seeded in deeper tissue generations, where the narrowness of the tube restricts the number of susceptible cells within a given distance of the infection. Infections do not reach a single peak in deeper tissue generations, since the infection must travel substantial distances to spread between the branches of a given tissue generation. The upper– and lowermost branches of Generation 5, for example, are only connected via the top of the tree in Generation 1. Infections seeded in Generation 1, however, have the unique property that, as the infection front spreads down the branching tree, it reaches all the branches of each tissue generation more or less simultaneously. As such, there is little difference in the dynamics of infection in each tissue generation and the overall dynamics are relatively flat.

This finding explains the shape and position of the peak of the infected proportion time series, but does not account for why infections seeded on the left edge of Generation 1 should be slower than those seeded in Generation 2. To explore this in more depth, we ran more simulations where the single initially infected cell was placed at different depths (counting from left to right) along the branching tree in finer resolution. We report the time of the peak infected proportion along with the time to 50% infection for the ten simulations in each case along with the means, and plot our results in Figure 4. Figure 4 confirms an overall increasing trend in both the time to the infected peak and the time to 50% infection as the depth of the seed cell increases from Generation 2 onwards, however, it also shows that infections seeded further left of this reach a peak significantly later. In this case, the infection is seeded close to the top of the tree where we have imposed a no–flux boundary. The infection front, therefore, can effectively only spread in one direction (down the tree), slowing the rate of infection. We note that while estimates of the time to 50% infection are very consistent between iterations of the model regardless of the seed cell position, the time of the infection peak is subject to considerable noise when seeded close to the top of the tree. This is because in these cases as the infection front progresses down the tree, there is little variation in the target cell availability near the front, hence the infected proportion time series is flat and the peak is subject to substantial noise. These two metrics also substantially diverge for seed cell positions in the first two tissue generations, indicating that in these cases, the infected proportion time series is significantly biased. Here, the infection reaches a peak early (as the infection front spreads through Generation 1), then slows down as it spreads through the remainder of the branching tree, reaching 50% infection of the tissue later.

3.3 Tissue geometry and aspect ratio influence the fate of competing viral lineages

Recent genetic surveys of the within-host viral populations of individuals infected with influenza have shown a remarkably low degree of genetic diversity [23]. The same study also found that single nucleotide variants — the majority of which were synonymous — from early samples were very likely to be absent from samples taken later in infection. These findings suggest that genetic diversity of the within–host influenza virus population is rapidly lost, and dominated by stochastic rather than evolutionary processes. We sought to use our model to explore the dynamics of competing viral lineages within the tissue. We did so by introducing multiple, colour–coded viral populations. We used this approach — which was inspired by experimental work by Fukuyama and colleagues, who studied mice co–infected with a suite of differently-coloured influenza reporter viruses [14] — in an earlier work, as it has the additional benefit of aiding in visualising the spread of infection from multiple foci [30]. This method was implemented by seeding infection with several cells, each infected with viruses of different lineages, and determining the lineage associated with each newly-infected cell by tracking the viral and cell population associated with each lineage. All lineages were assumed to follow the same mechanics and have the same parameters. For full details of implementation, refer to Section 2.

As a starting point, we sought to compare the lineage dynamics on tubes of varying circumferences. We did so by constructing tubes of length 500 cells and different choices of circumference $C$. As in an earlier result, we ran simulations of our model by seeding infection with four initially infected cells randomly placed on the left edge of the tube, with the exception that these cells are now assumed to be infected by different viral lineages. As a control, we also compared each case with a “naive” approach, where we constructed an approximately square tissue with the same number of cells, and applied toroidal periodic boundary conditions. Note that we define a tissue as “square” by having equally many cells in width as in length, however, due to the hexagonal packing of the cells, the physical dimensions of the sheet are rectangular. In the absence of tissue edges in this case, the four seed cells were placed randomly in the toroid. This control case is representative of the standard modelling approach in the literature where the geometry of the lung is not considered. Schematics and simulations on different geometries are presented in Figure 5. We show schematics of the tube and toroid geometries in Figure 5(a)–(b). We compare the results obtained for a narrow tube (with a circumference of 8 cells) to those obtained for a wide tube (circumference of 64 cells) — along with their corresponding toroids — in Figure 5(c)–(j).

In Figure 5(c)–(d), we plot the infected proportion time series for ten simulations on the narrow and wide tube and their corresponding toroids. We also indicate the time taken to reach 50% infection of the sheet. Note that we use this metric as opposed to the time to infection peak due to the very flat and noisy trajectories of the tube time series. Figure 5(c)–(d) show that, compared to the toroid, infection of the tube is substantially slower, and the proportion of cells infected at any given time is significantly lower. Also, while infection of the larger toroid (corresponding to the $C=64$ tube) is much slower than for the smaller one (corresponding to the $C=8$ tube) — owing to the proportionally smaller initial infected population — the time scale of infection on the narrow and wide tubes is very similar (within 1%). This is despite the fact that there are eight times more cells in the wide tube case. This finding suggests that the rate of infection along the tube depends primarily on the length of the tube, and that spread along the circumference of the tube is comparatively fast (at least for circumferences up to 64 cells).

In order to visualise the spatial structure of infected cell populations of different viral lineages, we took snapshots of the final state of the tissue for different simulations and coloured each cell by the viral lineage that infected it. In Figure 5(e)–(f) we show representative snapshots of the unrolled tube tissue for the narrow and wide tubes respectively, and in Figure 5(g)–(h) we show representative snapshots of the corresponding toroids. Collectively, Figure 5(e)–(h) show that, while there is qualitatively little difference in the structure of infected populations in the toroid instances, there is a sharp difference between the two tube snapshots. In the snapshot of the narrow tube, the lineages are well-mixed near the left edge of the tube (where infection was seeded) but the orange and green populations rapidly disappear, and the purple cells also become more sparse and are also lost approximately halfway down the tube, leaving only the yellow population. In the wide tube case, however, while the green lineage is only found in the first third of the tube, all other lineages are found throughout the tube, albeit in unequal proportions. By comparison, in both toroid cases, all four viral lineages coexist in relatively similar amounts.

In order to quantify these observations, we ran 100 simulations on both the narrow and wide tube and their corresponding toroids and reported the distribution of proportions of the tissue infected by each lineage. We show these results as histograms in Figure 5(i)–(j). Figure 5(i)–(j) shows two key results. Firstly, the distributions for the toroids and tubes are qualitatively completely different. For the tubes, lineages are by far most likely to infect <10% of the tissue, with higher proportions markedly less likely, demonstrating dominance by a single lineage (or a small number of lineages). By contrast, for the toroids, density is centred at 20–30%, corresponding to a relatively even split of the tissue between the four lineages. This distribution is more tightly peaked for the larger toroid due to the larger cell population. This finding suggests that taking the typical approach to multicellular modelling, that is, in taking the cell sheet to be approximately square and with toroidal boundaries, is insufficient for capturing the lineage dynamics seen on more realistic lung tissue geometry. Figure 5(i)–(j) also show a small but not insignificant probability of a single lineage dominating the cell population (>80%) in the narrower tube. This is not observed for the wide tube.

3.4 Branching of lung tissue promotes loss of within-host viral diversity and spatial isolation of viral lineages

Having established the results in the previous section for a single tube, we next sought to explore how viral lineage dynamics were impacted by the inclusion of branching structure in the tissue geometry. To do so, we first explored how the lineage dynamics on the branching tree were impacted based on whether the infection started in the upper or lower branches. We conducted a set of 100 simulations on a branching tree where, as with the tubes, infection was initiated with four infected cells of different lineages, randomly placed on the left, open edge of the tree (left edge of Generation 1). Then, for comparison, we also carried out 100 simulations where the seed cells were randomly placed within a single branch on the right, branching edge of the tree (right edge of Generation 5). These two scenarios are chosen to represent, in an abstract way, infections of the upper and lower respiratory tract, respectively. As in Section 3.3, we again use a branching tree with maximum circumference of 64 cells and five generations of length 100 cells, such that the overall dimensions of the tissue are the same as that of the wide tube, as previously defined.Simulations on tube and branching geometries and analyses are shown in Figure 6.

We constructed a histogram of the proportion of the tissue infected by each lineage across many simulations, shown in Figure 6(a). Figure 6(a) shows a marked difference in the outcomes of viral lineages for the two different source scenarios. While the distribution for the case where the seed cells are placed on the open edge qualitatively resembles that of the wide tube in Figure 5(j), the distribution for the case where infection is seeded on the branching edge resembles a far more extreme version of that of the narrow tube (Figure 5(i)). When infection is seeded on the branching edge, lineages are fated to either infect virtually none (0–10%) or virtually all (90–100%) of the tissue.

In Figure 6(e)–(f) we show representative snapshots of the unrolled branching tissue for both cases of the source cell positions, along with the tube snapshots from the previous section for comparison (Figure 6(c)–(d)). Figure 6(b) gives a schematic of these simulations. Figure 6(e)–(f) provide clear visual evidence for the sharp contrast between the lineage dynamics of infection on the branching tree starting in the upper or lower branches. These snapshots also suggest that the tissue is more likely to be dominated by a single lineage in the branching tree case (seeded on the branching edge) compared to the narrow tube case due to the larger overall cell population. As can be seen in Figure 6(f), by the time the infection spreads the length of the tree and reaches Generation 1, there is essentially only one lineage left but the entire upper half of the tree remaining to be infected, resulting in that lineage infecting a far greater overall proportion of the total cell population.

As another descriptive tool for analysing the lineage dynamics of infections on varying tissue geometries, we next developed a notion of lineage “extinction” if a lineage is not found beyond some specified cut-off depth into the tissue (we use a cut-off of 300 cells). We denote by $P_{\text{\tiny ext}}$ the probability that, in a given geometry, a given lineage will go extinct. Figure 6(b) shows a schematic of this construction. Clearly, extinction as we have described it relies on a single direction of viral invasion, and is therefore not well defined for infection of the branching tube from the branched edge. In the narrow tube, we computed $P_{\text{\tiny ext}}=0.633$, meaning that a given lineage was more likely than not to go extinct, whereas for the wide tube, we found $P_{\text{\tiny ext}}=0.245$, meaning extinction only occurred around 25% of the time. Interestingly, on the branching tree (with source on the open edge), the probability of extinction was slightly reduced at $P_{\text{\tiny ext}}=0.200$, despite containing the same number of cells as the wide tube but with the inclusion of narrow branches within its geometry. This finding results from the fact that the branching structure acts to keep lineages separate, which slightly reduces the chance of their being crowded out. This can be seen in the lineage structure at the right end of both the wide tube and the branching tube (with source on the open edge). While in the tube case the remaining lineages are well mixed across the tube’s circumference, in the branching case, each of the narrow terminal branches is almost entirely uniform in lineage. When infection is seeded in these narrow branches, however, the branching structure does not act to separate and protect lineages since, as was observed in the narrow tube (Figure 6(c)), lineages are rapidly crowded out before they have a chance to be separated and lineage diversity is quickly lost.

To further quantify the difference in lineage structure on the branching tree (seeded on the open edge) and the wide tube, we analysed properties of the lineage dynamics along the depth of the tissue. To do so, for a given depth $d$ cells along the tissue, we examined the band of cells with depth $d-b/2$ to $d+b/2$ (where $b$ is the bandwidth, taken to be 10 cells). We sketch this computation process in Figure 7(a). For each band of cells analysed, we computed both the number of lineages found in the band, and a clustering metric, $\kappa_{\text{\tiny lin}}$, as a measure of the extent to which lineages are clustered together. This metric, which is similar to one which we developed elsewhere [30], is computed on a cell population after an infection is complete and simply reports the mean proportion of a cell’s neighbours which are of the same lineage as the cell. We show a schematic of this calculation in Figure Figure 7(b). We plot both the number of lineages and the clustering metric $\kappa_{\text{\tiny lin}}$ as a function of band depth $d$ in Figure 7(c)–(d). In both cases we indicate the mean along with the 10th and 90th percentiles. Figure 7(c)–(d) both show that for the first approximately 200 cells of depth, the two geometries are relatively indistinguishable, and the difference in the final number of lineages remaining is not substantial compared to the noise of the system. However, for $\kappa_{\text{\tiny lin}}$, we observe that, as well as there being a notable divergence in the means for the branching tree and tube case, there is also a substantial difference in the amount of variation in each case. While $\kappa_{\text{\tiny lin}}$ is an extremely noisy metric towards the right end of the tissue, it is very consistently close to 1 for the branching tube. This confirms our observation from the snapshots that while lineages are somewhat mixed in the tube case, the branches in Generation 5 of the branching geometry generally contain only a single lineage. Note that the reason for the small range of $\kappa_{\text{\tiny lin}}$ values (around [0.7,1]) is a result of the high proportion of infections arising from cell–to–cell infection, leading to a high probability of same-lineage neighbour cells. For a different ratio of cell–to–cell to cell-free infections, it is possible that the difference in $\kappa_{\text{\tiny lin}}$ trajectories for the branching tree and tube geometries would be more pronounced.

3.5 Immune response depends on the location of infection within a branching structure

Our simulations so far have not explicitly considered the role of the immune response. Resolution of infection in our model is controlled only by the availability of the susceptible cell population, a defining feature of target cell–limited models. However, actual respiratory infections in vivo do not result in the destruction of the entire lung epithelium. Instead, only reasonably localised regions are affected, such as the upper airways or the alveolar passages, before the host immune response or drug intervention results in the clearing of infection [15]. We sought to investigate the behaviour of a basic model of immune activity, and how it might interact with infections seeded in different regions of the branching tree geometry.

We introduce the following simple model for the immune response, which is presented and investigated in Figure 8. We assume some threshold proportion, $F^{\text{\tiny thresh}}$, such that, once the cumulative infected proportion of the tissue reaches $F^{\text{\tiny thresh}}$, the host immune response is triggered. At this point, there is a delay of $t_{\text{\tiny act}}$ hours while the immune response activates, then, once the response is active, we assume that the infection is immediately and totally cleared. At this point, we compute the final cumulative infected proportion $F^{\infty}$ as a measure of the cumulative damage to the tissue. We show a schematic of this construction in Figure 8(a).

As in the previous result, we sought to apply our simple model of the immune response to infections of the branching tree beginning from either the open left edge of Generation 1, or from the highly branched right edge of Generation 5. We recall that these two scenarios are chosen to represent upper– and lower–respiratory infections in vivo, respectively. In each scenario, as in Figure 3, we initiate infection with a single infected cell randomly placed on the specified edge of the tissue. We plot representative behaviour of both infection scenarios in Figure 8(b). Figure 8(b) shows the cumulative infected proportion over time for infections seeded on either the open or the branched edge. In each case we plot ten representative trajectories and the mean trajectory. The dynamics of infections progressing from either end of the tissue geometry are substantially different from each other. For infections initiated on the open edge, the cumulative infected proportion grows roughly linearly throughout, whereas for infections spreading from the branched edge, initial growth is slow and exponential, followed by a change in convexity at a cumulative infected proportion of around 0.7. There is little variation between individual simulations of the model.

We began by testing the final cumulative damage to the tissue, $F^{\infty}$, for each of these infection seed positions for varying values of the detection threshold, $F^{\text{\tiny thresh}}$. We ran 100 simulations of the model for infections of the branching tree seeded either on the open or the branched edge, and applied our immune model for a range of $F^{\text{\tiny thresh}}$ values between 0.1 and 0.5. Here, we have kept $t_{\text{\tiny act}}$ fixed at 30h. We constructed violin plots for the distribution of $F^{\infty}$ values at each $F^{\text{\tiny thresh}}$ value, and annotated these with the mean $F^{\infty}$ value. We compared this to the $F^{\infty}$ computed on the mean trajectory for both source positions, which we also show on the same axes. That is, we compare the mean of the $F^{\infty}$ on the individual trajectories and compare this to the $F^{\infty}$ on the mean trajectory. We show these results for both the open and branched edge source in Figure 8(c) and (d), respectively. Figure 8(c)–(d) show that cumulative damage increases roughly linearly on this range of detection thresholds for the open edge case, while for the branched edge case, damage increases at a sublinear rate. This corresponds to the different curvature in the cumulative infected curves for the two scenarios up to around 70% infected. Infections seeded on the branched edge are also subject to greater variation in the final damage to the tissue between iterations of the model, especially when $F^{\text{\tiny thresh}}$ takes a value around 0.2. Finally, we also note that while the mean $F^{\infty}$ and the $F^{\infty}$ of the mean trajectory agree well for the open edge source case, where the overall variation in final cumulative damage is minimal, there is a notable difference between the two for the branched edge source case. Here, especially for $F^{\text{\tiny thresh}}\in[0.2,0.4]$, the mean trajectory sustains a substantially lower final proportion of damage compared to the average damage in the individual runs. This difference in $F^{\infty}$ suggests that stochastic effects in individual runs of the model have a substantial effect on the outcome of the immune response, according to our model, and that averaged, deterministic models (as others have used in spatial viral dynamics [7, 6]) may not fully capture the observed dynamics.

Next, we decided to examine how the interplay of the two immune parameters, $F^{\text{\tiny thresh}}$ and $t_{\text{\tiny act}}$, would influence the final cumulative damage to the two infection seed positions. To do this, we computed the final infected proportion, $F^{\infty}$, on each of the 100 simulations mentioned above over a wide range of values for $F^{\text{\tiny thresh}}$ and $t_{\text{\tiny act}}$. As a control, we applied the same method to infections on tubes the same length as the branching tree and with the same circumference as either the widest branch (64 cells) or the narrowest branch (4 cells) of the tree. These tissues represent a geometry where the tube that the infection started in continues out to the length of the full branching tree, but without any of the branching structure (a tube of circumference either 64 cells on the open edge, or 4 cells on the branched edge). We show a schematic of this setup in Figure 9(a)–(b).

In Figure 9(c)–(d) we show a heatmap of $F^{\infty}$ values on a range of values for $F^{\text{\tiny thresh}}$ and $t_{\text{\tiny act}}$. We also display the $F^{\infty}=0.5$ contour in both cases, which we use as a rough threshold for whether the host survives the infection. Figure 9(c)–(d) show that, predictably, when the immune response is both sensitive and quick to activate, very little damage is sustained, whereas when the immune response requires a higher threshold of detection and more time to become active, the tissue is mostly destroyed. The two cases differ in the shape of the contours between these extremes.

In order to compare these two scenarios, we plotted the $F^{\infty}=0.5$ contours from both cases on the same axes in Figure 9(e). We also include the $F^{\infty}=0.5$ contours on the same region of $(t_{\text{\tiny act}},F^{\text{\tiny thresh}})$ parameter space for the wide and narrow tubes as defined above (that is, tubes of the same length as the branching tree and circumference equal to that of the widest or narrowest branch of the tree respectively). Figure 9(e) shows that the curves for the two source positions on the branching tree intersect, enclosing two distinct regions. This overlap is not a special property of the $F^{\infty}=0.5$ contour; other contours exhibit similar patterns (see Figure 9 and Supplementary Figure S2). In Region I of Figure 9(e), the immune response is triggered only at a high threshold, but activates quickly. In this case, infections seeded in Generation 1 (upper respiratory) are more effectively cleared by the immune response. In Region II, where the immune detection threshold is very low but there is a substantial lag in its activation, infections seeded in Generation 5 (lower respiratory) are more effectively targeted. In the curves for the two control cases — the wide and narrow tubes — we do not observe this overlap, indicating that this is a unique property of the branching tree. Our results indicate that even with an extremely simple immune response, we see heterogeneous responses to infections in the upper and lower regions of the branching tree based on the parameters of the immune response.