Generalizing the Wells–Riley Infection Probability: A Superstatistical Scheme for Indoor Infection Risk Estimation

Recent evidence supports that air is the main transmission pathway of the recently identified SARS-CoV-2 coronavirus that causes COVID-19 disease. Estimating the infection risk associated with an indoor space remains an open problem due to insufficient data concerning COVID-19 outbreaks, as well as, methodological challenges arising from cases where environmental (i.e., out-of-host) and immunological (i.e., within-host) heterogeneities cannot be neglected. This work addresses these issues by introducing a generalization of the elementary Wells-Riley infection probability model. To this end, we adopted a superstatistical approach where the exposure rate parameter is gamma-distributed across subvolumes of the indoor space. This enabled us to construct a susceptible (S)–exposed (E)–infected (I) dynamics model where the Tsallis entropic index q quantifies the degree of departure from a well-mixed (i.e., homogeneous) indoor-air-environment state. A cumulative-dose mechanism is employed to describe infection activation in relation to a host’s immunological profile. We corroborate that the six-foot rule cannot guarantee the biosafety of susceptible occupants, even for exposure times as short as 15 min. Overall, our work seeks to provide a minimal (in terms of the size of the parameter space) framework for more realistic indoor SEI dynamics explorations while highlighting their Tsallisian entropic origin and the crucial yet elusive role that the innate immune system can play in shaping them. This may be useful for scientists and decision makers interested in probing different indoor biosafety protocols more thoroughly and comprehensively, thus motivating the use of nonadditive entropies in the emerging field of indoor space epidemiology.


Introduction
As a consequence of the COVID-19 pandemic, humanity has faced an unprecedented crisis affecting every sphere of society, and causing a considerable health and economic burden. Unsurprisingly, several mitigation strategies have been proposed and implemented with the ultimate goal of controlling and, if possible, preventing the transmission of SARS-CoV-2 strains. Growing evidence suggests that air is the main pathway through which SARS-CoV-2 is transmitted [1][2][3][4][5][6][7][8]. Infection by SARS-CoV-2 is, thus, much more likely to occur indoors than outdoors. In particular, an infectious occupant can spread virus-containing aerosol particles (VCAPs) into the air via exhalation, resulting in the infection of susceptible occupants according to three main scenarios [9]: (a) the short-range airborne transmission scenario [10], where the spreader and the susceptible individual are in geometrical proximity; (b) the shared-room airborne scenario, where the spreader and susceptible individual are sharing the same indoor space, thus breathing from and exhaling into the same air container; (c) the longer-distance airborne transmission [11] scenario, where the spreader and the susceptible are not in geometrical proximity (i.e., they either do not share the same room or are far apart in an ample indoor space). Realization of any of these scenarios can trigger COVID-19 outbreaks of varying epidemiological magnitude depending on complexly interwoven biological (e.g., viral infectivity and the host's immunological preparedness) and nonbiological (e.g., indoor space geometry and ventilation flow) factors. For example, there is a growing consensus that micrometer-sized VCAPs might underpin the transmission dynamics of so-called "superspreading events" [3,12,13], thus catalyzing the spread of SARS-CoV-2 in communities [8,14,15]. Although more than three years have passed since the beginning of the COVID-19 pandemic, the emergence of highly mutated and transmissible Omicron subvariants, such as XBB.1.5, poses significant threats to public health [16].
Mathematical modelling and simulation approaches to investigating airborne transmission indoors ultimately aim at developing comprehensive, quantitative methods for indoor infection-risk estimation (IIRE). The emergency of the COVID-19 pandemic accelerated multidisciplinary scientific efforts and provided a wealth of mathematical modelling approaches operating at various levels of abstraction, thus also claiming different degrees of biological plausibility. At a coarse-grained scale, SIR-like models can be insightful for understanding indoor infection dynamics at the occupant-group level (e.g., see [17][18][19][20]). On the other hand, refining the spatiotemporal scale leads to particle-based models attempting to shed light onto the transport mechanisms underlying airborne transmission and viral accumulation in the human body (e.g., see [21,22]). Typically, what all of these models have in common is an-at least-implicit connection between the indoor concentration (densities) of VCAPs and an infection rate parameter established via fluid mechanics. To this end, computational fluid dynamics (CFD) have been extensively applied to the study of airborne SARS-CoV-2 transmission pathways in different indoor environments (e.g., see [23][24][25][26]).
Due to the multiscale nature of the IIRE procedure, a crucial element of any modelling approach is the statistical frame upon which it relies to construct probability measures. Commonly, IIREs are obtained by either employing case-specific modifications of the classical Wells-Riley infection probability (WRIP) [27] or even developing entirely novel probabilistic approaches (e.g., see [28]). The WRIP's primary assumption, inherited from Poissonian statistics, is that of a homogeneous (i.e., well-mixed) indoor air environment, implying that the transmission range is predominantly short, and exposure events are probabilistically independent of each other. Stated differently, the WRIP is based on the idea that, under steady-state conditions, indoor-air-environment-property gradients (in short, steady-state-invariant gradients), which are the underlying cause of observed statistical distances between a well-mixed and a not well-mixed, i.e., heterogeneous, VCAP spatial configuration, can be neglected. The realization, however, that omitting effects of steady-state-invariant gradients may return underrated IIREs in the vicinity of an infectious source [29] has led many researchers to reconsidering the applicability of the WRIP scheme by carefully localising it (e.g., see [25,26,28,30,31]). In practice, this approach can yield satisfactory IIREs at locations of high epidemiological interest, e.g., the breathing zone of susceptible individuals without, however, providing any systematic way to integrate local WRIPs into a nonlocal (i.e., macroscopic) measure for evaluating the biosafety of the indoor air environment as a whole. An additional layer of complexity that, to the best of our knowledge, remains unexplored, mainly due to insufficient knowledge of the innate immune system dynamics [32], can be added here by considering the possibility of heterogeneous within-host responses to inhaled VCAPs.
In this work, we present a superstatistical [33] solution to the problem of integrating local WRIPs into epidemiologically relevant macroscopic measures in terms of a gamma mixture model encoding spatial fluctuations of the exposure rate parameter in an arbitrary indoor space volume. The decisive step is to conceptualize a spatial-epidemiology model where susceptible occupants may receive VCAPs via distinct (i.e., noninteracting) pathways under the influence of indoor-air-environment stochasticities. Accordingly, we consider local exposure rates to reflect the joint yet probabilistically independent action of intrinsically stochastic airborne transmissions occurring simultaneously, but via different spatial routes. The core assumption accompanying this kind of spatial thinking is that fast-and slow-dynamics timescales co-exist and are well-separated: in a steady-state indoor air environment, fast dynamics is given by the emission and inhalation of VCAPs, while slow dynamics describes how an indoor-air-environment-property-gradient field may change. This implies that macroscopic changes in the VCAP (spatial) distribution occurring during a predetermined exposure period are forbidden as long as the steady-state assumption is satisfied. Our modelling procedures suggest a frameshift from Poissonian to Paretian statistics, thus leading directly to a q-exponential WRIP, with q denoting the Tsallis entropic index [34] and admitting the interpretation of the degree of heterogeneity associated with a given steady-state indoor air environment. From a spatial-epidemiology viewpoint, q > 1 signals a nonvanishing transmission range, thus enabling construction of a distance-sensitive susceptible (S)-exposed (E)-infected (I) dynamical model where the magnitude of the rate at which occupants are transferred from the E subgroup to the I subgroup depends on immunological traits. Specifically, the Richards growth model [35] is employed to describe the dynamical relationship among viral accumulation, infection activation, and a hypothetical innate-immune-system defence mechanism orchestrated by neutralising antibodies (NAbs), and potentially enhanced by the interferon (IFN) system (for reviews covering the crucial roles that NAbs and the IFN system play in disrupting SARS-CoV-2 pathogenesis, see [36,37], respectively). This allows for investigating the interplay between out-of-host and within-host heterogeneities in a single model that, as we show, provides a simple tool for evaluating distance-based mitigation strategies such as the six-foot rule.

Preliminaries
We consider an enclosed space of rectangular volume V (m 3 ) occupied by N randomly mixed individuals split into a group of susceptible individuals of size S and infectious individuals of size F. The indoor air environment was assumed to relax into a steady state. Infectious occupants act as virus spreaders by emitting VCAPs into V. Susceptible occupants can be exposed to the virus via airborne transmission, i.e., by inhaling air samples from v br < V containing VCAPs, with v br denoting the breathing zone volume, i.e., the air volume surrounding a susceptible occupant and determining their epidemiological status. The dynamics of VCAP emission and inhalation are paced by τ rel , denoting the fastdynamics timescale. The breathing-cycle period gives the magnitude of τ rel , which is ≈3 (s). τ rel also determines the time it takes for the local equilibrium density of the indoor air particles (thus also of VCAPs) to be restored. Accordingly, perturbations in the local density of VCAPs are expected to be damped out very quickly. The time over which macroscopic environmental changes can occur, i.e., the slow-dynamics timescale, is denoted with T . The magnitude of T determines the period over which an indoor-air-environment-propertygradient field might change. It is required that T /τ rel 1. The epidemiological status of susceptible occupants may be probed at any time t ∈ [t 0 , t ], where τ := t − t 0 (h) is the occupancy time, i.e., the total time that N occupants spend in V. It is required that τ < T so that the steady-state assumption is satisfied at any t . The average maximum time that the breathing zone of a susceptible occupant is contaminated with at least one VCAP gives the exposure time; let this be denoted with τ ≤ τ. The average maximum distance over which a VCAP can be transported during τ rel delimits the transmission range; let this be denoted with ξ (m). ξ plays a similar role throughout this work to the one that the correlation length plays in superstatistical applications (e.g., see [38]). For simplicity, it is assumed that each VCAP contains an equal number of virions. An explicit connection with epidemiology is obtained via the notion of infectious quantum (IQ) (plural form abbreviation: IQa) [27] defined here as the critical number of VCAPs, Φ crit (VCAPs) ≡ 1 (IQ), that, once deposited in the body, is expected to activate an infection. Table 1 aims to ease the reader by presenting key parameters and highlighting their interdependencies. Table 1. We consider two parameter sets, namely, the out-of-host set {N, F, V, v br , τ, r, w, W, ξ, τ } and the within-host set {Φ crit , ζ, β * }. Out-of-host parameters r, w, and W denote the volumetric inhalation rate, the VCAP exhalation rate, and the ventilation rate, respectively, and are introduced in Section 2.2.1. Within-host parameters ζ and β * codetermine a susceptible occupant's immunological profile and are formally introduced in Section 2.3. λ, ρ, κ, and µ may be considered summary parameters, as they are expressed as combinations of out-of-and within-host parameters. λ and ρ represent the exposure rate parameter and VCAP density, respectively, and are introduced and reinterpreted in Section 2.2.1 and Section 2.2.2, respectively. κ denotes the number of local air environments (modelled as surrounding subvolumes of volumetric size ∝ ξ 3 ) determining a susceptible occupant's epidemiological status and is introduced in Section 2.2.2. µ gives the average time separating any pair of subsequent airborne transmission events and is introduced in Section 2.2.2.

Summary Param.
Out-of-Host Param. Within-Host Param.

Homogeneous Indoor Air Environment
If the indoor air environment is homogeneous (for a geometric description of a homogeneous indoor air environment, see Appendix A.1, a series of exposure events realised anywhere in V can be thought of as a Poisson process, where the probability that any pair of subsequent exposure events are timely separated by t = t − t 0 is given by the exponential probability density function (PDF): where λ (1/h) is a macroscopic rate designating the speed at which individuals are exposed to the virus via airborne transmission. It is common practice to consider substitution λ = r Fw WΦ crit , where w > 0 (VCAPs/h) is the rate at which VCAPs are emitted into V by an infectious occupant (i.e., VCAP exhalation rate), r > 0 (m 3 /h) is the rate at which a susceptible occupant breathes air from V (i.e., volumetric inhalation rate), and W > 0 (m 3 /h) is the rate at which clean air is supplied to V by a ventilation system [27]. The total number of VCAPs exhaled and inhaled during t by infectious and susceptible occupants is then given by Fwt and λt, respectively, and the steady-state VCAP density reads [27] from which it follows that The cumulative distribution function (CDF) associated with Equation (1) gives the classical WRIP [27]: which is assumed to serve as a good approximation for the ratio E S [27], i.e., P(t|λ) ≈ E S , where E represents the number of susceptible occupants exposed to the virus after spending t time in V (in short, exposed occupants). For t = 1 and λ = 1, (4) returns an exposure risk of ≈ 63.2%, namely, ≈ 63.2% of susceptible occupants have been exposed to the virus.

Heterogeneous Indoor Air Environment
If the indoor air environment is heterogeneous (for a geometric description of a heterogeneous indoor air environment, see Appendix A.1, getting a similar closed-form expression for the WRIP to (4) requires updating our knowledge concerning the location of VCAPs in V. The first step in this direction is to refine the spatial resolution of the IIRE procedure on the basis of our knowledge concerning the value of ξ (to gain some insight on what the order of magnitude of ξ might be, see Appendix A.1, relationship (A2)). For this, let us assume that V can be partitioned into Ω ∈ Z + nonoverlapping cubic subvolumes where Ω is chosen, so that 3 √ v ∝ ξ (for an illustration, see Figure 1). With Equation (5) at hand, we can express K as K = ∑ Ω i k i , where k i gives the number of VCAPs suspended in the i-th subvolume. By doing so, we have silently introduced a random variable, namely, the random variable k accounting for fluctuations of the local VCAP number, i.e., of the number of VCAPs suspended in a v-sized subvolume. Accordingly, k i denotes the i-th realisation of k. No assumption concerning the PDF of k is made except by requiring that the mean value of k be given by where angular brackets · indicate that the mean value calculation was performed over Ω subvolumes. k can only weakly fluctuate over τ rel ; specifically, large-amplitude fluctuations of k are only allowed over the slow-dynamics timescale T , since steady-state indoor-airenvironment conditions are assumed. Given k, let us now also introduce a local version of Equation (3): where, like k i , λ i represents a realisation of a random variable, namely, of the exposure rate random variable, λ, which accounts for local fluctuations of the exposure rate parameter. Stated differently, λ is no longer a mere phenomenological construct (as it was considered in Section 2.2.1), but it has acquired a new, microscopic interpretation: it is distributed across v-sized subvolumes with some probability f (λ) that is shaped by steady-state-invariant gradients. Hence, like k, λ can only weakly fluctuate over τ rel , since large-amplitude fluctuations of λ are forbidden during T . In light of (7), obtaining a macroscopic estimation for exposure risk statistics requires the calculation of marginal probability: which returns the mean value of p(t|λ) over f (λ) for a given t, with f (λ)p(t|λ) denoting the joint probability (i.e., the probability for a pair of subsequent exposure events to be timely separated by t given a certain value of λ). The CDF corresponding to (8) is given by This serves as a generalisation of (4), in the sense that it quantifies the probability for a susceptible occupant to become exposed to the virus after spending t hours in a heterogeneous indoor air environment. A question that naturally arises is what an epidemiologically motivated choice for f (λ) could be. In what follows, we attempt to answer this question under the macroscopic constraint that the mean value of λ, λ , equals rρ/Φ crit , i.e., λ = rρ/Φ crit .
Plausibly, an exposure event can be thought of as the outcome of κ airborne transmissions, which do not necessarily result in the inhalation of the same amount of VCAPs. Accordingly, κ is defined as a dimensionless quantity describing the number of air samples of size v inhaled during τ , i.e., where ψ (m 3 /h) re-scales r with respect to the epidemiologically relevant parameters v and τ . Within a spatial-epidemiology context, κ is interpreted as the number of v-sized subvolumes surrounding a susceptible occupant, thus determining their epidemiological status (see Figure 1). Consequently, κv gives the volumetric size of an occupant's breathing zone, i.e., v br = κv = rτ . Continuing this line of thought, the κ-th surrounding subvolume is supposed to act as an airborne transmission pathway by facilitating routes through which suspended VCAPs can reach the nearest susceptible occupant. The rate at which VCAPs are transmitted via the κ-th subvolume to the nearest susceptible occupant is determined by a random variable; let it be denoted with x (1/h). We refer to x as the transmission rate, and we note that x essentially corresponds to the κ-th spatial component of λ, i.e., Introduction of x highlights the fact that the transmission of VCAPs is an intrinsically stochastic process that can occur over a wide range of timescales averaging 1/ x , where x > 0 denotes the mean value of x. Assuming that the κ-th transmission event is probabilistically independent of all the others (which, in turn, implies that the κ-th pathway is not interacting with any of the other κ − 1 pathways), and that the transmission of a small number of VCAPs is more likely than the transmission of a large number of VCAPs, the simplest function that can be chosen for describing the PDF of x across v is the following exponential: where parameter µ = 1 x (h) denotes the average time separating any two transmission events and is defined as follows: where the ratio gives the local IQ number (i.e., the average size for an IQa dose). Like λ, x can only weakly fluctuate over τ rel , so that changes in the functional shape of h(x|µ) during T are insignificant. From (12), it follows that x can be expressed as follows: Given (10) and (13), the simplest possible functional form that one can assign to f so that λ ∼ ∑ κ i=1 x i is guaranteed is that of a gamma distribution [39] of the form: with Γ(·) being the gamma function, denoting the mean value of λ (note that λ satisfies the initially imposed macroscopic constraint), and Var(λ) = κ µ 2 denoting the variance in λ. For κ = 1 (14) becomes identical with (11) implying that v br = v (see Figure 1). When (14) is substituted into (8), we get: (16) corresponding to a Pareto Type II distribution [40] with CDF: For t = 1, µ = 1, and κ = 1, (17) returns an exposure risk of 50% (i.e., 50% of susceptible occupants were exposed to the virus). To show that (17) serves as a generalisation of (4) one simply has to calculate P(t) while having ψ be vanishingly small, i.e., lim ψ→0 which is achieved by taking either v → 0 or τ → ∞. Taking (7) and (14) together implies that k is realised by a rescaled (by a factor of c = Φ crit v/r) gamma distribution, i.e., with the required mean value Schematic illustration of the proposed spatial-epidemiology model. We consider an indoor space of volume V occupied by one infector (i.e., F = 1) and two susceptible (i.e., S 0 = 2) individuals, shown in red and blue, respectively. Round dots represent steady-state VCAP densities in V. The classical WRIP scheme assumes that the indoor air environment is homogeneous, i.e., that VCAPs are roughly uniformly spaced in V (see subfigure on the left). On the other hand, the generalized WRIP scheme does not rely on the homogeneity assumption (e.g., as we can see in the subfigure on the right, VCAP density can be higher near an infectious source). To systematically capture deviations from homogeneity, we partition V into i = 1, 2, . . . , Ω = n16 subvolumes of size v = V/Ω, where n denotes the number of subvolume layers used to fill V. For clarity, we illustrate only the front layer containing the first 16 subvolumes (subvolume boundaries are highlighted in black). Supposing that the epidemiological status of susceptible occupants is determined by a single surrounding subvolume (i.e., if κ = 1), then the subvolumes indexed with i = 9 and i = 10 correspond to the breathing zones of susceptible occupants A and B with λ 9 = x 9 and λ 10 = x 10 , respectively, representing the values of the corresponding realizations of λ with λ ∼ Gamma(1, µ) = Exp(µ) (see Equations (11) and (14)).

Infection-Activation Considerations
Once inside the body, SARS-CoV-2 can enter cells located on the surface of the upper respiratory tract via binding to the angiotensin-converting enzyme 2 (ACE2) receptor [41]. High replication levels during the first hours following exposure are correlated with the risk of developing symptomatic disease [42]. Although biological details concerning infection activation remain largely unknown, knowledge gained from long-lasting superspreading events, such as the Skagit Valley Chorale choir practice [13], suggests that a prolonged exposure time increases the risk of symptomatic disease. The critical size of the cumulative viral dose that could activate an infection leading with high certainty to symptomatic disease is impossible to measure. Nevertheless, it was recently suggested that the interplay between the size of the cumulative viral dose and the efficiency of the innate immune system plays a crucial role in determining the course of infection and disease severity [43]. Anti-SARS-CoV-2 NAbs are at the frontlines of the innate immune system, since they can inhibit the binding of the virus to the ACE2 receptor, thus offering protective immunity against SARS-CoV-2 infection [36,44,45]. In this work, we assume that infection activation is a necessary but not sufficient condition for the development of symptomatic disease.
The average number of IQa inhaled after spending t hours in V, i.e., the average size for the cumulative IQa dose up to time t, can be estimated with ( * is a reminder that this equality holds only under steady-state indoor-air-environment conditions). Without loss of generality, a relationship between Φ and a hypothetical NAborchestrated antiviral defence mechanism can be obtained via the following generalised logistic differential equation: which is known as the Richards growth model [35], where R ∈ [R 0 , R Φ→∞ = 1) is a hypothetical infection-activation biomarker, R Φ→∞ denotes the upper asymptotic bound of R, R 0 > 0 is the initial condition for R, a > 0 is the rate of change of R with respect to Φ, and ζ > 0 is a dimensionless exponent accounting for host susceptibility. Concretely, ζ is considered inversely correlated with the neutralising capacity of NAbs. Hence, the larger ζ is, the smaller the neutralising capacity of NAbs is expected to be, which, in turn, implies a lower degree of protective immunity against SARS-CoV-2 infection. The solution of (21) reads A connection with the notion of IQ is established by requiring that the inflection point of (22) is equal to 1/ζ, i.e., Φ in f l designates the Φ value for which R attains its maximum value. Qualitatively, Φ in f l can be understood as the beginning (with respect to Φ) of R's asymptotic convergence towards R Φ→∞ . From a disease biology viewpoint, Φ in f l represents a threshold value of "no return" that, once surpassed, signals a high likelihood for infection activation. For ζ = 1, we have that Φ in f l = 1 = Φ crit , implying that one IQ suffices for infection activation. Accordingly, if ζ > 1 (ζ < 1), then the host is considered to exhibit low (high) NAbattributed preparedness since Φ in f l < Φ crit (Φ in f l > Φ crit ). Φ in f l can, thus, be understood as a personalised estimation for the number of IQa required to activate an infection. This can be tuned to match an occupant's immunological profile.

Indoor Infection Dynamics
Following [27], let us now claim that (17) = E S . Then, for the initial condition S 0 = N − F > 0, the decrease in the number of susceptible occupants (or the increase in the number of exposed occupants) under steady-state indoor-air-environment conditions is given by which for a sufficiently small time step dt becomes the solution of the following differential equation [46,47]: where α is the "effective" [46,47] exposure rate, and q is the Tsallis entropic index [34]. Equation (25) can be derived by maximising the Tsallis entropy functional [46,47]: subject to constraints L 1 = t t 0 S dt and L 2 = t t 0 tS q dt, where L 1 and L 2 are Lagrange multipliers with their values being manually adjusted so that the desired values for S 0 and λ , respectively, may be obtained [46]. By using (22), we may now extend (26) to account for an I subgroup representing exposed occupants who are expected to develop symptomatic disease. We consider R to determine the value of a hypothetical infection rate β (1/h), delimiting the speed at which occupants are transferred from the E subgroup to the I subgroup, i.e., where β * > 0 (1/h) rescales R and imposes an upper asymptotic bound on β. Generally, β * can be thought of as being inversely correlated with the rate at which successful IFNsystem-driven immune responses take place (for the crucial role that the IFN system can play during the first hours following infection activation, see [42]). Hence, its value determines whether an initially activated infection is sustained or not. Altogether, our modelling considerations lead us to the following set of ordinary differential equations (ODEs) describing the indoor infection dynamics of S 0 in V: Stability analysis of the system of ODEs described in (29) is trivial; for t → ∞, all occupants were expected to have been transferred to the I subgroup, i.e., (S t→∞ → 0, E t→∞ → 0, I t→∞ → N − F) is the unique equilibrium point globally attracting from within the positively invariant region {(S, E, I) The presented SEI model belongs, from a mathematical point of view, to the class of q-SEIR models recently introduced in [19].

Scrutinising the Generalised WRIP
To better understand the implications for indoor biosafety stemming from (17), we may focus on (18). We set r = 1, and considered the case where K VCAPs were suspended in V with λ = 1 so that Fw WΦ crit = ρ Φ crit = 1. With this choice of parameters, we gain some insight into how (17) approaches (4), and how f (λ) and g(k) behave for decreasing the transmission range ξ while keeping the exposure time τ (and, thus, also v br ) fixed.
As we can see in Figure 2a, for a vanishingly small ξ (i.e., v → 0), (4) approaches its upper asymptotic bound, which is given by (17). Simply put, the classical WRIP is an overestimation of the exposure risk justified on the basis that, in the absence of any substantial knowledge concerning fluctuations of λ, the worst-case scenario may be assumed, namely, that any realisation of λ would be approximately equal to λ . The necessary condition supporting this simplification is that the VCAP density is very large, i.e., K and V should be very large and small, respectively (see also Appendix A.1). On the other hand, the generalised WRIP returns a lower but more realistic estimation for the exposure risk on the basis of the expectation that f (λ) may be well-approximated in terms of a gamma distribution. In fact, (17) serves as a more "fair" IIRE, in the sense that the chance for an occupant to become exposed to the virus after spending one hour in V is the same as tossing an unbiased coin, namely, 50% (see Figure 2a). In particular, the difference between the classical and the generalised WRIP at t = 1 decreases from ≈ 0.13 to ≈ 0.01 if one considers a tenfold increase in the number of subvolumes Ω (see Figure 2b and the corresponding legend text).  (17) for µ = 1/ψ and κ = 1/ψ with τ = τ = 2 and v ∈ (0, 2] (i.e., ψ ∈ (0, 1]). The black arrow indicates that the generalised WRIP measure approaches the classical one as v decreases. The shaded area between the two curves visualises the image of the generalised WRIP function for ψ ∈ (0, 1]. An estimation of the exposure risk difference (ERD), i.e., the difference between the generalised and the classical WRIP measures, calculated by (1 + t µ ) −κ − exp(−t) for v ∈ (0, 2] and t = 1 (=60 min) is shown in the inset graph. Colorful markers in the inset graph indicate the v key values considered in (a). In the same inset graph, the ratio S/S q→1 is plotted, where S is calculated by using the discretised version of (27) and with S q→1 denoting the Boltzmann-Gibbs entropy (see Appendix A.3, Equations (A3) and (A4), respectively). (b), for a hypothetical volume of size V = 60, we plot i = 1, 2, . . . , 30 randomly-chosen realisations of λ for v → 0, v = 0.2, and v = 2 implying that Ω → ∞, Ω = 300, and Ω = 30, respectively, since v = V Ω (see (5) The information content of an S trace is measured in terms of q entropy S (see Equation (27)) that, as we can see in the inset of Figure 2a, is inversely correlated with ξ. The loss of S-trace-related information with increasing ξ accounts for the degree of viral dispersity in V quantified in terms of the reciprocal of κ, 1/κ ∝ ξ 3 v br . Intuitively, we may understand 1/κ as a rough indicator for the likelihood that a suspended VCAP misses its target due to either the smallness of v br or the largeness of ξ (or both). Of particular interest is the case where ξ is large since, as we show in Figure A1 found in Appendix A.2, it may support a contaminated-air-sharing scenario where a pair of susceptible occupants are competing for the same IQa dose, i.e., a suspended VCAP could potentially reach any two susceptible occupants during τ rel . Thus, one might anticipate that S-trace-related information losses associated with a composite epidemiological system A ⊕ B (where A and B may represent any two susceptible occupant subgroups (or "subsystems")) should be proportional to 1/κ. Indeed, because S is nonadditive, we have that

3, Equation (A3)), and the term S[A]S[B]
κ quantifies the losses in S-trace-related information due to potential realisation of a contaminated-air-sharing scenario involving A and B as a pair. This is summarised in terms of the Tsallis entropic index q: for v → 0 (i.e., 1/κ → 0), we have that q → 1, indicating that indoor-air-environment homogeneity is restored. In turn, this implies that the information content of S is maximised, i.e., the Boltzmann-Gibbs entropic functional is recovered (see Appendix A.3, Equation (A4), and Figure 2b), and the likelihood of pairwisely sharing contaminated air is negligible.
To gain more insight into how steady-state-invariant gradients shape the statistics of λ and k, the asymptotes of f (λ) and g(k), respectively, are deduced. First, for v 1, we have that since λ ∝ v (because λ is macroscopically constrained by construction) and Var(λ) ∝ v, and k ∝ v and Var(k) ∝ v 3 apply, respectively. Inequality (a) in (31) implies that f (λ) resembles a Gaussian distribution of the form [50]: (see Figure 2c, green distribution), and, eventually, approaches the Dirac delta distribution located at λ , i.e., is the asymptote of f (λ) (see Figure 2c, blue distribution). Equation (33) sets the basis for constructing classical WRIPs (see Figure 2c). In fact, plugging f v→0 (λ) into (9) gives (4). Let us now turn our attention to Inequality (b) in (31). Its main implication is the same as previously: g(k) can be approximated with the Gaussian distribution of the form [50]: (see Figure 3, green distribution), and, eventually, approaches the Dirac delta distribution located at zero, i.e., g v→0 (k) = δ(k) (35) is the asymptote of g(k) (see Figure 3, blue distribution). Equation (35) tells us that V is partitioned into an infinite number of subvolumes, each containing an infinitesimally small number of VCAPs. This is because the mean value and variance of g(k) are proportional to v and v 3 , respectively, so that g(k) is translocated towards the origin while at the same time becoming increasingly narrow as v decreases for some finite value of K (see Figure 3).
For v → 0, it is, thus, almost certain to find an infinitesimally small amount of VCAPs anywhere in V as if VCAPs were molecules of a well-mixed gas [27].

Refining the IIRE
Of particular epidemiological interest is understanding how IIREs depend not only on indoor-air-environment properties, but also on personalised immunological traits. Towards this end, we may utilise the SEI model described by (29) as a simple tool to probe different scenarios. In Figure 4, we demonstrate how our model refines the IIRE procedure in the sense that the classical WRIP is now represented in terms of two additive components incorporating out-of-host and within-host information, namely, E and I, respectively, so that E + I → S 0 (1 − exp(− λ t)) for v → 0. Normalising SEI traces over S 0 turns them into personalised biosafety scores returning the t-dependent probabilities {s := S/S 0 , e := E/S 0 , i := I/S 0 } of escaping exposure, being exposed, and getting infected (i.e., developing symptomatic disease sometime in the near future), respectively. Decreasing the value of the host susceptibility parameter ζ translocates the I trace towards the origin, thus decelerating the transfer of occupants from the E-subgroup to the I-subgroup (see Figure 4 and the corresponding legend text) since infection activation is efficiently suppressed due to a high degree of protective immunity. This reflects the action of NAbs and is manifested as a decelerated increase in a hypothetical infection-activation biomarker, R (see (22)), induced by translocating inflection point Φ in f l away from the origin as ζ decreases (see the inset graph in Figure 4a). As one might expect, the rate of IFN-system-driven successful responses, β * , sets the speed at which occupants are transferred from the E subgroup to the I subgroup (compare Figure 4a and Figure 4b). In principle, however, the value of β * is unknown in the IIRE procedure; it may exhibit nontrivial time dependencies over viral replication dynamics and the host's immunological profile.

Evaluation of the Six-Foot Rule
We now demonstrate how our modelling procedures can aid in designing and testing distance-based mitigation strategies specifically targeting indoor spaces. We focused on how changes in the {κ, ζ, β * } triplet affect SEI dynamics under different steady-state VCAP density conditions. As a concrete example, we scrutinise the effectiveness of the six-foot rule, which dictates that the maximum number of occupants should be N max = length · width/d sa f e , V = length · width · height, where d sa f e = 1.8 (≈ 6 ft) is the minimal distance to be kept at all times between any two occupants to guarantee biosafety. To cover the worst-case scenario, we require that (a) the average maximum distance over which virions can be transmitted is determined by the minimal distance separating any pair of occupants, i.e., ξ = d sa f e , and (b) the breathing zones of susceptible occupants are uninterruptedly contaminated, i.e., we set τ = τ, so that ψ was minimised. On this basis, we proceeded with collecting i τ := i(t = t ) values for τ = 0.25, 0.5, 0.75, 1, and r = 0.54, 1.38, 3.30 (breathing rates values correspond to the mean values reported for the activities of standing, light exercise, and heavy exercise, respectively [51]), and analyse their functional dependence over ρ for different values of {ζ, β * }. First, we focus on how i τ traces behave as functions of ρ for ζ ∈ [ζ min , ζ max ] and β * = 1. We rely on the existence of an inflection point ρ in f l marking the location where the change from convex to concave in an i τ trace takes place for increasing ρ. For small values of ζ (strong protective immunity), ρ in f l is translocated away from the origin, thus decelerating the increase in i τ (see Figure 5a,c,e). On the other hand, for large values of ζ (weak protective immunity), i τ exhibits a steep increase for ρ < ρ in f l with ρ in f l being located very close to the origin (see Figure 5a,c,e). For ρ > ρ in f l , a saturation i τ plateau is gradually formed due to epidemiological spatiotemporal constraints imposed by the magnitude of ψ (see Figure 5a,c,e and the corresponding inset graphs). Specifically, there is a large value of ρ, let this be ρ * > ρ in f l , for which i τ traces converge (i.e., the range within which i τ;ρ≥ρ * values lie tends to be vanishingly small) irrespective of what the value of ζ is (see Figure 5a,c,e, and the corresponding inset graphs). This regime emerges when local VCAP densities are so high that protective immunity is inadequate in counteracting the viral threat, i.e., k i Φ in f l ; consequently, the influence of ζ on the SEI dynamics is rendered marginal. In computational practice, ρ * is obtained by finding the value of ρ for which the first-order differences of i τ drop below a threshold (for details, see the legend of Figure 5). (e), same as (a), but for r = 3.3. Inset graphs in (a,c,e) illustrate how the first-order differences trace ∆ ρ i τ = i τ;ρ+∆ρ − i τ;ρ "flattens" for increasing ρ. The following computational criterion was used to decide whether a ∆ ρ i τ -trace has "flattened": find ρ * so that ∆ ρ i τ < , = 1 × 10 −4 , is satisfied for any ρ ≥ ρ * . Once all ρ * values had been gathered for a specific r-value, we found their maximum, {ρ * }. ρ * max is provided here for convenience: it serves as a baseline value for probing the behaviour of i τ -values with respect to 1/β * . Round and square markers in (a,c,e) highlight the values of i τ at the inflection point ρ in f l for ζ = ζ min and ζ = ζ max , respectively. For clarity, only round markers are used in the inset graphs, highlighting the maxima of ∆ ρ i τ for ζ = ζ min . (b) We plot the value of i τ obtained for ρ = nρ * max , n = 1, 2, 3, . . . , 10, ρ * max = 384 versus 1/β * for r = 0.54, ζ ∈ [ζ min , ζ max ], and τ = 0.25, 0.5, 0.75, 1. (d), same as (b), but for r = 1.38 with ρ * max = 354. (f), same as (b), but for r = 3.3 with ρ * max = 245. Different colour intensities in (b,d,f) were used to visualize i τ -values in ascending n-order. Black traces appearing in (b) are particularly interesting here, as they indicate that for exposure times as short as 15 min, i τ can take values as high as ≈ 0.2 if the host's innate immune system falls short in counteracting the viral threat, i.e., if 1/β * is very small.
Lastly, we evaluate the performance of the six-foot rule by ensuring again that our calculations covered the worst-case scenario. Namely, we required that ρ ≥ ρ * and β * → ∞, implying that both NAb-orchestrated and IFN-system-driven lines of defence cannot provide sufficient protection against a viral threat. Figure 5b shows that, if N max occupants spend 15, 30, 45, and 60 min in V standing, then i τ could attain values as high as ≈0.2, ≈0.35, ≈0.5, and ≈0.6, respectively. Crucially, this finding shows that the efficiency of the six-foot rule drops significantly after 15 min, since susceptible occupants with a weak immune system might face as high a symptomatic-disease risk as ≈ 20% (notice the behaviour of the black traces shown in Figure 5b for very small 1/β * ). If N max occupants spend time in V while engaging in light exercise activities, then i τ could potentially exceed 0.4, 0.6, 0.8, and 0.9, respectively (see Figure 5d). Moreover, if N max occupants engage in heavy exercise activities, then i τ can attain values larger than 0.7 and 0.9 for τ = 0.25 and τ = 0.5, 0.75, 1, respectively, as we can see in Figure 5f. The last two cases suggest that the six-foot rule cannot guarantee the biosafety of immunologically weak susceptible occupants in indoor spaces where exercise activities are undertaken, even if the exposure time is shorter than 30 min. Obviously, for β * → 0, we have that i τ → 0, since the transfer from the E to the I subgroup is suppressed as if the IFN-system-driven response would successfully disrupt viral replication and eliminate the virus (see Figure 5b,d,f).

Discussion
Given the current rate of biosphere degradation, respiratory viruses of zoonotic origin capable of spreading via air, such as SARS-CoV-2, are expected not only to emerge more often, but also to carry an unprecedentedly high epidemic potential. Self-evidently, the probability of infection increases significantly indoors due to the very nature of the built environment, that is, enclosure. Thus, the biosafety of societies whose members spend most of their time indoors is likely compromised.
In this work, we presented a minimal (in terms of parameter space size) ODE-based model for describing indoor exposure to and potentially also infection by an airbornetransmitted virus. The initial motivation for this work was to construct a modelling framework that could stretch beyond the idea of environmental and biological homogeneity within the context of indoor air biosafety. We achieved this by generalising the WRIP on the basis of a κ-pathway spatial model that treats the breathing zone as a stochastic VCAP transmitter. Following this line of thought, we deduced that the most general PDF that could be used to statistically describe the spatial fluctuations of the exposure rate parameter and VCAP number is that of a gamma distribution. This pointed towards a Tsallisian entropic origin of transmission dynamics, with q measuring the departure from the homogeneous case. The connection between Tsallis entropy and superstatistics is usually established via χ 2 distribution (e.g., see [33]) which is nothing but a special case of the gamma distribution. Lastly, by extending exposure dynamics to account for the possibility of infection activation in relation to innate-immune-system defence mechanisms, we propose that i (defined in Section 3.2) can be used as a probability measure for estimating the risk for developing symptomatic COVID-19 disease. Although Tsallis entropy-based modelling approaches have proven useful in understanding the spread of SARS-CoV-2 among the general population (e.g., see [19,52]), the model curated in this work serves as a first attempt to zoom in on a single indoor space. We emphasise that even for a 15 min exposure under the six-foot-rule guideline, a symptomatic-disease risk as high as 20% might apply.
Let us now discuss some of the limitations of this work. First, theoretical predictions concerning the VCAP distribution remain to be verified (or dismissed) by CFD simulations. In practice, the values of parameters {κ, µ} can be estimated by analysing steady-state VCAP distributions obtained from CFD experiments by employing standard distribution parameter estimation procedures. If {r, w, W} are known, then one can investigate which combinations of v and τ offer the best description for the experimental PDF of k. Inevitably, rigorous definitions for τ and ξ then have to be provided in accordance with the simulation details. Second, very little is known about how the innate immune system reacts to SARS-CoV-2. Therefore, the presented extension of exposure dynamics based on parameter pair {ζ, β * } serves as a rather general starting point for investigating the role of immunologically heterogeneous responses to SARS-CoV-2 and might be updated in the future as our knowledge increases. In particular, future studies could explore the possibility of β * := β * (t) in order to probe the relationship among IFN-system-driven immune responses and SARS-CoV-2 replication dynamics.
Overall, our work demonstrates how superstatistics and related q-entropies can open up possibilities for systematically obtaining exposure and symptomatic-disease risk estimations in heterogeneous indoor air environments. In turn, this allows for minimising the number of parameters used when constructing personalised biosafety scores such as i. Therefore, our work provides a recipe for incorporating VCAP-related spatial information into simple ODE-based models for indoor epidemiological scenario explorations.

Conclusions
At the macroscopic level, the classical Wells-Riley infection probability results in an overrated exposure risk estimation. Locally, however, the situation is different: a classical Wells-Riley infection probability can either over-or under-estimate the associated risks depending on whether the corresponding realisation of exposure rate parameter λ i is larger or smaller, respectively, than its mean value λ .
The Tsallis entropic functional can serve as an information-theoretical starting point for exploring SEI dynamics in heterogeneous indoor air environments. To probe the geometry underlying a steady-state VCAP distribution, one may adopt a moments-expansion scheme borrowed from molecular field theory [53]. In fact, momentsexpansion schemes serve as simple, yet, informative tools for studying the spatial organisation of heterogeneous molecular systems [54,55]. Imposing periodic boundary conditions over V, we utilise the following geometrical nonuniformity index: where p i is a t-dependent vector from the coordinate system origin to the location of the i-th VCAP, d i,j is a t-dependent vector difference between the i-th and j-th VCAP locations, h i is the first-order geometric moment expanded about the i-th VCAP location calculated and normalised over 0 < K ≤ K − 1 nearest-neighboring VCAP locations, · is the Euclidean norm, and · t returns the mean value over all t instances (i.e., ensemble average). h i serves here as a t-dependent measure of the displacement tendency (or "imbalance" [53]) of the i-th VCAP about its current location, capturing the magnitude and direction of local VCAP density fluctuations. If VCAPs tend to be uniformly spaced in V (which is a plausible scenario if K and V are very large and small, respectively), then we expect that ||h i || ≈ 0 for K would exceed a lower-value threshold value. In turn, this implies that γ → 0, thus making the case for a homogeneous indoor air environment. On the other hand, for heterogeneous indoor air environments, γ is larger than zero, letting its maximum value delimited by the diagonal of V be D, i.e., γ ∈ [0, D). Given γ, nothing can really be said about ξ except that ξ ∝ γ.
Proportional relationship (A2) implies that a nonvanishing value for γ sets the ground for long-range airborne transmission since the displacement of VCAPs is possible. On the other hand, for γ → 0, the range of airborne transmission is expected to be predominantly short; large-magnitude displacements of VCAPs in homogeneous indoor air environments are unlikely to take place.
Appendix A.2. Schematic Illustration of the Contaminated-Air-Sharing Scenario Figure A1. Investigating contaminated-air-sharing scenarios. We consider the same setup as that in Figure 1, but for V partitioned in two different ways. For simplicity, we assume that v = ξ 3 . The susceptible occupants' breathing zones are highlighted and outlined in light grey and blue, respectively. On the left subfigure, we partitioned V into i = 1, 2, . . . , Ω = n4 subvolumes of size v = V/Ω, so that v br = v/4, implying that ξ > 3 √ v br . The susceptible occupants' breathing zones were embedded within the subvolume indexed with i = 3. Because ξ > 3 √ v br , the epidemiological status of A and B may be influenced by any of the surrounding subvolumes i = 1, 2, 3, 4. In fact, during τ rel , a VCAP suspended within any of i = 1, 2, 3, 4 could reach A or B. Hence, it is likely that A and B are sharing contaminated air; imagine a VCAP originating from any of i = 1, 2, 3, 4 being displaced towards and inhaled by either A or B during τ rel . On the right subfigure, we partitioned V into i = 1, 2, . . . , Ω = n64 subvolumes of size v = V/Ω, so that v br = 4v, implying that ξ < 3 √ v br .

Appendix A.3. Discretisation of the Tsallis Entropic Functional
The following discretised version of (27) was used: with s i := S(t i ), t i=1,2,... ,W = t 0 + i−1 W−1 (t − t 0 ), where p i ∝ (1 + t µ ) −κ (17) = 1 − P(t) is a time-step-specific probability of escaping exposure (i.e., not inhaling any VCAPs) after spending t hours in V. ∑ W i=1 p i = 1 with W being identified as the number of accessible microstates. Accordingly, a subgroup of susceptible occupants of size smaller than S 0 may be thought of as a subsystem visiting the i-th microstate with some probability p i . The probabilistic independence of any pair of subsystems is guaranteed by model construction, since probabilistically independent airborne transmissions are assumed throughout. For q → 1, (A3) gives us the Boltzmann-Gibbs entropy, i.e.,