original model estimates the expected number of infections as a function of variables for ventilation rate, exposure time, and number of carriers. It has also been used to study the effectiveness of methods to mitigate the risks of transmission in confined spaces, such as increasing ventilation capacity [9] .
Other researchers developed a variation of the original Wells-Riley model to account for the impact of mask-wearing. Fennelly et al. [16] used a mask-wearing-specific Wells-Riley model to conclude that mask-wearing can be useful for lowering the risk of transmission. Dai et al. [17] applied the Wells-Riley model to study the transmission risk of COVID-19 under various ventilation levels and in different indoor environments, including aircraft cabins, classrooms, offices, and buses. They concluded that the risk of airborne transmission in these indoor environments is not negligible.
Some researchers used the Wells-Riley model to explore the transmission risk in transportation. Ko et al. [18] , Chen et al. [19] , and Furuya [20] applied the model to assess the risk in commercial airline and train trips. Andrews et al. [21] conducted a complete risk analysis for various means of public transportation, including bus, rail, taxi, etc. and concluded that tuberculosis transmission may occur in South Africa's public transportation.
The various studies above, although used the Wells-Riley model for risk analysis of indoor environments, they treated the public transportation vehicle as a close-to-static indoor environment, similar to classrooms or offices, with a fixed number of people and infectors inside. The exposure time is also assumed to be constant.
However, urban public transportation vehicles travel through different areas and keep picking up and dropping off passengers in a more dynamic manner. The number of people on board varies spatially and temporally from section to section. Previous studies do not capture the dynamic and spatial characteristics of public transportation demand.
This paper aims to model the risk of transmission of airborne diseases in public transportation systems using the Wells-Riley equation and a microscopic urban rail simulation platform. The model takes into account the spatio-temporal characteristics of the demand using the system, service delivery characteristics, as well as mitigation measures related to passenger mask-wearing behavior and vehicle ventilation performance.
The paper is organized as follows. Section 2 discusses the derivation of a Wells-Riley based model incorporating the demand and supply dynamics in public transportation systems. Section 3 presents a case study that evaluates the effectiveness of different mitigation strategies for various demand levels and disease characteristics. The case study also explores the impact of system operating characteristics (such as frequency and reliability) on the risk. Section 4 concludes the paper.
The Wells-Riley model was proposed by Riley et al. [13] to evaluate the probability of infection in indoor premises, considering the intake dose of airborne pathogens in terms of the number of quanta. The equation is shown below.
where, is the probability of infection, the number of infectors (carriers), the breathing rate per person (m 3 /hour), the quanta generation rate (1/hour), the exposure time length (hour), the room ventilation rate of clean air (m 3 /hour). It is worth noticing that , the quanta generation rate (1/hour), is not a physical unit. Instead, it captures the number of infectious particles that are exhaled to the environment and also the infectivity of the particles. The value of can be back-calculated from epidemiological studies of outbreak cases using equation (1) . The Wells-Riley equation assumes a steady state of infectious particles in the environment and that the indoor air is well-mixed. The biological decay of the airborne pathogens is assumed negligible. The model has been used to evaluate risks of transmission in various indoor environments, including hospitals, office rooms, classrooms, transit vehicles, trains, and airplane cabins [13, 15, 16, 17, 18, 19, 20, 21] .
The original model (equation (1)) has been modified in the literature to incorporate additional factors. For respiratory diseases, the use of respirator or mask can reduce the number of inhaled infectious particles. The term (− !"#$ % ) in equation (1) is a measure of the intake dose. Fennelly et al. [16] developed a form of Wells-Riley equation considering mask-wearing to capture the impact of masks on the exhaled infectious dose. Similarly, Furuya [20] captures the filtration effect of respirators, including surgical mask, as inhaling protection, by adding a discounting factor towards , the breathing rate per person (m 3 /hour).
In the case of both infectors and susceptible persons wearing a mask, the probability of infection can be expressed as:
Wearing masks protects through filtering both inward (inhaling) and outward (exhaling) air movements. The infectious particles exhaled by infectors are filtered by a certain rate, + , the particle penetration rate of respirator, + = 1 represents a non-effective mask with all infectious particles from the infector exhaled into the air. + = 0 means 100% particle blockage. On the other hand, the particles inhaled by susceptible persons could be filtered by a certain rate, + ( + = 1 means completely ineffective mask and + = 0 100% particle blockage).
The various studies into the transmission risk in public transportation vehicles, neglect the passenger flow dynamics in terms of their origin and destination stations. They assume a fixed number of people inside those areas for a given time period. This ignores passenger flow characteristics, including spatial and temporal patterns in terms of number of boarding, alighting passengers, and the number of infectors and their distribution along the line based on their Origin-Destination (OD) pair. The approach proposed in the next section aims at relaxing the assumptions of earlier studies.
We illustrate the approach using the following example ( Figure 1 ) of a subway line with 4 stations. There are three passengers in the system. Passenger (pax) 1 is susceptible, while passengers 2 and 3 are infected (carriers). The train is moving from station 1 to station 4. Passenger 1 is traveling from station 2 to station 4, passenger 2 from station 1 to station 3, and passenger 3 from station 2 to station 4.
Passenger 1's risk of being infected comes from their interactions with passengers 2 and 3 when they share their transit rides in the same vehicle. The probability of not getting infected by pax 2, based on equation (1) is:
where ,-is the time for the train to travel from station 2 to 3, which is the time pax 1 and pax 2 spend together in the same vehicle. Similarly, pax 1's probability of not infected by pax 3 is:
where ,. is the time for the train to travel from station 2 to 4. Thus, the overall probability of not being infected is:
Equation (5) indicates that the infection risk of a passenger is a function of the load of infected passengers onboard each segment of their trip. Hence, it can be easily generalized for any passenger traveling on OD pair ij, assuming that the number of infected passengers on any segment of the trip is known.
Therefore, the probability of infected, 34 , is given by:
where 34 is the probability of infection for a passenger traveling on OD pair ij, 1 the number of infected carriers on segment s, 1 travel time on segment s. Given the probability 34 , the distribution of the number of infections among passengers traveling on OD ij follows the binomial distribution and the expected number of infections, 34 , is
where 34 is the OD flow from i to j (susceptible). The model can be easily extended to accommodate mask-wearing behavior. Let + be the fraction of mask-wearing passengers, 3 the fraction of infectors (carriers), and 1 the number of boarding passengers at station .
For each OD pair , passengers boarding train can be divided into 3 subgroups: 1. The probability E,F,5
VVVVVVVVVVV of a susceptible person traveling from i to j on train k and wearing a mask not infected is calculated by combining equations (6), (7), (10) and (11).
Similarly,
where E,F,5
VVVVVVVVVVVVVVVV is the probability that a non-masked susceptible passenger traveling from station to station on train is not infected. + as stated earlier is the effectiveness of a mask. Consequently, the probability of infection for a passenger traveling from i to j is given by:
3,4,5
3,4,5
Using equations (9), (14) , and (15), the expected number of infected passengers traveling on OD pair ij and train k is given by: 
The expected number of infections, 5 , onboard train k is :
The expected number of infections per trip at the OD level is:
where K is the number of train trips.
The total expected number of infections r at the system level is given by:
Another metric of interest is the probability of infection taking into account (weighted by) the distribution of the trips in the system. For a random passenger using train k the probability of infection 5 is given by:
where 5 is the total number of passengers on train k. At the system level the overall probability of infection is:
where is the total number of passengers in the system.
These metrics, and , are useful measures of risk and can be used to assess the impact of operating characteristics, as well as the effectiveness of different mitigation strategies related to mask wearing, ventilation, etc., on a consistent basis.
The above model is combined with models simulating the operations of the system to calculate the associated risk for a given OD demand, under different operating practices and mask-wearing behavior, as Figure 2 shows. OD data inferred from smart card transactions (see for example [22] ) is input to the subway line simulation.
The line simulation model generates the inputs to the risk model at the individual trip level. This includes train loads and passenger journey times between stations. Given the fraction of infectors in the general passenger population (could vary by station) and the fraction of mask-wearing passengers, the model generates the inputs required for the calculation of the infection risk at the train, OD, and system levels for different groups of passengers.
The objective of the case study is to use the methodology proposed in the previous section to evaluate the transmission risk under alternative operating strategies and various assumptions about infectiousness, ventilation, and passenger behavior in terms of mask wearing using a specific heavy rail line. The analysis also explores the effectiveness of different transmission risk mitigation strategies and the sensitivity of transmission risk to different demand levels, different infectiousness levels, and carrier rates in the population. A major subway line in one of the metropolitan areas in the U.S. provides the background for the case study.
To apply the framework described in Figure 2 , the model derived in the previous section is combined with an urban rail simulation model (SimMETRO) to assess the risk of transmission. SimMETRO is a microscopic, agent-based stochastic simulation model [23, 24] . It is detailed and designed for heavy rail system performance analysis, operations planning, signaling system evaluation, real-time control strategies evaluation and refinement, and operator-in-the-loop training. It models the network geometry, signaling system, including block design and speed code, providing speed commands for trains to follow. The input includes train schedules and dispatching and rolling stock characteristics. Passenger demand can be modeled at various levels of detail, with dynamic origin-destination (OD) flows (e.g. in 15minute intervals) at the most detailed level. A detailed dwell time model is utilized to determine a train's dwell times at stations as a function of the number of boarding and alighting passengers.
Various operating strategies can be tested and evaluated in detail by the model under different signal system configurations, passenger demand, schedules, and rolling stock characteristics. Detailed data at the train and passenger levels are generated by the simulation model, such as train runtimes, headways, dwell times, train load at each segment, individual passenger journey times, wait times, and denied boarding due to capacity constraints. The simulation model has been validated in previous studies [24] . Table 1 summarizes the base parameter values for the simulation experiments.
Base assumptions 1
Mask effectiveness in terms of particle penetration rate for exhaling (inhaling) + = 0.5 ( + = 0.5), as used in [25] , to model normal quality face coverings.
2 value = 270/hour, back-calculated using the method mentioned in [9] and [15] . Two COVID-19 outbreak cases [10] in public transportation vehicles are used to derive the value. 3 Ventilation rate of clean air in a cabin of a subway train = 1958 m 3 /hour [26] . The air exchange by door-opening at stations is neglected.
7 Risk due to waiting at stations is not considered. 8 The time period for the analysis is the PM peak direction, 4:30 to 6:30 PM 9
The regular frequency of service and schedule is used (no adjustments for changes in demand)
A number of scenarios were generated using the various factors that impact transmission. Other factors are at their base values (Table 1) . Carrier rate is the same at all stations.
The above scenarios examine the effect of each factor independently. Scenarios with the combined impact of various factors were also evaluated under different demand levels.
1. Mask-wearing behavior (0% -100%) and ventilation rates (70%-140%). 2. Mask-wearing behavior (0% -100%) and infectiousness of disease (150-270 quanta/hour).
The risk for a given scenario is evaluated for demand levels from 10% to 100% in 10% increments, to capture the demand changes during the various phases of a pandemic.
A. Base case The base case assumes a pre COVID-19 demand level (OD flows in 15 minute intervals) using the corresponding schedule (during the afternoon peak period). It is assumed that 0% of the passengers wear a mask, there is a 2% carrier rate, the infectiousness level is defined by q = 270 quanta/hour, and the ventilation system operates at its design capacity. Figure 3 shows the heatmap of the spatial distribution of risk 34 (using equation (18)), in terms of the average number of infections per train trip, for different OD pairs. Stations 6-8 correspond to the downtown area and OD groups which go from the downtown area to the end of the line are having relatively higher numbers of secondary infections. These high-risk areas could be locations for entry control, mask distribution, or stricter social distancing.
We also analyzed the contribution of each station to the overall risk. The results show that the busiest stations are not necessarily the ones that contribute the most. Total entries do not always determine the risk order. For example, in the above network, the highest ranked stations in terms of total entries are stations 8, 9, 10 and 6, but the highest ranked stations in terms of their contribution to the overall risk are 8, 6, 9, 10. Figure 4 summarizes the impact on risk at the system level (equation (19) ) of various factors at different demand levels (x-axis). The y-axis measures the risk as the ratio to the base case risk. Figure 4a shows the impact on risk of mask-wearing under various demand levels. When the demand levels are lower than 30% the risk is relatively low, and the differences among the various maskwearing cases are small. As the demand keeps increasing from 30%, the risk increases more significantly. When passengers wear masks, the risk level increase is not as sensitive compared to the non-masked cases. The sensitivity is lower as the fraction of user wearing masks increases. With 60% of passengers wearing masks, the risk level decrease from the base case is approximately 50% at 100% demand. With 100% of passengers wearing masks, the risk level decrease from the base case is approximately 73%. The mitigation effect of mask-wearing is significant. Figure 4b shows the impact of ventilation. When the demand levels are lower than 30% the risk is relatively low and the impact of the different ventilation levels is small. With demand above 30%, the risk increases more significantly. Increasing ventilation rates makes the risk level less sensitive to demand increases. Furthermore, decreasing ventilation by 20% increases risk by about 27%. On the other hand, a 20% increase in ventilation reduces risk it by only 15%. Improvements in poor ventilation can have larger impact than improvements where ventilation is satisfactory. Figure 4c shows the impact of different infectiousness levels. The behavior is similar as in Figures  4a and 4b. Figure 4d shows the impact of different carrier rates. Again, the risk is relatively low and not sensitive to the rate of infectors when demand levels are lower than 30%. At demand levels higher than 30%, the risk increases significantly.
In general, when the demand level exceeds 70%, the risk increases more rapidly, the reason for this are discussed in the next section D (operational characteristics).
We also explore the impact of combinations of different factors. The scenarios explored in the analysis are summarized in Table 2 . In the discussion that follows the risk is expressed as the ratio to the base case system risk (calculated using equation (19)). Figure 5a shows the risk heatmap for Scenario 1 at 100% and 60% demand levels, under combinations of ventilation and mask-wearing factors. Lower ventilation and mask-wearing rates result in higher transmission risk. However, differences exist between the two demand levels. At 100% demand, risk levels vary with changes in ventilation rates as well as with changes in mask-wearing popularity. With 70% ventilation rate and 0% mask-wearing rate, the risk is 60% higher than with 100% ventilation. With 100% mask-wearing rate, the risk levels are relatively low regardless of ventilation rates. At 60% demand and mask-wearing rate above 90%, the overall risk levels remain as low as 10% of the base case for all ventilation rates. On the other hand, if no one is wearing a mask, the transmission risk changes significantly as a function of the ventilation rate. The results confirm the importance of mask-wearing in mitigating the transmission risk in public transportation. Even with under-performing ventilation, the transmission risk remains relatively low if more than 90% of the passengers wear masks. Figure 5b shows the results for Scenario 2. The heatmaps illustrate the risk under combinations of mask-wearing rates and infectiousness levels. With 100% demand, risk levels vary with changes in infectiousness as well as with changes in mask-wearing rates. With infectiousness of 270 quanta per hour, the transmission risk with 0% mask-wearing is more than 2 times higher than with 100% mask-wearing. When mask-wearing is 100%, the risk levels are relatively low for all infectiousness levels. At 60% demand, and mask-wearing rate higher than 80%, the risk levels remain as low as 10% of the base case for all infectiousness values. On the other hand, if the mask-wearing rate is 0%, the transmission risk changes significantly with changes in the infectiousness levels. Figure 5c summarizes the results for Scenario 3. The risk heatmaps for 100% and 60% demand under different combinations of ventilation rates and infectiousness levels indicate that lower ventilation rates and higher infectiousness levels result, as expected, in higher transmission risk. With 100% demand and infectiousness level of 270 quanta per hour, the transmission risk is 40% higher when ventilation performs at 70% compared to when it operates at normal levels. When the ventilation rate is 140%, the risk levels are relatively low, but the risk is still sensitive to an increase in infectiousness. At 60% demand, even with ventilation rate of 140%, the overall risk level is still relatively sensitive to the increase in infectiousness. The results highlight the effectiveness of ventilation in mitigating the transmission risk in public transportation operations. However, even at high ventilation rates the risk is not negligible.
The heatmaps on Figure 5 provide some other insights. The blue line in each group is the boundary below which the risk is more than 50% lower, compared to the base case for 100% demand, and more than 90% lower for 60% demand. Comparing the areas under the blue line, it is clear that mask-wearing is more effective than ventilation (Figure 5b and 5c) . Furthermore, reduction of the risk by more than 50% of the base case is attainable with a number of combinations of mask-wearing and ventilation (or reduction of more than 90% in the case of 60% demand).
It is worth noticing that in general, in all the scenarios discussed above (for example, section B), the risk levels increase more significantly after the demand goes above 70%. This is the result of load and travel time changes as the demand increases. Figure 6a shows the boxplots for runtimes from the terminal station to the peak load station under demand level ranging from 50% to 90%. Figure 6b shows the corresponding boxplots of passenger load per train at the peak load point. The mean and variability of runtimes as well as loads increase significantly as demand increases. The system experiences congestion and delays, and trains are more crowded. The increase in the number of passengers onboard, as well as the increase in runtimes and decrease of reliability, result in the risk increasing at a higher rate for demand levels above 70%. To further explore the impact of service reliability on risk, we compare the risk under two alternative schedules, one with more regular dispatching and uniform scheduled headways during the entire peak period, the other with relatively uneven headways. Both schedules have the same average headway of about 4.3 minutes. The headway coefficient of variation (cv) of the variable schedule is 0.3 and the cv for the regular 0.2. The other settings for the experiment are the same as in Table 1 . The demand is 100% and mask-wearing is at 0%. Figure 7 illustrates the passenger load under the two schedules. The mean loads per train are similar. However, the more uniform schedule results in less variability in train loads. The 90 th percentile load for the uniform schedule is about 9% lower than the other alternative.
The variability in loads impacts the distribution of the risk. We use the risk associated with individual trains 5 (equation (17)), the expected infections of a single train trip, to measure the impact. The results show that the more variable schedule results in 5 = 8.2 infections while the regular schedule in 5 = 7.8. The standard deviation for the variable schedule's 5 is 6.0, and the regular schedule 4.8. Schedule regularity can help reduce the risk level, but more importantly it reduces the extreme cases. Figure 8 shows the histogram of the probability of infection by a random passenger associated with individual trains under the two schedules (equation (20)). With the regular schedule, the mean and standard deviation of this probability is 0.0071 and 0.0046 respectively, while with the more variable schedule, the mean and standard deviation of this probability is 0.0076 and 0.0061 respectively. More importantly the probability of infection of the worst performing train is lower for the regular schedule (0.023 compared to 0.033), a 30% decrease.
While it is important for agencies to adjust their service frequency to meet passenger demand, during a pandemic it is critical to find a balance between operating costs and passenger safety. We assume a demand level at 70% of the base case. Under normal conditions, agencies would adjust the scheduled frequency accordingly. However, under a pandemic, reducing headways may increase infection risk. We hence, explore the impact of frequency of service on risk for a given demand level. The remaining settings for the experiment are the same as Table 1 . Figure 9 shows the infection probability for a random passenger in the system under 7 different headway scenarios (equation (21)). As the headway increases from 4 to 4.5 minutes, the infection probability increases slightly. For headways exceeding 5.5 minutes, the infection probability increases sharply. Increase of headways causes increase in runtimes, dwell times, train load, and the deterioration of service reliability, so the infection probability increases accordingly. The paper explores the transmission risk of airborne viruses through the use of public transportation, in particular subway systems. The Wells-Riley model is utilized to develop an overall risk metric as a function of OD flows, actual operations, and factors such as mask-wearing, ventilation, infectiousness, and virus carrier rate. The model is integrated with a microscopic simulation model of subway operations (SimMETRO) and uses trip by trip data on passenger and train movements as inputs. The risk metric captures the spatio-temporal characteristics of demand and operations and provides means to evaluate the effectiveness of various strategies to reduce the associated risk. The paper explores, using a case study with an actual subway system, the sensitivity of overall risk to various mitigation strategies, including operational strategies, mask-wearing, ventilation rates, as well as factors such as carrier rates in the population, and infectiousness of the virus. The results show that mask-wearing and ventilation are effective mitigation methods under various demand levels, infectiousness, and carrier rates. Mask-wearing is more effective at mitigating risks. Impacts from operations and service frequency are also evaluated, emphasizing the importance of maintaining reliable, regular operations at high frequencies in order to lower transmission risks. The spatial pattern of transmission risk is also explored, identifying high-risk locations.