Abstract
The global devastation of the COVID-19 pandemic has led to calls for a revolution in heating, ventilation, and air conditioning (HVAC) systems to improve indoor air quality (IAQ) [1, 2, 3], due to the dominant role of airborne transmission in disease spread [4, 5, 6]. While simple guidelines have recently been suggested to improve IAQ mainly by increasing ventilation and filtration [7, 8], this goal must be achieved in an energy-efficient and economical manner and include all air cleaning mechanisms. Here, we develop a simple protocol to directly, quantitatively, and optimally control transmission risk while minimizing energy cost. We collect a large dataset of HVAC and IAQ measurements in buildings and show how models of infectious aerosol dynamics and HVAC operation can be combined with sensor data to predict transmission risk and energy consumption. Using this data, we also verify that a simple safety guideline is able to limit transmission risk in full data-driven simulations and thus may be used to guide public health policy. Our results provide a comprehensive framework for quantitative control of transmission risk using all available air cleaning mechanisms in an indoor space while minimizing energy costs to aid in the design and automated operation of healthy, energy-efficient buildings.
Introduction
The COVID-19 pandemic disrupted the global economy and caused the worldwide shut-down of many public and private buildings essential for daily life, including schools, gyms, religious centers, and offices [9]. At first, public health guidance focused on limiting transmission from fomites and exhaled large droplets, via surface disinfection and social distancing (such as the 6 foot rule), respectively [10]. As the pandemic continued, however, it was recognized that a dominant mode of transmission of COVID-19 is through virus-laden exhaled aerosol droplets, which are small enough to remain suspended in the air for minutes to hours and become well-mixed across indoor rooms [4, 5, 6, 11, 12, 13], so the recommended mitigation strategies shifted from social distancing to masking and improved ventilation and filtration of indoor air [14, 15, 16]. Notably, a collection of leading experts in respiratory disease transmission and building science called for a “paradigm shift” in the design and operation of indoor air control systems, analogous to historical efforts to reduce pathogen transmission through food and water sources [2], as healthy indoor air is becoming recognized as a fundamental human need [1].
In this work, we propose a physics-based, data-driven strategy to achieve this paradigm shift in healthy buildings, which integrates mathematical models of airborne disease transmission with available building data streams, and apply it to data collected from multiple buildings and indoor space types. By combining airflow measurements from building management systems (BMS) with CO2 concentration data and other measures of indoor air quality (IAQ) obtained from portable sensors, the underlying physics-based models are calibrated and used to simulate the transmission risk and energy consumption for each indoor space, as operated. Using our workflow, public health officials can quantitatively assess mitigation strategies for indoor airborne disease transmission and recommend novel building control protocols to limit transmission while minimizing energy costs. Such quantitative analysis and automated building controls is necessary for preventing the next pandemic and minimizing widespread shutdowns during future pandemics.
Experiment
We collected data from diverse indoor spaces on the Massachusetts Institute of Technology (MIT) campus for several weeks in April 2022. The monitored spaces are all served by heating, ventilation, and air-conditioning (HVAC) systems, which include sensors to measure and record total supply air flow, outdoor air flow (ventilation), and supply air temperature. Temporary in-room Kaiterra sensors were also deployed to collect additional measurements relevant to IAQ, including temperature, relative humidity, CO2 concentration, total volatile organic compounds (TVOC), and size-resolved particle concentrations. Only the first three measurements are used in this study, but the entire dataset is publicly available (SI 8). In addition, we collected information about each room, including floor area, ceiling height, use case, design occupancy, and HVAC configuration (SI 2).
The overall workflow for this study is illustrated in Fig. 1. The key idea is that all data streams are integrated with physics-based models to predict and control transmission risk and energy consumption. This data fusion ensures that all space-specific variables in the model are known or accurately estimated, thus allowing different spaces to be compared with confidence.
Diagram of overall workflow. Portable indoor air quality (IAQ) sensors are placed in each monitored room, and a more durable sensor is placed on the roof to record outdoor conditions. These sensor measurements are then combined with time-series data from the heating, ventilation, and air conditioning (HVAC) system and basic properties of each room (area, ceiling height, design occupancy, etc.) for analysis.
Example time series data for four of the monitored rooms are shown in Fig. 2 and exhibit variations on the scales of days, hours, and minutes. The CO2 and humidity measurements come from the in-zone IAQ sensors, while outdoor-air flow measurements come from the BMS or are estimated from CO2 measurements. Of these measurements, CO2 concentration is the most strongly varying, as it is driven primarily by room occupancy. The peak values for Classroom 3 are significantly higher than for the other rooms, as it is does not have a forced supply of outdoor air provided by the HVAC system and is thus only naturally ventilated.
Plot of sample time series data from the study period at multiple time scales (days, hours, and minutes). The shaded region in the first (or second) column shows the time range covered in the second (or third) column. Missing data points have been filled in via interpolation.
Theory
Simple mass-balance models [17, 18] have been used for decades to study airborne transmission via aerosols and have been successfully applied to explain transmission in prior diseases [19, 20] and COVID-19 [21, 22, 23, 14, 24]. Each room is assumed to be well-mixed, such that the particle concentration can be treated as uniform throughout the room [17, 25]. This assumption is shown to produce high-quality occupancy estimates, indicating sufficient accuracy for our purposes (SI 1). The time-evolution of infectious pathogen concentration per droplet size can thus be expressed as a partial differential equation, which can be solved numerically (SI 3).
Following [14], the model can also be approximated analytically to derive a “safety guideline” that bounds the indoor reproductive number, ℛin, defined as the expected number of transmissions if an infector were present for a time τ in a given room:
where Cq represents the infectious quanta concentration in exhaled air (SI 5); Qb is the occupants’ breathing rate; pm is a mask penetration factor of aerosols; V is the volume of the room; Ns is the number of susceptible occupants; and ϵ is the desired risk tolerance. λEOA, the “equivalent outdoor air” supply rate, lumps together all infectious particle removal mechanisms including ventilation, filtration, sedimentation, and disinfection [26]. EOA quantifies each removal mechanism in terms of volumetric flow of outdoor-air ventilation that would lead to an equivalent removal rate of infectious particles, thus facilitating comparisons among disparate processes. Droplet size dependencies are integrated out by defining an effective droplet size
(SI 4). Typical ranges for these parameters and the values used in this study are provided in Supporting Information (SI 9).
Simply put, ℛin is proportional to the product of susceptible occupants and the time spent, divided by the EOA provided to the room. Thus, any holistic risk assessment and transmission mitigation strategy must consider all three dimensions. For example, simply mandating a maximum occupancy in an indoor space may not adequately reduce transmissions if the occupants spend a large amount of time in the space. Alternatively, occupancy limits may be unnecessary, if an appropriate amount of EOA is delivered to the space.
CO2 is often used as an indicator of transmission risk [27, 28, 29, 8], but we stress that it is only a partial proxy that requires care to interpret, especially when comparing across spaces with significantly different HVAC systems. To derive a CO2-based guideline, a pseudo-steady analysis of the dynamical model for CO2 concentration can be combined with Equation (1) [28], which incorporates the differences between infectious particle dynamics and CO2 dynamics stemming from sources of EOA beyond ventilation (SI 5).
Short-range respiratory flows also contribute to the risk of disease transmission in indoor spaces [30, 31, 32] and should be compared with the risk of long-range airborne transmission in any safety guideline [14]. Here, we conservatively estimate that for most of the spaces monitored, short-range effects account for less than 5% of the expected transmissions caused by long-range mixing (Fig. S25), using an experimentally validated model [33] of respiratory turbulent plumes of warm exhaled air rising by natural convection [34]. The only exception is for rooms where close-range face-to-face contact between occupants is common, where short-range risk may reach up to 30% of the long-range risk. However, the flows responsible for short-range transmission can be eliminated by requiring occupants to wear masks [34].
Results
Infection risk
The workflow can be assessed by calculating the transmission risk in the different indoor spaces using the collected data along with the full dynamical model and the pseudo-steady approximation. A key factor affecting the transmission rate is the time-varying occupancy in each space, which is then used to estimate the numbers of susceptible, Ns, and infectious, Ni, occupants. There are various occupant-counting technologies available, which utilize images, sound, or IAQ data, such as CO2 concentration, temperature, and humidity [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. We estimate occupancy by modifying previous methods that solve the full inverse problem of the CO2 concentration dynamical model and data, with a novel extension to simultaneously estimate ventilation rates (SI 1) if not measured. This method protects occupant privacy and requires no training data, so it can be easily implemented in new spaces.
The transmission rate values predicted by the full model and pseudo-steady approximation are in excellent agreement, and where they differ, the pseudo-steady model produces more conservative estimates (Fig. 3). We assess these models in two ways. In Fig. 3, we extract random segments of data from Classroom 2 (during nominally occupied hours) and plot points for each segment on a plane with axes for average occupancy and time. The color of each point corresponds to the event reproductive number calculated from the full dynamical model. Since EOA delivery is essentially constant for this space, the pseudo-steady model predicts that ℛin ∝ Occupancy × Time, which is the same trend exhibited by the color of the points in Fig. 3. This validates the use of the safety guideline, Eq. (1), to limit occupancy and time such that ℛin is below a given tolerance, as indicated by the dashed curve, which can be shifted by altering the amount of EOA provided in the room.
Transmission predictions for selected rooms in the study. Left: Scatter plot of the reproductive number for randomly chosen time periods from the full dynamical simulations, found to be in good agreement with the safety guideline [14] from the pseudo-steady formula, Eq. (1). Right: Distributions of transmission rates throughout the study period computed using full model simulation and pseudo-steady approximation, again showing good agreement. The distributions are weighted by occupancy and thus predict the expected number of transmissions if one occupant were to be infectious for one hour. Green curves are the distributions from the full dynamic model, with green lines showing the minimum, median, and maximum values. Gray curves are the distributions from the pseudo-steady model, with grey lines showing the same three statistics. Black dots show expected values assuming minimum ventilation rates and occupant density per ASHRAE standard 62.1 [49].
To assess other spaces, Fig. 3 shows distributions of transmission rates across all the monitored spaces. These points are based on a 5-minute sample rate for both models, and the distributions are weighted by the number of occupants within each point. The black dots indicate the corresponding steady-state transmission rate for a space of that type with baseline ventilation rates and occupant density, per ASHRAE standards. We see that, in almost all spaces, the worst-case transmission rate is below the baseline value as expected, since MIT buildings were deliberately operated with extra ventilation during the monitoring period to limit COVID-19 transmission. Again, the distribution of transmission rates calculated from the pseudo-steady model (gray curves) closely matches the distribution of transmission rates calculated from the full dynamic model (green curves). The median values of these distributions are generally within 1% agreement, and the maximum values differ by less that 10%.
The primary outlier from this trend is Classroom 3*, which has no mechanical ventilation, resulting in significantly smaller EOA delivery than in the other spaces. As a result, when a large number of people enter this room, the infectious particle concentration takes longer to approach the pseudo-steady values. As such, the pseudosteady model predicts a more conservative, higher transmission rate. When a large number of people leave those spaces, the opposite transient effect occurs, but since the occupancy is generally lower as people are leaving the space, those events are weighted less in the final distribution. The net result is that the pseudo-steady model predictions are slightly conservative.
The results shown for Classrooms 3 and 3* are for the same space, occupancy profile, and ventilation rate, but different levels of filtration. We assume an active filter delivering 4.5 ACH of HEPA filtration in Classroom 3 and inactive filter in Classrom 3*. When the filter is inactive, there is roughly a threefold increase in median transmission rate. We highlight this distinction because measured CO2 concentrations would be exactly the same for the two scenarios, since filtration is a form of EOA that does not impact CO2, which thus cannot be used by itself to assess transmission risk. Instead, safety guidelines [14, 28], which incorporate the differences between total EOA and ventilation (SI 5), based on a fixed ℛin tolerance, are more appropriate to account for differences in indoor spaces.
Energy and control analysis
While addressing public health concerns, there are still many opportunities to reduce energy consumption in buildings, which account for 18% of total energy use in the United States [50]. Long term building operation must balance airborne transmission risk and IAQ with energy consumption. Specifically, given a desired level of total transmission risk, buildings should operate to achieve that risk as efficiently as possible, taking advantage of all available mechanisms of infectious particle mitigation.
Formulating all removal processes in terms of EOA provides a common basis to compare various technologies in terms of cost per ACH of EOA. Under this lens, filtration, which can either be provided by in-room air cleaners or recirculated supply air, is often a much more energy-efficient source of EOA than ventilation. As an alternative, ultraviolet (UV) light can provide significant EOA by eradicating any infectious material within particles [51, 52, 53, 54] and can be installed in an upper-room configuration with shielding or as “far-UV” that is not harmful to occupants [55]. If properly installed, such systems can deliver EOA even more efficiently than filtration-based sources [56]. Thus, to optimize energy efficiency, it is necessary to consider all of these options.
Unfortunately, the primary source of EOA for many rooms is in fact ventilation. Considering the large variance in transmission risk for spaces with strongly time-varying occupancy, however, it is possible to significantly reduce energy costs by limiting extra ventilation to periods of high occupancy. Here, we analyze and compare demand-controlled ventilation (DCV) and transmission-controlled ventilation (TCV) operation modes. DCV is a feedback control mechanism implemented in most modern HVAC systems, which adjusts ventilation rates in real time to maintain a setpoint of CO2 concentration. This method thus does not consider any other EOA sources, and thus the mapping from CO2 setpoint to transmission risk can vary strongly from space to space. By contrast, we propose TCV as a novel operating mode to maintain a transmission-rate setpoint by interfacing with HVAC and accounting for other sources of EOA, all of which impact airborne disease transmission.
To quantitatively assess the inherent tradeoff between energy consumption and transmission risk, we estimate the daily energy cost for each room as operated and under various hypothetical scenarios using standard thermodynamic and equipment modeling procedures from our previous work [57, 58]. For the monitored rooms, the primary cost driver is the energy required to heat the outdoor air up to its supply temperature.
We refer to the actual operation during the monitoring period as the “Baseline” scenario. The hypothetical scenarios considered for each room are as follows:
Curtailed: ventilation is supplied at the same rate as observed in the data only during nominal occupied hours, assumed to be 8 am through 10 pm.
In-Zone Filtration: in addition to the schedule change in the “Curtailed” scenario, 2 ACH of in-room filtration is provided via standalone air cleaners active during occupied hours.
In-Zone Far UV: in addition to the schedule change in the “Curtailed” scenario, upper-room far UVC lamps are activated during occupied hours so as to provide 5 ACH of EOA.
ASHRAE Minimum: ventilation follows the “Curtailed” schedule and is further adjusted to provide the minimum amount of ventilation required in each space per ASHRAE standard 62.1 [49].
Demand Controlled (DCV): ventilation is provided by a standard demand control algorithm for a given CO2 concentration setpoint.
Transmission Controlled (TCV): ventilation is provided by a modified algorithm that provides enough ventilation to operate below a given transmission risk as calculated by the pseudo-steady model.
After estimating the time-varying ventilation that would be provided by each hypothetical strategy, transmission risk and energy cost can be calculated using the modeling approach discussed previously. These results are shown for three representative spaces in Fig. 4. Similar plots for other spaces are provided in the Supporting Information.
Summary of energy versus transmission rate tradeoffs for hypothetical ventilation scenarios. Points and shaded regions show means and joint standard deviations for daily values within the study period.
In all spaces in Fig. 4, we see there are significant opportunities to reduce energy consumption without large changes in the average and spread of the transmission rate. Simply curtailing ventilation during nighttime unoccupied hours cuts energy consumption roughly in half, with only a slight increase in transmission rate due to a small number of after-hours gatherings in that room. Adding in-room filtration can reduce transmission below baseline values while still providing significant reduction in energy cost. Based on current experimental data [55], far UV disinfection is even cheaper than in-room filtration. The ASHRAE minimum protocol further reduces energy consumption but increases the mean and spread of the transmission rate. However, by applying some of the more advanced control algorithms, transmission rate can be maintained near or below a desired threshold while maximizing energy savings. In particular, DCV at 800 ppm achieves minimum energy cost with transmission risk, while the novel TCV strategy at 0.05 per infector·h forgoes some of the energy savings to achieve further reduction in observed transmission rate. Note that the specific setpoints of these two strategies could be adjusted up or down to further tune the tradeoff. Finally, combining TCV with in-room far UV can achieve the same average and lower spread of transmission rate as the conservative baseline schedule with up to a tenfold decrease in energy cost.
Overall, these results illustrate that advanced control strategies and alternative sources of EOA can be employed to provide similar or better expected transmission rate while significantly reducing energy costs compared to constantly operating at high ventilation rates. We see that the TCV strategies deliver Pareto-optimal performance, which demonstrates that the pseudo-steady transmission model is sufficiently accurate to achieve its control objectives while remaining mathematically simple enough to integrate into existing HVAC control logic. Such TCV systems could interface with other sources of EOA, such as filtration and UV disinfection, to simultaneously control transmission rates and minimize energy consumption by prioritizing lower energy sources.
Conclusion
Our framework successfully combines data-streams from sensors with accurate physical models of infectious aerosols to predict airborne disease transmission rates in diverse indoor spaces. We have demonstrated the possibility of transmission-controlled ventilation by implementing our models in HVAC control logic. We have also provided a design framework for future buildings, which quantifies the tradeoff between different modes of clean air delivery.
Our results can also inform public health guidance, using data from real buildings. We have validated the use of the simple safety guideline, Eq. (1), to limit infection risk in different classes of indoor spaces, rather than strict occupancy limits [14]. We have shown that the underlying pseudo-steady approximation is consistent with full, dynamical simulations, and whenever small discrepancies arise, the guideline always provides a more conservative estimate of the risk. For normal occupancy in the monitored spaces, we also predict that short-range transmission via respiratory flows can be neglected (compared to the long-range airborne transmission) without imposing physical distance limits.
Controlling disease transmission must also be considered within the context of broader societal needs, such as minimizing energy usage, pollution, and carbon emissions. Our framework is able to quantify and optimize these tradeoffs, as the various sources of EOA are all incorporated, including their different energy requirements. Enabling real-time control of infection risk, prioritized against energy consumption, is a critical first step in the paradigm shift toward more healthy, energy efficient buildings [2].
Methods
Occupancy Estimation
A key input to the transmission-rate model is the (time-varying) occupancy in each space, which is then used to estimate the number of susceptible, Ns, and infectious, Ni, occupants. Our proposed approach to estimate occupancy and ventilation rates proceeds via parameter estimation applied to the following ODE model:
The state variable is the CO2 concentration
, while V is the space volume, Qb ≈ 0.6 m3/h is the occupant breathing rate,
ppm is the exhaled-breath excess CO2 concentration,
ppm is the outdoor-air CO2 concentration, and ka(t) = Qa(t)/V is the (possibly time-varying) ventilation rate. The values to be inferred are the timevarying occupancy Nt(t) and also the ventilation rate ka if it is not measured. A key benefit of the proposed approach is that it requires no prior training or actual occupancy counts, and the room volume is the only parameter that has to be specified. Exact details on how we solve the inverse problem to infer Nt(t) and ka are provided in SI 1.
Full Dynamical Model
The mass balance for infectious particles in a room results in the following partial differential equation model for the time-evolution of infectious pathogen concentration, C(r, t), per droplet size in a room of volume V and area A [14]:
where Ni are the number of infectors present in the room exhaling infectious droplets with rate P. Infectious droplets are removed through ventilation (Qa), filtrationin the recirculated airflow (pf (r)Qr), sedimentation (υs (r)A), deactivation (λυ (r)V), and the action of disinfection devices (∑d pd(r)Qd). All removal mechanisms can be expressed as rates, λa = Qa/V, λf (r) = pf (r)Qr /V, λs (r) = υs (r)A/V, and λd (r) = ∑d pd (r)Qd /V, and lumped into a single parameter the describes the supply of “equivalent outdoor air” (EOA) delivered to the space, λEOA = λa + λf + λs + λυ + λd [26].
The production rate, P, an be further expressed as P(r) = Qbnd(r) Vd(r) pm(r) cυ(r), where Qb is the breathing flow rate of the individuals; nd(r) is the number density of pathogens per volume of breath, which is known to vary with factors that include respiratory activity, time since infection, etc.; Vd(r) is the volume of the aerosol droplets; pm is a mask penetration factor which accounts for the proportion of pathogen that may be filtered out by the mask (where a value of 1 means all pathogen escapes the mask and 0 means all pathogen is filtered by the mask) [59]; and cυ(r) is the pathogen concentration in the droplets.
According to this model, the steady-state value of pathogen concentration if one infector is present is .
CO2-based Safety Guideline
The CO2-based safety guideline can be derived by considering the steady-state CO2 concentration to the dynamical model given in Equation (2). The steady-state solution is
We can rearrange Equation (4) and combine with Equation (1) to arrive at the CO2-based safety guideline
where we have assumed that Nt /Ns Ni ≈ 1.
Short-Range Transmission Risk
Estimates of short-range transmission rates can be derived from the theory of turbulent jets [14]. This analysis predicts that the concentration of infectious particles in the jets of infectors’ exhaled breath decays as 1/x where x is horizontal distance. A key deficiency of this model is that it does not account for the buoyancy of exhaled breath, which causes it to quickly rise out of the breathing zone of a potential susceptible. Thus, rather than use the turbulent jet models directly, we instead employ an empirical model derived from the experimental results of [33]. In that paper, the authors calculate a “susceptibility index” defined as ϵ := C(x)/C∞ where C(x) is the infectious particle concentration at horizontal distance x from the mouth of the infector, and C∞ is the background room concentration. For our purposes, we assume the concentration within the jet follows the model
where k is an unknown constant to be determined. Assuming C∞ = CbQb/V λEOA follows the pseudo-steady well-mixed model, we can derive the relationship
To quantify the short-range transmission rate, we use the model
in which the new parameter pshort represents the probability that a susceptible is directly within the short-range plume exhaled by each infector. More information about how we determine this parameter is provided in SI 6.
Transmission-Controlled Ventilation
Given that our primary goal is to control the transmission risk in each room, we propose using the reproductive number directly as a controlled variable, rather than using CO2 as in DCV.
To implement this control strategy, we of course first need to evaluate the current transmission rate . The pseudo-steady model gives
, in which we have removed some extra factors for brevity. The value of λEOAV can be calculated using flow measurements and filtration parameters for the BMS-provided clean air and the humidity measurements and physics-based models for the deposition and deactivation components of EOA. To estimate Ns, the CO2 generation rate
can be calculated from successive measurements of
and Qa = λaV. We could then calculate Ni := max(Nt − 1, 0), although we propose using Ni ≈ Nt to add a slight degree of robustness for small rooms. We thus arrive at the formula
which can be evaluated by the BMS.
To define the action of the controller, we thus take a transmission-rate setpoint and invert the previous formula to find the corresponding EOA setpoint
From this value, the BMS can adjust its various setpoints to deliver the required amount of EOA. In cases where the BMS can control multiple sources of EOA (e.g., ventilation, filtration via recirculation, and possibly in-zone disinfection devices), some form of prioritization would be needed, for example selecting in order of increasing energy consumption. More information on the corresponding control logic is provided in SI 7.
Data availability
All data used in this study is available at https://github.com/acoh64/MIT-JCI-IAQ-HVAC.
Code availability
Please contact M.J.R. at michael.james.risbeck{at}jci.com for code.
Author Contributions
M.Z.B. and Y.M.L. conceived the study. M.J.R. wrote the code and performed the simulations. M.J.R., J.D.D., A.E.C., Z.J., C.F., K.B., J.D., M.T., and L.B. set up the sensors and planned experiments. M.J.R. and A.E.C. prepared the manuscript. All authors discussed the results and edited the manuscript.
Competing interests
The authors declare no competing interests.
Supplementary Information
Supplementary Information is available for this work.
Acknowledgments
A.E.C. was supported by the Department of Defense (DoD) through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program. Partial funding for this study was provided by Johnson Controls, and portions of the research are subject to patents pending.
References
- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].↵
- [10].↵
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].↵
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].↵
- [51].↵
- [52].↵
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].↵
- [58].↵
- [59].↵