We present VIRIS, an agent-based modelling framework and simulator that combines airborne viral transmission, people movement and detailed architectural design. Below, we present the model, assumptions, parameter values, and its numerical implementation.
Assumptions
An infectious individual is assumed to be a moving point source of infectious aerosols. We assume that the aerosols are emitted with zero velocity, a good assumption, since we do not consider coughing or sneezing21. In reality, particle size varies, following a probability density function22,23; for simplicity, and similarly to other works4,6,7, we assume that the aerosols can be modelled using a single representative size, with a constant gravitational settling rate24, \(\gamma\). Furthermore, we are justified to neglect large droplets, since evaporation happens very quickly (seconds to minutes)25,26.
Experimental studies27,28 show that the virus deactivation rate depends on a variety of factors including relative humidity, temperature and sunlight. From data27,28, we assume a constant virus deactivation rate, \(\beta\), as reported by van Doramalen et al.29; thus is an accurate assumption for most indoor settings (temperature around 20 \(^\circ \)C and relative humidity below 60%).
The infectious aerosols are not well-mixed in an indoor space30,31; they assume a spatiotemporal distribution. Here, we focus on single storey buildings, but our methodology can be extended to multistorey buildings. The indoor space, which includes furniture, is bounded by walls or doors. Furthermore, without loss of generality, we assume perfect reflection of aerosols off walls4. Furniture that is higher than a person’s height are also modelled as walls; otherwise we assume that furniture has no effect on viral transmission, since the airborne particles can flow over them. We assume that the infectious aerosols are well-mixed along the room height; this allows to reduce the three-dimensional (3D) model to a two-dimensional (2D) model of viral transmission.
Poor ventilation is found in many diverse indoor settings (for example, schools32, care homes33, supermarkets34, courtrooms12, and restaurants6), usually when windows are closed and fresh air is mainly provided by a mechanical ventilation system ineffectively. For these poorly ventilated spaces, we can neglect the advective transport of aerosols due to air flows30,31. However, in well ventilated spaces, air flows play a more significant role in viral transmission10,14. Therefore, applying our model to well ventilated settings requires further validation with real-world data. We represent ventilation as a first-order removal of aerosols using a constant air exchange rate, \(\lambda\). To quantify ventilation quality, we use air changes per hour (ACH), which is widely used by building engineers and space managers (ACH is converted to \(\lambda\) via \(1\ \textrmACH \approx 2.78\times 10^-4\) air changes per second).
Following the Wells–Riley ansatz17,35, we, thus, model the first-order removal of infectious aerosols with rate \(\kappa\), which is a sum of the virus deactivation rate, gravitational settling rate, and ventilation rate, respectively, that is
$$\beginaligned \kappa =\beta +\gamma +\lambda . \endaligned$$
(1)
All parameter values are given in Table 1.
We approximate the diffusion coefficient of aerosols using the turbulent eddy diffusivity, an empirical function that considers local turbulent mixing36:
$$\beginaligned D = 0.8\lambda V^2/3, \endaligned$$
(2)
where V is the room volume. We assume a constant diffusion coefficient in each scenario we consider below, determined using the volume of the largest room. Movement of individuals might enhance the advective and diffusive transport of aerosols, as shown in the experiments of Mingotti et al. 37, which would then cause a faster spread of the aerosols. Also, the effects of walking speed, crowd density and human-to-source distance on the diffusion coefficient have been investigated in the experiments of Lim et al. 38; some of these effects could be incorporated in future work.
Each individual has several attributes: (I) whether they are infectious or susceptible, (II) whether they are a superspreader, (III) if they wear a mask, (IV) their walking speed, v, (V) their movement schedule. Each individual follows a schedule consisting of a list of events (e.g. “watching TV” or “attending a lecture”), and each event is characterised by a location (x, y) in a 2D plane, a start time, an end time, and an activity (resting, talking, talking loudly, exercising, and exercising intensively). Furthermore, each activity is characterised by a breathing rate and an aerosol emission rate; see Table 1.
We assume that an individual is either stationary or moving with a constant speed from an event to the next event, and that each individual always moves ahead of the end time of their current event to arrive exactly at the start time of the next event. We assume that each individual moves along the shortest path between the previous and next event locations (assuming they know the space well), and cannot move through domain boundaries and obstacles, such as walls and furniture. The trajectory \((x_i(t),y_i(t))\) of the ith individual is computed in “Numerical implementation”. Moreover, the infection risk of a susceptible individual depends on the total number of inhaled infectious aerosols (dose)17,39; we assume that each individual needs the same dose to get infected.
Modelling viral aerosols using a reaction-diffusion equation
Based on the modelling assumptions above, a reaction-diffusion equation governs the concentration of infectious aerosols, C(x, y, t), in the indoor space, which is our domain \(\Omega\):
$$\beginaligned \frac\partial C\partial t=D\nabla ^2C+S-(\kappa +R_I)C, \endaligned$$
(3)
where S is the aerosol source term due to the aerosol emissions from the infectious individuals, \(R_I\) is the removal rate of infectious aerosols due to the breathing of all individuals, and \(\nabla =(\partial /\partial x,\partial /\partial y)\), since we reduce the model to 2D. Perfect reflection of aerosols off walls corresponds to the Neumann boundary conditions \(\partial C/\partial n=0\) on \(\partial \Omega\), where n is the outward normal direction.
For m infectious individuals, the aerosol source is
$$\beginaligned S(x,y,t) = \sum _i=1^m \frac(1-\eta _i)R_ih\delta (x-x_i(t))\delta (y-y_i(t)), \endaligned$$
(4)
where \(\eta\) is the mask efficiency, R is the infectious aerosol emission rate, h is the room height, \(\delta\) is the Dirac delta function representing a point source, and the subscript i denotes the parameter value for the ith agent. All n individuals (both infectious and susceptible) due to breathing correspond to a point sink for the infectious aerosols
$$\beginaligned R_I(x,y,t) = \sum _i=1^n (1-\eta _i)\rho _i\delta (x-x_i(t))\delta (y-y_i(t)), \endaligned$$
(5)
where \(\rho\) is the breathing rate. The infectious aerosol emission rate, R, is calculated by \(R=\mu R_t\), where \(\mu\) is the viral load, which depends on the virus type and whether the infectious agent is a superspreader, and \(R_t\) is the aerosol emission rate, which depends on the activity. The viral load is estimated via
$$\beginaligned \mu =\frac\pi 6d_p^3c_v, \endaligned$$
(6)
where \(d_p\) is the average aerosol diameter, and \(c_v\) is the viral load concentration (copies/mL)40.
Spatiotemporal infection risk
Following the Wells-Riley model17,35, we assume the infection risk of a susceptible individual is an exponential function of the dose, d, the number of inhaled infectious aerosols as follows,
$$\beginaligned P(d)=1-e^-Id, \endaligned$$
(7)
where I is infectibility, which depends on the median infectious dose \(d_m\) via \(I=-\ln (0.5)/d_m\). The median infectious dose is the number of viral copies required to infect 50% of the individuals; the median dose depends on the disease41.
The total number of inhaled infectious aerosols can be calculated by
$$\beginaligned d(t)=\int _0^t(1-\eta )\rho C(x(\tau ),y(\tau ),\tau )d\tau . \endaligned$$
(8)
We model viral transmission using a discrete S-I-R (Susceptible-Infectious-Removed) model42, where ‘R’ refers to the removal of individuals from the space. Removal takes place in some settings, e.g. in a supermarket, where infected individuals enter and then leave the space. Note that ‘R’ can also refer to recovery in some models43. However, since our simulations cover less than one day, recovery is not applicable here.
Parameter values
Table 1 summarises the parameters used in the simulations, for the SARS-CoV-2 virus. Note that our model can easily simulate other airborne diseases (influenza, MERS and SARS among others) by changing the values of parameters \(\beta\), \(\mu\), \(c_v\) and \(d_m\). The value of \(c_v\) for the SARS-CoV-2 virus ranges from \(10^2\) to \(10^11\) copies/mL40; the range of \(10^8\)–\(10^9\) copies/mL is most commonly used in modelling studies12,44,45,46. Here, we take \(c_v=10^9\ \mathrm copies/mL\) (\(\mu =6.5\%\)) for a “normal” infectious individual and \(c_v=5\times 10^9\ \mathrm copies/mL\) (\(\mu =32.5\%\)) for a superspreader4; this was validated with data from a courtroom superspreader outbreak12.
For the SARS-CoV-2 virus, the median infectious dose has been found to be 52 and 256 for seroconversion and fever, respectively, using an animal model39. Following Lau et al. 4, we use \(d_m=100\). The values of \(c_v\) and \(d_m\) will be further justified using uncertainty quantification in “Courtroom”. We consider three types of masks: cotton, surgical and N95, with efficiency 50%, 60% and 90%, respectively47,48. The breathing rate and the aerosol emission rate of an individual is based on their activity, which has been provided in Table 1. When an individual is moving from one event to another, they are assumed to walk with speed 1.5 m/s49. The parameter values for walking and exercising are assumed to be the same.
Numerical implementation
The domain \(\Omega\) and the architectural design are constructed using topologicpy18,19, which has been extended to incorporate people movement for the purposes of this project. Topologicpy is a Python-based spatial modelling and topological analysis library that supports the conceptual design of architectures; details are given in a companion paper20. The first step is to create cuboids to represent rooms and corridors. The surfaces of the cuboids represent walls, floors or ceilings. Furniture is also modelled by cuboids. To simulate the movement of individuals, we generate a navigation graph, G, for the indoor space20. The vertices of this graph are the locations of all possible events, and the edges connecting the vertices represent all possible pathways an individual can move on. For each two neighbouring events, we compute on G the shortest path between the event locations20. Each individual moves with constant walking speed, and arrive exactly at the start time of the next event.
The reaction-diffusion equation (3) is solved on \(\Omega\) using a finite element method (FEM)54. We generate a mesh for \(\Omega\) using the default Frontal-Delaunay algorithm in the gmsh library55. Then we discretise the time derivative using a backward Euler scheme, discretise the space using piecewise linear elements and compute the FEM problem at each time step using the Python-based open-source FEM library scikit-fem56. Note that using scikit-fem enables smooth integration with topologicpy in a lightweight fashion, since scikit-fem is smaller than 1MB in size and only uses NumPy57 and SciPy58 libraries, which are used in topologicpy19. VIRIS runs very fast; for example, for tens of individuals with a one-day schedule, it only takes seconds to compute the concentration and the infection risk on a standard laptop. All Python scripts can be found in a GitHub repository at
link