Application of mathematical models in combination with Monte Carlo simulation for prediction of isoflurane concentration in an operation room theater.

Applicability of two mathematical models in inhalation exposure prediction (well mixed room and near field-far field model) were validated against standard sampling method in one operation room for isoflurane. Ninety six air samples were collected from near and far field of the room and quantified by gas chromatography-flame ionization detector. Isoflurane concentration was also predicted by the models. Monte Carlo simulation was used to incorporate the role of parameters variability. The models relatively gave more conservative results than the measurements. There was no significant difference between the models and direct measurements results. There was no difference between the concentration prediction of well mixed room model and near field far field model. It suggests that the dispersion regime in room was close to well mixed situation. Direct sampling showed that the exposure in the same room for same type of operation could be up to 17 times variable which can be incorporated by Monte Carlo simulation. Mathematical models are valuable option for prediction of exposure in operation rooms. Our results also suggest that incorporating the role of parameters variability by conducting Monte Carlo simulation can enhance the strength of prediction in occupational hygiene decision making.

Introduction

Exposure and risk assessment are the main roles of occupational hygienists. In this case, inhalational exposure assessment is a vital and dominant step in human exposure assessment1). Nowadays, there are many validated and ready to use guidelines and procedures for sampling, analysis and quantification of airborne chemical hazards. These are mainly based on classic sampling and analyze methods for inhalational exposure assessment2,3,4,5). Most of these techniques need organic solvents, sampling train, and expert workforce. Therefore, by considering the problems such as shortage of budget and lack of good trained specialists, it is difficult to obtain good and even reliable results with these strategies6, 7). The range of scientific critics about the application of these methods, as an option of assessment of human inhalation exposure, is so frank that sometimes the “air sampling based methods” is called “the traditional methods”7). In addition to the above mentioned drawbacks, the results obtained by these methods did not consider the variability nature of the exposures. Some studies have shown that job exposure data suffer from between and within worker variability which have stem in analytical and environmental sources. Nicas et al. claimed that the role of environmental variability is much more than the role of analytical variability and should be taken into account in exposure assessment studies8). Exposure assessment based on single point sampling in its traditional form could not cover the variable nature of the exposure data.

During last two decades a large number of studies explored the applicability of exposure modeling techniques in the field of inhalational exposure assessment9,10,11,12). On this basis, application of mathematical models along with probabilistic approaches, such as Monte Carlo (MC) simulation, can quantify the role of the variability in the model outputs. Application of the mathematical modeling along with MC simulation can be useful in determining variability in exposure estimations13). Data obtained from these models can be used in retrospective epidemiologic studies and assessing prior exposures in the cases which no sufficient data are available14).

Among different models which have been used in this field, physico-chemical models such as zero ventilation, well − mixed room, two zone model, and eddy diffusion model are the most popular1, 10). Zero ventilation model, the simplest one, usually used as a first tier in exposure modeling. However, it usually overestimates the exposure intensity in ventilated rooms15). But, the results of this model can be used as a guide in decision making processes. Well–Mixed Room (WMR) model is located in the next tire and assumes that the contaminant mixed completely with the air and then distributed evenly in the space. WMR model underestimates the concentration in the locations near the pollution source. Because of its simplicity and easy parameter calculation, it is used frequently in rough exposure estimation. Due to these limitations, other more complicated models such as two zones near field- far field (NF-FF) and eddy diffusion model are also developed and applied in exposure estimations1, 10). In NF-FF model, the room volume is divided into the two zones, and the concentration calculations performed separately for each zone. Despite successful application of these models in evaluation of exposure intensity in different scenarios and processes, the validity of them is also under doubt and should be further evaluated in different processes. On this basis, we conducted this study to examine the applicability of mathematical exposure prediction models in the prediction of inhalational anesthetic concentrations in operation rooms.

Halogenated inhalational anesthetics are used routinely in many hospitals and veterinary clinics, among them isoflurane is one of the most recognized because of its low price and low toxicity3, 16).

Operation room medical staff could be exposed to this compound during the operations. The main source of exposure to this compound is the inhalational machine leak, spills and exhalation of the patient. In this study, we at first quantified the generation rate of isoflurane in one operation room and then applied it in the mathematical models for quantification of exposure. MC simulation was employed to involve the role of variability of exposure parameters in the modeling. Results of modeling were compared with the results obtained by direct sampling.

Subjects and Methods

All measurements and modelings were performed in an Ear, Nose & Throat (ENT) surgery room (Fig. 1). The operation room was equipped with patient bed, anesthetic delivery system and some other devices. There were four diffusers and four exhausts in the room walls and ceiling. The exhaust tube of anesthetic delivery system was placed inside the lower exhaust window. Tracer gas dilution is the most popular method to obtain the real ventilation flow and air change per hour (ACH) in a specific space. Due to some restrictions in this study we did not use this method and the real air flow in the room and thereby ACH number was calculated based on the documents of Heating, Ventilation, and Air-Conditioning (HVAC)system. Dry Kata thermometer with cooling rang of 35–38°C (CASSELLA No. M112002) was used in order to measure air velocity in the operating room.

Fig. 1.

Operation room, lay-out of ventilation system and equipments.

Air sampling was performed according to OSHA 103 standard method by charcoal tube (Anasorb 747, SKC Inc., Eighty Four, PA, USA) and pocket pump (Pocket Pump 210-1002TX, SKC) at 200 ml/min17). Samples were taken simultaneously in NF and FF in every half an hour. The samples were extracted chemically by 0.5 ml carbon disulfide (Merck, 99.5%) and analyzed by gas chromatography flame ionization detector (GC-FID) (Varian 3400 GC, Varian, Walnut Creek, CA) equipped with 10%OV-101 CWHP 80/100 (2m×1.8”ss) column. GC oven was programmed isothermally at 45°C. Injector was set at 180°C and detector temperature was also 180°C. Nitrogen with 20 ml/min flow rate was used as a carrier gas.

Crystal Ball 11.1.1.1.00 (Oracle, Redwood Shores, CA, USA) was used as a simulation tool. SPSS software package version 16 (SPSS, Inc., Chicago, IL), was used for statistical tests.

Predictive models

In this study the applicability of two mathematical concentration prediction models includes WMR model and two zones NF-FF model were evaluated10, 11, 18, 19). These models derived based on mass balance law. Some other studies also examined the applicability of these models as a predictive tool in exposure estimation10, 11, 20). But it is necessary to validate their output in more specific exposure scenarios and conditions.

The WMR model was used at first step for determination of concentration in the operation room. The concentration at time t (Ct) was calculated according to simplified equation of well mixed room WMR model with a constant emission rate (equation 1).

Equation (1)

Where: G is a generation rate (mg/min), Q: ventilation rate (m3/min), Cin: pollutant concentration in inlet flow (mg/m3), V: room volume (m3), t: time (minute) and C0 is the initial concentration at time t=0. The effect of wall loss (kl) in this study was considered negligible and therefore not evaluated.

The NF-FF model is a modified WMR model which assumes that the pollutant concentration is not same in all points of a room. It assumes that the concentration in the locations near the emission source or pollutant application source (near field) is more than the other locations (far field) in the room.

According to the NF-FF model, NF is the area which consists of the pollution source and the breathing zone of a worker. In this study, the anesthetic delivery system was considered as a pollution source in operating room. The anesthetic delivery system was near the patient’s bed. Surgeons and the operation room personnel were also beside the patient’s bed. Therefore 1.2 meter was considered as a radius of a NF hemisphere. NF volume thus was obtained by the equation 2.

Equation (2)

Where: VN: NF volume (m3) and RN: hemisphere radius (m).

Determination of the exact value of β is one of the problematic steps in the NF-FF model. In this study, β was measured according to the procedure described by Nicas10) by equation 3.

Equation (3)β = 0.5 (FSA × S)

Another important challenge in the application of mathematic models is to determine actual generation rate of pollutant (mg/min). In this study the generation rate was considered constant and the models were solved with considering constant generation rate.

Distributions were selected for desired parameters according to the published data and/or the data from the field measurements (Table 1 ). Easy fit software (Mathwave Tech.) was used for distribution fit tests.

Table 1.

Parameters distribution used in deterministic and stochastic modeling

Parameter	Symbol	Unit	Type	Distribution	Reference
Flow Rate	Q	m3/min	triangular	Max=16.1Min=13.1Mode=14.6	documentations
Generation Rate	G	mg/min	lognormal	µ=142.67SD=92.23	Sampling
Velocity	V	m/min	normal	Min=33.79Max=55	field measurement
Near Field Volume	VN	m3	–	3.62	field observation
Far Field Volume	VF	m3	–	110.48	field observation
Initial Concentration	C0	mg/m3	lognormal	µ=0.68SD=0.48	Sampling

Results

Input parameters and distributions selection

For initial concentration (C0), distribution parameters were selected by direct sampling from the operation room before beginning of daily operations. The room was supplied by fresh air from four air diffuser. Therefore, the concentration of isoflurane at the inlet flow (Cin) was considered zero. The room dimensions were measured and the volume of room was calculated. The volume of different equipments and personnel in the room also were considered in the calculation of a volume. One week observation showed that the maximum and minimum number of persons in operation room was 7 and 4 respectively. However the equipments were more or less same. By considering the standard volume of human body as 0.07 m3, calculations showed that the personnel presence in the room has negligible effect on the room volume (less than 1%). Therefore, the role of room occupation by personnel was not considered in the selection of distribution parameter.

Uniform distribution (minimum: 33.79 m/min, maximum: 55 m/min) was selected for air velocity based on 50measurements by Kata thermometer.

General ventilation rate was determined according to HAVC plots, documentations of the civil department and past measurements. The room was designed based on 8 ACH. It was in agreement with air flow measurements, therefore this value selected as an input to the models.

For quantification of generation rate, isoflurane concentration in the operation room was measured in two consecutive weeks and then G was calculated according to the equation 421).

Equation (4)G = Q × C

Where G is generation rate (mg/minute), Q is the ventilation flow rate of the room (m3/min) and C is the concentration of isoflurane (mg/m3) in the outgoing flow. Results showed that the generation rate has lognormal distribution with µ=142.67 mg/min and SD=92.23.

Modeling versus measurements

The fixed values of the model parameters were used in deterministic calculation of concentration by both models. Mean generation rate (G) obtained by field measurements was used in deterministic calculations. Figure 2 shows the time profile of isoflurane concentration obtained by application of WMR and NF-FF models in deterministic mode.

Fig. 2.

Time profile of isoflurane concentration by WMR and NF-FF models in deterministic mode.

Results showed that in WMR model (in deterministic mode), the concentration reaches steady state (Css) after 21 minutes. The Css time derived from NF-FF model by use of mean generation rate was 51 and 68 minutes for NF and FF respectively. ANOVA test showed that there is significant difference between Css derived from 3 outputs. Tukey’s post hoc multiple comparisons test was used to evaluate the difference between the concentration predictions between two models. However, the NF concentration was significantly higher than the FF values but there was no significant difference between the WMR Css and FF Css. By considering 0.5 ppmv as an exposure limit for isoflurane (in the cases of application with nitrous oxide); according to deterministic equations, the concentration will reach to this value after 8 minutes in WMR model. It was 4 minutes for NF and 6.5 minutes for FF region.

Unfortunately, deterministic approach does not incorporate the role of input parameters variability. MC simulation was used to incorporate the role of parameters variability in the concentration predictions. This approach in stochastic modeling uses repeated random sampling from input distributions to construct an output distribution.MC simulations performed according to U.S. Environmental Protection Agency (EPA)22) and Burmaster and Anderson13). Deterministic formula of WMR and NF-FF models were defined in the spreadsheet program and appropriate distributions were defined for the variables of interest. Simulation was performed with 104 iterations with MC sampling method. Table 2 also shows the TWA and Css results obtained by stochastic modeling and their proportion which is above the exposure limit.

Table 2.

TWA and Css results obtained by stochastic modeling (mg/m3)

Model	Mean			Percentage aboveexposure limit
	TWA	Css		TWA	Css
WMR	4.11	4.1		46.11	46.77
NF	4.8	4.89		53.98	55.34
FF	4.08	4.17		43.05	44.55

As can be seen from Table 2, mean of TWA exposure intensity in both models are higher than the exposure limit. Results show that NF part of NF-FF model is more conservative than the WMR model. However the statistical analysis showed that there is no significant difference between the NF-FF model and WMR results.

Air sampling

Totally 96 air samples were collected from operation room atmosphere. Forty six samples were collected in NF area and the remaining from FF area. Table 3 shows the descriptive results of the air samples in the operation room.

Table 3.

Results of air sampling in Near Field (NF) and Far Field (FF) region of operation room (n=96) (mg/m3)

Day	Location	Min	Max	Mean	SD
1	NF	0.495	0.655	0.575	0.08
	FF	0.340	0.873	0.545	0.286
2	NF	0.152	9.185	2.646	3.08
	FF	0.122	8.342	2.157	2.78
3	NF	1.405	9.33	4.36	3.26
	FF	0.66	8.65	2.36	3.52
4	NF	0.87	3.22	1.72	1.02
	FF	0.34	2.31	1.09	0.74
5	NF	0.86	8.26	3.73	3.97
	FF	0.71	6.75	3.03	3.25
6	NF	0.58	4.71	3.75	1.59
	FF	0.41	4.48	3.316	1.47
7	NF	0.48	14.26	7.063	6.14
	FF	0.35	13.36	6.48	5.69
8	NF	0.79	8.38	6.00	3.33
	FF	0.71	8.19	4.79	2.57
9	NF	0.95	19.60	10.41	9.33
	FF	0.72	18.23	9.05	8.78
10	NF	0.65	12.65	7.65	4.49
	FF	0.33	5.18	3.20	1.77

As can be seen from Table 3, the mean observed concentration of isoflurane in NF is higher than the mean value observed in FF. Theoretically there is perfect linear relationship (r2=1) between NF and FF measurements according to deterministic equations. Results of Spearman’s correlation showed there is also significant correlation between the NF and FF results of field measurements (Spearman’s rho=0.789, p<0.001). About 47.9% of NF measurements and 33.3% of FF measurements were above the isoflurane exposure limit (3.77 mg/m3).

By taking arithmetic mean from Cnf and Cff, and set as a mean concentration of room, it is possible to compare the pooled concentration in this way with the WMR results. Simple comparisons of the results obtained by stochastic modeling and direct sampling shows that the model predictions are comparable with the direct measurements results. However, Wilcoxon rank test was performed to investigate the statistical difference between WMR model and air sampling results and the results showed that there is no difference between the WMR prediction and the air sampling results (p-value=0.218).The predicted concentration by the near field far field model was compared by t-test to detect the difference between the model prediction and the standard method measurements. Results showed that there is no difference between the results predicted by NF-FF model and the results obtained by field measurements. According to Table 2, in all models more than 43% of predicted values were higher than the exposure limit. This portion is higher than the values obtained by direct measurements. But because of high variability in the field observations, it should be interpreted carefully.

The measurements in 3rd, 4th, 7th, and 10th day were performed for same type of operation. Results of direct sampling show huge variability in these measurements (about 17 times).

A simple Sensitivity analysis for the isoflurane concentration prediction was performed for the applied models by tornado chart method in Crystal ball. Testing range of variables were selected in 10% to 90% of their distributions. All stochastically defined parameters were included in sensitivity analysis. Effects of inputs were determined on predicted concentration by considering output concentration as a decision variable. The sensitivity was assessed by calculation of variables contribution to variance. Results showed that in the both models, the generation rate of pollutant has the most contribution to the concentration of the analyte (87%). However the concentration also depends on flow rate in less extent (8%).

Discussion

Despite wide application of mathematical modeling in other branches of science, it is not fully involved in occupational hygiene decision makings7). In addition, most of available models in occupational exposure assessment need to be validated against standard instruments like direct sampling. In this study, mathematical modeling in combination with stochastic simulation of the models, leads to acceptable prediction of isoflurane concentration in comparison with direct sampling results. Results of air velocity measurements in this study were different from our prior measurements in another hospital in Iran3, 23). Therefore, we propose to perform separate air velocity measurements for each specific location for modeling purposes.

With considering input parameters variability, application of single point measurement results in a specific exposure scenario is not feasible for other situations. We found that the intensity of exposure for the specific operation in the specific operation room may be hugely variable in different days. It seems that different generation rates of isoflurane because of variable amount of isoflurane administration according to patient’s body weight, age and so one lead to this variability. By considering the results of sensitivity analysis; it also revealed that the generation rate has the most important role in the variability of the results. However, application of stochastic modeling can incorporate the role of this variability in occupational hygienist’s decision making. Most of other studies performed in the field of occupational inhalation exposure modeling have also suggested that incorporation of variability in the predictive models can enhance the quality of exposure predictions20, 24, 25).

At first glance, results of NF-FF model and direct sampling suggests that the room complies with tow zone model of pollutant dispersion. However, further analysis showed that there is no significant difference between the mean NF and FF measurements. It can be interpret by redefining the model of pollutant dispersion in the room. We think high value of Q and β in combination with low volume of the room, lead to fast dispersion of pollutant in the room and change the pattern of dispersion to WMR. Plisco and Spencer26) also claimed that high value of β in comparison with Q, and small volume of a room can lead to this condition. In general with an increase in interzonal airflow (β), the mass transfer between NF and FF increased and the NF concentration became same as FF concentration. Factors such as room size and ample of ventilation (about 8 ACH) can lead to well mixing regime of dispersion. Robbins et al.20) also concluded that area concentration of benzene in far field and near field does not have significant difference; however they found that for personal samples it is significant difference between far field and near field results.

We found that the isoflurane generation rate may be up to 17 times variable during various ENT operations. This suggests that there is great variability in the exposure data of the operation room personnel even in the same class of operations. Results of sensitivity analysis suggesting that the isoflurane generation rate is the most influential parameter in concentration build up. Other factors such as Q (in the other word, the ACH) did not have strong effect on concentration. This is in agreement with Frey et al. study27).

Conclusion

Results showed that stochastic application of WMR and NF-FF model is a good screening tool in occupational hygiene decision making process. WMR model results were not different from the results obtained by the NF-FF model. The values calculated by WMR model were always lower than those measured directly by air sampling. However, considering variability in the air samples; there was no statistical difference between the measured results and those predicted by the models. In this study the exact value of anesthetic drug used in the machine was used for calculation of generation rate which yields more conservative results. Further studies should be conducted to examine the exact generation rate mechanism of the anesthetic pollutants in operation rooms. We also suggest further studies for determination of interzonal airflow and comparison of the results with available procedures.