A mesoscale agent based modeling framework for flow-mediated infection transmission in indoor occupied spaces.

The ongoing Covid-19 pandemic, and its associated public health and socioeconomic burden, has reaffirmed the necessity for a comprehensive understanding of flow-mediated infection transmission in occupied indoor spaces. This is an inherently multiscale problem, and suitable investigation approaches that can enable evidence-based decision-making for infection control strategies, interventions, and policies; will need to account for flow physics, and occupant behavior. Here, we present a mesoscale infection transmission model for human occupied indoor spaces, by integrating an agent-based human interaction model with a flow physics model for respiratory droplet dynamics and transport. We outline the mathematical and algorithmic details of the modeling framework, and demonstrate its validity using two simple simulation scenarios that verify each of the major sub-models. We then present a detailed case-study of infection transmission in a model indoor space with 60 human occupants; using a systematic set of simulations representing various flow scenarios. Data from the simulations illustrate the utility and efficacy of the devised mesoscale model in resolving flow-mediated infection transmission; and elucidate key trends in infection transmission dynamics amongst the human occupants.

Introduction

Viral epidemic outbreaks have posed major public health concerns in recent times. The currently ongoing Coronavirus Disease 2019 (or Covid-19) pandemic caused by the SARS-CoV-2 virus (and its variants) is perhaps the most severe of these outbreaks, having led to multiple waves of cases and large-scale fatalities, and a massive worldwide public health and socioeconomic burden [1]. Other recent known outbreaks of global significance include the SARS outbreak in the early 2000s, the H1N1 influenza outbreak in 2009, the MERS-CoV outbreak originating in 2012, and the Ebola virus outbreak in 2014 [2], [3], [4]. When counted with multiple dispersed influenza outbreaks, a significant number of these diseases involve respiratory system infections with documented airborne transmission modes. Mitigation of severe public health impacts of such outbreaks necessitates a comprehensive understanding of infection transmission and control measures; with a focus on airflow-mediated transmission. This is particularly critical for planning operations in human occupied indoor spaces due to the high potential of infection spread, as reported in a large number of studies [5], [6], [7], [8], [9]. Transmission of infections like SARS-CoV-2 occurs through respiratory ejecta particles, which can be mediated by transport of such particles by background air flow — leading to fomite transmission modes via direct contact; short-range transmission modes via large respiratory droplets; and long-range transmission modes via aerosolized small/micro-particles [10], [11], [12]. For indoor occupied spaces, this is further influenced by factors such as ventilation, filtration, and occupancy density [13]. Together, these aspects motivate the need to advance our understanding of flow-mediated infection transmission in human occupied spaces. With rapidly increasing high-performance computational and digital data-processing capabilities, computational models can provide a viable approach to address this need, by enabling integration of flow physics models with indoor operations data and disease pathogen biology data [14].

The underlying phenomena for flow-mediated infection transmission in such scenarios comprise an inherently multiscale problem [15], [16], which is manifested in corresponding approaches for investigating transmission and control. At the very largest scales, compartmental epidemiological models such as SIR (acronym for Susceptible–Infected–Recovered) models are employed to study population level transmission dynamics and control intervention strategies [17], [18], [19]. The subsequent coarse scale models commonly investigate infection transmission at larger building scales, using statistical infection transmission models based on well-mixed airflow assumption [20], [21], [22]. No human occupancy or interaction variables are accounted for in detail at such coarse-scale models. Additionally, the well-mixed airflow assumption is not true in most realistic operational scenarios, and often airflow patterns themselves can yield non-intuitive infection transmission patterns [9], [23]. Conversely, at the very finest scales of phenomenon being modeled, detailed respiratory flow and airflow physics including pathogen-laden respiratory droplet and ejecta formation [24], [25] are incorporated. Recently reported methods in this category also include fully resolved fine scale modeling of fluid–particle interactions, thermal effects, and detailed respiratory particle physics [26], [27], [28], [29]. However, these models are computationally expensive, requiring multiple cores of computing resources and long compute times for realistic indoor operational scenarios [29]. While they are highly effective in enabling a deeper understanding of disease transmission mechanics, they remain intractable for policy and operation decisions pertaining to infection control strategies and interventions. Computational complexity and expense, and need for specialized compute resources, are identified challenges that can prevent adoption and utility of models for infection control and policy decisions [30], [31]. This motivates the need for mesoscopic modeling approaches of low computational expense, at the intermediate scales where human occupants and their interactions as well as flow-mediated transmission phenomena are accounted for. Several existing studies have documented active agent-based models for systems comprising socially interacting human individuals. Broadly, active agent-based models have a rich history in describing interacting organism behavior. Popular versions include variants of the Vicsek model [32], [33] for active organisms, which have been extensively used to study herding and flocking behavior [34], swarm behavior and dynamics [35], collective motion in cellular colonies [36], and other forms of collective dynamics and decision making processes [37]. Similar approaches have been extended to human interaction dynamics in form of the Social Force Model [38], the modified Headed Social Force Model [39], and several variations thereof [40], which have been extensively used to study pedestrian dynamics and interacting social behavior [41]. Specifically, for infection transmission applications, there have been several attempts using active agent-based models to investigate dynamic transmission and exposure risk: for example during air-travel and airplane boarding and deplaning [42], [43]; in a crowded supermarket [44]; and in pedestrian queues [45]; amongst others. In a series of recent works [26], [27], [28], [29] agent based human dynamics models have been coupled with fully resolved fluid–particle interaction models, involving two-way coupling of flow variables and particle transport variables between fluid grid and human agents. These models have been employed to study infection transmission in a variety of scenarios including sneezing in an airplane cabin [27], transmission on a subway station [27], and transmission in a fast-food restaurant [29].

Existing agent-interaction based models used in studying infection transmission either: (a) do not incorporate flow physics and respiratory particle/droplet transport information; or (b) involves fully resolved fine-scale flow-particle and agent interactions, leading to substantial computational expense and complexity. Approaches in category (a) remain insufficient at properly resolving flow-mediated infection transmission — which is critical for respiratory infections such as the Covid-19 pandemic [10], [12]. Non-uniformity in flow-patterns in indoor airflow is widely acknowledged; yet most agent-based simulations incorporate probabilistic infection risk models which rely heavily on the assumption of well-mixed airflow [45]. Consequently, several critical factors such as ventilation, air changes per hour (ACH), airflow isolation strategies, occupancy density and occupancy patterns cannot be comprehensively accounted for in infection transmission and control investigations, thereby potentially leading to sub-optimal infection control strategies and interventions. Conversely, approaches in category (b) above, lead to simulations spanning many millions of degrees of freedom, requiring multiple hours or days of simulation times on many-core machines, for a few minutes of physical time for realistic indoor occupied spaces. These factors further substantiate the need for models at the meso-scale, where model complexity and computational expense remains reasonable for adoption for control interventions, operational decisions, and policy decisions that may involve rapid assessment of multiple ‘what if’ scenarios.

Here, we address this aforementioned modeling need by devising a mesoscopic computational modeling framework that integrates active-agent based human occupant interaction models with a customized semi-analytical flow physics based respiratory particle transport model; and, thereby, enables characterization of flow-mediated transmission amongst socially interacting occupants in an indoor space. For the agent-based approach, we use the original formulation of the Social Force Model. Here, we focus on demonstrating the model development in Sections 2, 3; outlining algorithmic details in Section 4; demonstrating validated implementation in Sections 5.1, 5.2; and illustrating the utility and efficacy of our model in characterizing flow-mediated infection transmission using a case-study in Section 6. We will not outline details of obtaining full-scale CFD models of building airflow as function of ventilation to keep with the scope of this study; and refer to our prior works [9], [23] for further details on this separate aspect. All models have been implemented in custom computer codes developed in Python using the scientific computing stack (NumPy, SciPy); and agent-data post-processing have been conducted using custom scripts based on the open-source Visualization Toolkit (VTK) library.

Mathematical description of social force models

Basic premise of the model

We will state the mathematics of the social force models, following the formulations and details in [38], for a system of active human pedestrians (each indexed by ). We assume that each pedestrian agent has a desired velocity with which the agent intends to travel within the domain of movement. The fundamental premise of agent-based social interaction models is that changes in the desired velocity are influenced by a combined set of phenomenological forces: Specifically for the classical formulation of social force model [38], this includes a collection of forces that determine collision avoidance, social affinity, and navigation to a specific target. To describe these further, we introduce a second agent velocity – the actual velocity – which is the true velocity with which the agent physically moves, and therefore changes its location as: Assuming further that each agent can only achieve a maximum achievable speed , we define a function as follows: We then relate the actual and the desired velocities as follows: Specifically, if the actual maximum possible speed of motion is less than the desired velocity, then the agent can achieve their desired velocity; otherwise the actual velocity is a scaled version of their desired velocity. This accounts for constraints on motion based on realistic human locomotion biomechanics considerations. For the social force model as proposed in [38], the combined social force and the resulting change in agent position is stated as follows: In Eq. (5), denotes a directional force to maneuver in a desired direction; denotes a collision avoidance force between agents ; denotes a wall collision avoidance force between agent and wall ; and denotes a social affinity force between agent and agent or wall or spatial location denoted by . Formulation of each force contribution is outlined in the subsequent sections. The term is a stochastic term modeling fluctuations in agent dynamics. The stochastic fluctuation term is commonly modeled as a Gaussian process with a zero mean, similar in form to the fluctuations in Langevin dynamics or active-agent Vicsek type models [33], [46].

Directional forces for each agent

Each agent is assumed to navigate towards a sequence of destinations with and being the total number of target destinations that the agent will be likely navigating. Using this, we define a desired navigation vector for an agent from its current location to the ’th destination point as follows: Assuming that the agent desires to get to with velocity , but is only able to navigate as per the actual velocity , we can state a direction force as follows: with being a relaxation time to account for pedestrian navigation velocity changes due to any acceleration or deceleration originating from social interactions. Here, physically plays the same role as particle momentum response time in particle-laden flows and discrete element methods [47], [48]; meaning that a low will enable agent to quickly approach the desired velocity; while higher will make agent motion more sluggish.

Pedestrian collision avoidance

When distance between two pedestrian agents becomes closer, territorial behavior of agents will lead to momentum changes to ensure collision avoidance. This is accounted for using a potential based repulsion force based on a repulsion potential function defined as follows: In reality, this force needs to be scaled based on agent direction , as the change in momentum will differ when a neighboring agent is located right in front of agent compared to when located behind agent . This scaled force is stated as follows: where the scaling function is defined for any given pair of vectors and as follows: and the angle represents the pedestrian agent perception angle of sight.

Wall collision avoidance

Realistic scenarios of agent dynamics must account for collision avoidance with walls and boundaries of spaces, as well as other indoor structures that substantially modify the human occupied space (for example, large furniture, space barriers etc.) We follow a potential based repulsion force formulation to account for momentum changes due to this collision avoidance as well. For any boundary wall or structure wall we identify — the point on that is the closest to pedestrian . Using this point coordinate, we can now state the wall collision avoidance force based on a wall collision potential as follows: We note that the potential force formulations in Eqs. (12) and (14) have similarities to potential-based forces in molecular dynamics, Langevin dynamics, and discrete element mechanics models. However, for maintaining consistent collision avoidance, and ensuring no spurious acceleration is allowed, the potentials must have a monotonic behavior with approach distance for both cases.

Social affinity interactions

Human agent momentum changes may occur due to attraction or affinity towards other agents, leading to the formation of socially friendly groups of individuals, or affinity towards other non-human targets in the space (for example, displays in a museum etc.). To account for this, we employ a similar potential based approach as in Eqs. (12) and (14). Specifically, we state the attraction force based on an affinity potential between an agent and another agent or object as follows: and we scale this force based on location of the agent or object with respect to the agent (perception or line of sight) using the same scaling function as described in Eq. (12) as follows: where, retains the same definition as in Section 2.3. The affinity potential is allowed to vary over time for this model, to enable changes in varying affinity towards a specific individual agent or object of interest. Similar to , and , monotonic behavior with approach distance for the functional form of is desirable to avoid spurious oscillations in agent dynamics.

Flow-mediated infection transmission model

Semi-analytical model for droplet dynamics in airflow

In this computational framework, we account for the airflow-mediated multiphysics phenomena involved in respiratory droplet dynamics from infected host to susceptible individuals through a semi-analytical approach based on a set of coupled ordinary differential equations describing the dynamics of an individual evaporating droplet. The general approach is outlined conceptually in Fig. 1. We consider the dynamics of a respiratory droplet of diameter moving in indoor flow with background air density and viscosity , and velocity (interpolated at the individual droplet location). Assuming individual droplet size is small in comparison to larger indoor flow length-scales (that is, small particle slip-velocity Reynolds numbers, as well as small Biot numbers for temperature variations), to a first-order approximation we can write the coupled set of ODEs for individual droplet as follows: where Eq. (19) accounts for drag from the surrounding flow, gravity and buoyancy effects; and Eq. (21) describes droplet diameter change based on evaporative mass transfer theory as outlined in several previous works [47], [49]. While more detailed droplet dynamics and evaporation models can be formulated as outlined in several other studies [50], [51], [52], here the simplified equations above provide a robust approach to semi-analytically integrate flow-mediated droplet multiphysics into social interaction driven agent-based modeling. The parameter is derived from evaporative mass transfer considerations, and effectively determines the droplet desiccation lifetime to a certain final droplet diameter stated as follows: It can be shown that can in turn be more precisely stated in terms of specific droplet mass transfer quantities [47] as follows: where, is the mass diffusion coefficient of droplet vapor; is the mass fraction of vapor corresponding to the droplet liquid evaluated at the droplet surface; is the mass fraction of the same vapor at free-stream; is gas density at free-stream; and is the droplet liquid density; and is the droplet mass transfer Sherwood number, commonly expressed in terms of the Ranz–Marshall correlation stated as follows: for being the droplet slip-velocity Reynolds number, and being the droplet fluid Schmidt number. Due to the continuous reduction in droplet diameter, we must ensure that physically for all scenarios; and that at the limit of we end up at the tracer-limit for the droplet, such that droplet velocity will be the same as background flow velocity. In reality for respiratory droplets, the diameter will reduce to a minimum where droplet size does not reduce further, and instead a nucleus of proteins, salts, and other biochemical by products in respiratory ejecta is formed [12], [53]. This is illustrated in Fig. 1 where individual droplet undergoes a continuous drying phase until a transition time where and thereafter the droplet is in an advection phase. Restating Eqs. (19), (21) in terms of the momentum response time , we obtain the following: such that we get: Using the coordinate system notation as outlined in Fig. 1, we can now analytically integrate the stated ODEs for droplet velocity and trajectory, for the scenario where a human agent with their face located at height from the ground releases a respiratory droplet at a velocity in the direction of its own heading. At the transition time when the droplet reaches the nucleus diameter of , we assume that the droplet coordinates are given by and velocity . The resulting expressions for velocity components and trajectory coordinates in the height-wise cross-sectional plane of the human agent, for the drying phase during can be stated as follows (detailed derivation excluded): Following the same analytical derivation, but for the non-drying advection phase at and with , we obtain the expressions for and as follows: Eqs. (28)–(37) together constitute a procedure to obtain the droplet dynamics variables as a function of air-flow velocity , other flow-related variables, and human agent variables (such as and ). For the purpose of the rest of the computational framework description, we will identify this aggregate of integrated equations as an overall respiratory droplet physics model such that: . We use this notation to devise an infection exposure and transmission model next.

Fig. 1

A visual illustration of the key aspects of the proposed agent-based flow-mediated infection transmission model; identifying the flow physics aspects for droplet transport, and (in inset) the scheme used for agent-based exposure identification based on droplet dynamics. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Infection exposure and transmission potential

The fundamental premise of infection transmission in the scenarios discussed here is that a non-infected but susceptible individual can be infected by an infected individual when virus-laden respiratory carrier particles from the latter are exposed to the former. This typically constitutes two components. First, an effective contact between infected and susceptible agents is evaluated, and an equivalent exposure to infection is identified. Second, the contact and exposure can be used to further define an infection probability. Contact and exposure will further depend upon pathogen type. Certain infections (like Ebola [54]) require actual physical contact for transmission; some pathogens require short-range contacts for transmission; and some remain viable and suspended in air for sufficiently long durations to enable both short-range and long-range transmission (such as Covid-19 [10]). Existing studies typically define contact based on a simple cutoff distance thresholding between an infected agent and a susceptible agent. Here, we use the respiratory droplet physics to inform a flow-mediated exposure model. The schematic of this approach is outlined in Fig. 1. We consider an infected agent (marked in red), and two agents around them, with one susceptible agent (marked in green). We define a local airflow , and infected agent heading direction , and project along the vector joining the infected and susceptible agent locations to obtain and (orthogonal to such that is out of the plane of the agent). The respiratory droplet release velocity is then similarly decomposed as relative to the agent (that is, infected agent velocity is added to the droplet exit velocity). These inputs are fed into to obtain droplet location and distance from infected agent ; when the model is evaluated for a threshold integration time . The threshold time can be set based on inhalation–exhalation timescales of normal human breathing; and can be typically in the range of 2–4 s. If a susceptible agent at the instance of exposure evaluation is within this droplet traversal distance, we identify an exposure event, and update contact and exposure data for the susceptible agent. When multiple infected and susceptible agents are present, this computation can be conducted pairwise for each pair of susceptible and infected agents.

Algorithmic aspects for simulation framework

Structural boundaries

The agent dynamics simulations occur across computational domains that are 2-dimensional top-down views of the occupied spaces they are navigating. The boundaries of these spaces are composed using polyline segments, denoted in 2-dimensional space using the end coordinates of the segments: , for each wall or boundary object . The same is used for identifying any indoor barrier, furnishing, or structural boundary. Target destination points can then be set as points along the line segment connecting and : for example with being a location parameter. Additionally, closest point for determining wall collision avoidance for agent can be determined based on standard geometric calculations of distance of a point from a line segment [55]. For infection exposure estimations, it is important to identify whether an infected and a susceptible agent are on the same side or opposite sides of a boundary segment. This is identified, by checking whether the line segment joining the agents intersects a boundary segment; and by marking agents on same or opposite side based on whether intersection has occurred. We note that, it is possible for flow and transmission to occur across a boundary segment, especially if the segment is not an impervious structural barrier, if the segment has a window or ventilation passage, or if overall building ventilation can move flow across the segment. Lastly, the actual boundary wall segments can be either manually specified, or obtained directly from architectural drawings, floor plans, and room or building layout data, using techniques similar to those demonstrated in our prior work [9].

Occupancy pattern initialization

We use a discrete representation of indoor space occupancy based on our previous work [23], where occupancy is outlined in terms of vector comprising the location of each individual agent or occupant: with being the number of agents; denoting the interior domain of the occupied space; and denotes the parameter vector containing data on the domain geometry including built space layout, floor-plan, and dimensions. Once initialized with an occupancy pattern , the social force model will yield subsequent outputs of . Here, using the parameters in , we use a random sequential addition process [56], [57] to generate an initial configuration of agents represented geometrically as a circular disk with a specified radius of influence. The geometry parameters define rectangular rooms, corridor spaces etc. We note that random sequential addition is suitable and robust for initializing low to moderate occupancy densities (), which was the focus here for model development. For more dense and crowded occupancy scenarios, geometric tessellation based approaches [58] or dynamic algorithms such as Lubachevsky–Stillinger algorithms [59], [60] will be necessary.

Target and destination scheduling

Specification of a sequence of targets and destinations (as stated in Section 2.2) is a critical piece of the algorithmic outline; wherein a planned sequence of point coordinate inputs are needed. This sequence: (a) can be manually assigned for simple simulation scenarios; (b) can be obtained using detailed path planning simulations using techniques similar to Dijkstra’s algorithm [44], [61]; or (c) can be obtained from operational data for the indoor space. Here, to keep with the scope of method development, we employ a more direct approach based on agent location. Specifically, doors in the space layout are input as line segments similar to boundary segments in Section 4.1 above; and at each time closest point on a door segment from a given agent is computed, and updated into the series for . We note that, these line segments can be other locations in addition to doors as well; however, for our case study presented in Section 6 we use these segments only for doors since the scenario models all occupants exiting a building.

Numerical implementation

Here, we devise a generalized numerical implementation using the social force model position and velocity updates based on an explicit Runge–Kutta (RK) method with stages; and the overall algorithmic framework has been illustrated in Algorithm 1. Explicit RK methods are chosen here for simplicity of implementation, and ability to achieve customizable high-order accuracy. For model implementation for this study, we have chosen the popular 4-stage Runge–Kutta integrator (RK4). Specifically, since the social force dynamics does not depend upon background airflow information, the agent positions and velocities are updated first (using the Runge–Kutta integrator); followed by the interpolation of airflow velocities at agent locations; and subsequently, update of respiratory droplet model to evaluate contact and exposure of a susceptible agent to an infected agent. The choice of potentials , , and requires some numerical considerations, such that the forces remain monotonic in behavior to avoid artificial acceleration and deceleration during agent movement (See also Section 2). As per recommendations in prior studies [38], [39], these potentials were chosen as exponential functions. In addition, force computations for collision avoidance and affinity interactions, as well as exposure and contact for infection, will require pairwise checks between agents. For small number of agents (as is the case here) naive pairwise checks for distance between agents can be computed directly. However as the number of agents increases, classic nearest-neighbor based techniques such as cell-lists [46], [48] will need to be employed. Our proposed framework adds the contact and exposure time in a pairwise manner for each susceptible agent with respect to each infected agent (see, for example, lines 15–23 in Algorithm 1). We note that for cases where airflow velocity is obtained from building scale CFD, or grid-based measurement data, the interpolation procedure in line 17 of Algorithm 1 will involve a grid-based interpolation operation such as binning operations or trilinear interpolations in structured Cartesian grids and cell walking algorithms based on triangle barycentric coordinates [62] in unstructured triangular grids. Finally, the presence of the random fluctuation term needs to be accounted for in the integration scheme in a suitable manner. Stochastic multi-stage integration schemes based on Ito calculus have been proposed [63] in existing studies, and the algorithm described here (Algorithm 1) can be extended to implement such methods. Alternatively, a simplified approximation can involve scaling the fluctuation contribution at each stage of Runge–Kutta integration to ensure overall fluctuating displacement scales correctly with time-step size, such that mean square displacement is correctly predicted due to fluctuations.

Verification of model implementation

Verification of social force behavior

We first establish a verification of the social force model implementation with the algorithmic details presented here. For this, we simulate the self-organization behavior of pedestrians traversing along a walkway as detailed in Helbing’s study [38]. Specifically, pedestrians enter the walkway at random positions at both ends of a rectangular walkway, and move with desired velocities being a Gaussian random variable with mean 1.34 m/s and variance 0.068 m2s-2. Agent momentum relaxation time was set to 0.5 s. Collision avoidance potential was modeled as , with m2s-2, and m. Wall collision avoidance potential was modeled as , with m2s-2 and m. The social affinity potential was set to 0. Agent perception angle of sight was set to 100 degrees; with the scalar scaling constant for the function in the force equations being set to 0.5. Above a certain pedestrian occupancy density, spontaneous formation of lanes are observed as pedestrians try to maintain a uniform walking direction. Helbing’s [38] original study showed that the average number of lanes formed () will scale with the walkway width as . Results from our implementation are computed and compared against those in [38] in Fig. 2. It is observed that the predicted number of lanes formed obeys the same scaling relation with as reported in the original case-study. Additional visualization of the case of lane formation is presented in form of animations of pedestrian agents provided in the Supplementary Material.

Fig. 2

Validation of social force implementation using the lane formation examples stated in Helbing’s original social force model study [38]. Helbing’s data, and the proposed scaling relation are illustrated in yellow dots, and solid dashed lines respectively. Current model predictions, represented in blue triangles, line up in accordance with the prior data, thereby establishing validation of current model implementation. See also lane formation animations presented in Supplementary Material. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Verification of infection exposure model

The implementation of the proposed infection exposure and transmission model is validated using a simplified model scenario, involving a group of 100 individual susceptible agents arranged in accordance with a regular Cartesian grid (10 × 10), who are facing an individual agent who is infected (see illustration in Fig. 3). A constant airflow field is imposed, and exposure is evaluated based on the modeling approach outlined in Section 4, with any individual who is exposed being marked using a binary indicator to distinguish them from remaining susceptible population (marked in black and white as shown in all illustrations in Fig. 3). Several different parameter variations were tested, to verify whether the model captures the exposure behavior in a realistic manner. First, velocity direction was kept fixed in a direction directly facing the susceptible group from the infected individual; and velocity magnitudes were varied from 0.1 m/s to 1.0 m/s. The resulting exposure patterns for increasing velocity magnitudes are presented in Fig. 3, panel a. The observations indicate an expected trend where with increasing velocity more agents are exposed, and the exposure spreads to a greater distance away from the infected agent. Second, velocity direction and magnitude were both varied to determine whether exposure trends varied consistently, and resulting trends are illustrated in Fig. 3; panel b. Third, the droplet integration threshold time was varied between the typical range of s (see Section 3.2), and resulting trends are illustrated in Fig. 3, panel c. The illustrations indicate consistent behavior of the exposure trends with both parameter variations; with the integration time showing a significant effect in determining overall exposure patterns.

Fig. 3

An illustration supporting verification of our proposed infection exposure and transmission model. Exposed individuals are marked in black, while susceptible individuals are marked in white. Panel a. depicts results for exposure with varying airflow velocity magnitude, but fixed direction Panel b. depicts the same for varying velocity magnitude and direction. Panel c. depicts trends in exposure with increasing droplet integration threshold time .

Simulation case study for a model indoor space

Having validated the social force and the infection exposure model implementation, we present a simulation case-study to demonstrate our complete computational framework (as outlined in Algorithm 1). For this case-study we consider a representative model indoor space as illustrated in Fig. 4, left panel with four rooms connected by corridors towards exits at ends of hallways on both ends. Room dimensions are specified (in meters) in the floor plan layout, and locations of doors are identified in green rectangles. We simulate a specific social dynamics scenario where the building has occupants (15 occupants in each room) who are trying to exit through the corridor out of the building through one of the exit hallways. This model representation would, for example, correspond to a scenario where individuals are exiting an office building at the end of a business day (amongst other possible scenarios). We manually impose an airflow pattern as demonstrated in Fig. 4, right panel. We note that for more realistic simulation models, this airflow pattern can be obtained based on building scale airflow simulation computations or measured airflow data (see also Section 7); which we have not considered here for conciseness of scope of the study. Instead, the imposed simplified airflow pattern enables accounting for sufficient spatial heterogeneity to demonstrate utility and efficacy of our computational approach.

Fig. 4

Illustration of the simulation set-up for the case-study outlined in Section 5.2. The indoor space layout and floorplan is shown on the left, with dimensions of each room stated. Simulation set-up comprised 15 agents in each room 1–4, with a total agents. Manually specified airflow pattern through the space is depicted on the right, with representing a fixed velocity magnitude throughout, that is varied in the simulation study as outlined in Section 5.2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For the purpose of this study, a series of four different airflow velocity magnitudes ( as per Fig. 4, right panel) are considered: 0.1 m/s; 0.2 m/s; 0.5 m/s; and 1.0 m/s. Two different scenarios are explored. In the first scenario, the walls are treated as impermeable structural barriers that partition the indoor space, and airflow does not transmit across the walls. In the second scenario, the walls are treated as structural partitions that guide and direct human occupancy, but are pervious to airflow across (that is, either through windows or open passage of air, or through air transfer across ventilation pathways). To demonstrate that our model captures the influence of flow-mediated transmission as compared to transmission without considering airflow; two additional simulation scenarios are considered. In these, the droplet flow-physics is not integrated, and contact and exposure is defined simply based on a cutoff separation distance threshold of 1.5 m; and 2.0 m as informed by many prior studies on Covid-19 transmission (especially, the commonly used 6 ft or 1.8 m social distancing recommendations) [10], [64]. The combination of these led to 12 different simulation cases; each of which was numerically integrated using a time-step of 0.25 s until the time that all agents effectively reached the exit doorways. For each case, the potentials and were chosen as described in Section 5.1, while affinity potential was set to 0. The initial droplet diameter was set to to , and set to , to mimic normal respiratory ejected particles (no coughing and sneezing); the drying coefficient was set to 2/sec; the droplet integration threshold set to 4.0 s; and the initial droplet release velocity at 40.0 m/s. The resulting agent dynamics, and exposure time for infection were analyzed to demonstrate utility and efficacy of our modeling framework.

Fig. 5 depicts a successive series of 5 snapshots (panels a–e) representing the collective dynamics of the human agents who are moving from their respective rooms towards the exit doorways. Panel f. in Fig. 5 represents the effective emergent flow pathlines for the agent collective dynamics, which illustrates the spontaneous formation of an ordered exit flow pattern (indicated by a clustering of pathlines near doors of each room). We note that the number density of agents considered here is much lower than other evacuation scenario simulations reported in literature, and hence the typically observed jamming transition of active agents is not observed here [33], [65]. In Fig. 6, we map out the exposure time for all agents when considering one single infected agent located initially in room 1 (as marked in the figure), along each pathline trajectory. Three scenarios are shown here: panel a. depicts the case for airflow 0.5 m/s with structural walls fully impervious to any flow or transmission across them; panel b. depicts the case for airflow 0.5 m/s with airflow passage and transmission across the wall boundary lines; and panel c. depicts the case for no airflow, but a constant threshold radius of 2.0 m (a little over the 6 feet/1.8 meters physical distance recommendation). The effect of airflow is observed in terms of a greater spatial extent of exposure in panel b., followed by panel a., and panel c. respectively; while panel c. shows the greatest magnitude of exposure owing to continuous contacts within the fixed radius threshold irrespective of airflow. Consideration of separation distance alone without airflow, leads to a spatial localization of estimated exposure to infection, thereby poorly resolving long-range infection transmission, as well as short-range transmission depending upon the complexity of airflow patterns. This is further illustrated in Fig. 7, where we present a statistical distribution (including the mean values) of exposure for all agents for all the 12 simulation scenarios considered. Panel a. presents the exposure data samples for all agents from scenarios with no flow or transmission across structural walls; and panel b. presents the exposure data for all agents for scenarios where flow/transmission across the boundary lines was allowed. In both cases, mean exposure generally increased with increasing magnitudes; with panel a. depicting a narrow banded sample distribution and lower mean, owing to the fact that only local flow effects are considered (no transmission across walls). Conversely, panel b. depicts a wider banded sample distribution, and higher mean compared to panel a. A greater extent of airflow transports infectious particles to more susceptible agents within the critical threshold time seconds, leading to higher extent of long-range transmission; while the constant cutoff case simply captures the localized short-range effects. The simulation data also provide some insights into the dynamics of transmission and exposure; as illustrated in the variation of exposure over time for all agents (and the mean exposure across all agents) in Fig. 8. The left column presents all cases with changing flow magnitude for no transmission across walls; while the right column presents the cases with transmission across walls (with magnitudes increasing from top to bottom panel). In all cases, exposure trends show an initial rapidly rising phase, followed by a slow varying phase reaching nearly a plateau. As expected, the obtained exposure dynamics trends is significantly mediated by airflow; as well as the nature of human agent dynamics. Additional visualization of agent dynamics, and exposure behavior is presented in form of a series of animations provided in the Supplementary Material.

Fig. 5

Illustration of interacting human agent dynamics through the occupied space for the simulation case-study. Panels a.-e. depict the agent positions at successive instances in time, as the agents navigate from their respective rooms to the exit doorways. Panel f. presents the emergent flow pathlines from the collective agent dynamics. The infected agent is marked in red square on panels a. and f. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

A comparison of exposure for flow velocity at 0.5 m/s without any transmission considered across boundary wall segments (panel a.); with transmission considered across boundary wall segments (panel b.); and for the scenario where exposure is estimated without considering airflow but with a fixed separation distance of 2.0 m (panel c.). Exposure is measured in terms of effective exposure duration, which is in the unit of time. The infected/index agent is marked in red square on all panels. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7

Comparison of agent exposure statistics for all human occupant agents across all 12 simulation scenarios outlined in Section 5.2. The labels denote the cases where airflow magnitude is varied as 0.1 m/s, 0.2 m/s, 0.5 m/s, and 1.0 m/s respectively. Labels and represent scenarios where separation distance threshold of 1.5 m and 2.0 m were employed to identify exposure. Cases without any transmission across wall boundary segments are on the left, while those with transmission across are on the right.

Fig. 8

A compilation of exposure dynamics over time, for all 8 simulation scenarios with airflow based exposure models. Left column represents all cases without any transmission across boundary wall segments ; while right column represents all cases with transmission across boundary wall segments . Airflow velocity magnitudes increase from top to bottom. Exposure is measured in terms of effective exposure duration, which is in the unit of time. Individual agent data is show in black, while agent ensemble averaged exposure dynamics curve is shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion

Here, we have outlined the computational framework for modeling flow-mediated infection transmission in an indoor space amongst socially interacting human occupants. To the best of our knowledge, this is amongst the first instances of mesoscale transmission models that integrate flow physics and respiratory droplet transport considerations with active agent-based model for human occupant behavior in a computationally inexpensive manner. While we have demonstrated the model development and application for representative case-studies; several factors and assumptions need further discussion. First, we have not discussed here the computation of an actual infection probability or risk, but instead presented our analysis in terms of contact and exposure as surrogate parameters that correlate with infection. In practice, the contact and exposure data can be further used to estimate more detailed infection probabilities based on Wells–Riley equations and dose–response theory as proposed in [21], [22], [45], [66], and based on realistic pathogen-specific data such as viral shedding, viral load, and survivability [67], [68], [69]. Additionally, all agents are assumed to be equally susceptible here. In practice, agents may have varying degrees of susceptibility. For example, agents recently infected may have acquired immunity reducing susceptibility; similarly for vaccinated and masked individuals. This would require further advancement of the current model, where exposure leads to infection probabilities only for those agents marked susceptible; similar to the premise of SIR-type compartmental epidemiological models. Second, both droplet evaporation model and droplet dynamics model can be further advanced by incorporating higher-resolution respiratory droplet physics models [51], [52] to enable more accurate representation of respiratory droplet phenomena.

A discussion on what the airflow field represents in a top-down planar view of the room is also important in this context. Indoor airflow patterns are inherently 3-dimensional, and room height can factor in based on aspects such as location of inlet and exhaust vents, doors, and windows [9], [23]. We note that averaging of the flow-field along agent height, for a representative mean height, will yield a physically representative depth-averaged flow-field that can be used as an input to the modeling framework with spatiotemporal interpolation techniques as identified in Section 4.4. Such approaches have been demonstrated in other well resolved coupled flow physics models [28], [29], where flow variables are interpolated at agent location on a planar grid. Such a description will capture the dominant large scale flow features in the room which will have the leading influence in determining the flow-mediated transmission dynamics. Avenues to integrate finer-scale flow features such as the human convective boundary layer [70], [71], and effect of local flow diverters (fans, blowers etc.) need to be further investigated for continuing advancement of the model framework. Additionally, the current model computes exposure based on droplet distance for all droplets, irrespective of where their vertical location (depth-wise in the room, as per Section 3.2) is. This enables a sufficiently conservative exposure estimate, as droplets that may have fallen below the face location of susceptible individuals are included in exposure, allowing for some re-suspension of these droplets in the depth-wise direction due to strong flow patterns in/out of the simulation plane [23]. Future advancements of this mesoscale modeling approach will comprise additional investigations on these aspects, and systematic integration with complex building scale CFD data.

Another key aspect pertains to simulation parameters and their estimation. The framework encompasses several parameters and model choices, including: potential choices for , , and , and their associated parameters and constants; relaxation time; model constants for the social force model; droplet initial velocity and size; integration time etc. Existing studies with social force simulations have demonstrated the feasibility of conducting large-scale parameter sweep studies [43]; and integrated simulation predictions with real-world human movement data and camera or closed-circuit images to do inverse parameter estimates [72], [73]. Similar approaches can be employed using this modeling framework as well, given the relatively moderate computational cost for these models even as the number of agents increase. In the context of infection transmission, model predictions can also be integrated with contact tracing data [74], [75] to identify model parameter combinations. Lastly, parameter selection process can also be set up using surrogate neural network based optimization models for a given building scenario, based on infection records for occupants of that building. In summary, despite the wide range of model parameters, several avenues exist for handling the parameter estimation problem; and will comprise another area of future model development and investigation.

Concluding remarks

While the Covid-19 pandemic has had an unprecedented impact on the world population, similar viral outbreaks have occurred before, and are likely to recur in the future to varying extents. The need for suitable quantitative models for understanding flow-mediated infection transmission; critically evaluating common assumptions such as 6-feet social distancing and well-mixed airflow; and planning control strategies and interventions remains timely and critical. We have devised and outlined a mesoscale model for flow-mediated infection transmission amongst socially interacting human occupants in an indoor space. This is amongst the first attempts to pair active agent-based models with flow physics models for respiratory infections. We have outlined various numerical and algorithmic details, and demonstrated model utility using a representative case-study. Data from the case-study clearly elucidates that accounting for airflow can lead to different infection transmission patterns compared to cases where no airflow information is considered; and that out model can capture this flow-mediated modality suitably. We emphasize the importance of analyzing flow-mediated infection transmission as an inherently multiscale problem; and in that context, our mesoscale model can complement finer scale respiratory droplet and airflow investigations, and larger scale epidemiological models. Continued investigations can enable using such models for devising scientific evidence-based guidelines for infection prevention and control during operations in facilities with human occupied indoor spaces.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.